[bknr-cvs] r2371 - branches/bos/projects/bos/m2

hhubner at common-lisp.net hhubner at common-lisp.net
Mon Jan 21 11:10:28 UTC 2008


Author: hhubner
Date: Mon Jan 21 06:10:27 2008
New Revision: 2371

Added:
   branches/bos/projects/bos/m2/geo-utm.lisp
   branches/bos/projects/bos/m2/geometry.lisp
Modified:
   branches/bos/projects/bos/m2/bos.m2.asd
   branches/bos/projects/bos/m2/map.lisp
   branches/bos/projects/bos/m2/packages.lisp
Log:
Application-agnostic geometry functions should be put into geometry.lisp
now, point-in-polygon-p has been moved.
Add Common Lisp version of the Javascript UTM<->Lat/Lon converter, does
not yet work.



Modified: branches/bos/projects/bos/m2/bos.m2.asd
==============================================================================
--- branches/bos/projects/bos/m2/bos.m2.asd	(original)
+++ branches/bos/projects/bos/m2/bos.m2.asd	Mon Jan 21 06:10:27 2008
@@ -3,19 +3,25 @@
 (asdf:defsystem :bos.m2
   :depends-on (:bknr :bknr-modules :net.post-office :cl-mime :iconv :kmrcl :iterate :arnesi)
   :components ((:file "packages")
+	       (:file "geo-utm" :depends-on ("packages"))
+	       (:file "geometry" :depends-on ("packages"))
 	       (:file "config" :depends-on ("packages"))
 	       (:file "utils" :depends-on ("config"))
 	       (:file "news" :depends-on ("poi"))
 	       (:file "tiled-index" :depends-on ("config"))
 	       (:file "mail-generator" :depends-on ("config"))
 	       (:file "make-certificate" :depends-on ("config"))
-	       (:file "m2" :depends-on ("tiled-index" "utils" "make-certificate" "mail-generator"))
+	       (:file "m2" :depends-on ("tiled-index"
+					"utils"
+					"make-certificate"
+					"mail-generator"
+					"geo-utm"))
 	       (:file "contract-expiry" :depends-on ("m2"))
 	       (:file "allocation" :depends-on ("m2"))
-	       (:file "allocation-cache" :depends-on ("packages"))
+	       (:file "allocation-cache" :depends-on ("packages" "geometry"))
 	       (:file "poi" :depends-on ("utils" "allocation"))
 	       (:file "bitmap" :depends-on ("allocation"))
 	       (:file "import" :depends-on ("m2"))
-	       (:file "map" :depends-on ("m2" "allocation"))
+	       (:file "map" :depends-on ("m2" "allocation" "geometry"))
 	       (:file "export" :depends-on ("m2"))
 	       (:file "cert-daemon" :depends-on ("config"))))

Added: branches/bos/projects/bos/m2/geo-utm.lisp
==============================================================================
--- (empty file)
+++ branches/bos/projects/bos/m2/geo-utm.lisp	Mon Jan 21 06:10:27 2008
@@ -0,0 +1,152 @@
+(in-package :geo-utm)
+
+;; Converted from Javascript
+
+;; Origin: http:;;home.hiwaay.net/~taylorc/toolbox/geography/geoutm.html
+;; Copyright 1997-1998 by Charles L. Taylor
+
+;; Page says:
+
+;; <P>Programmers: The JavaScript source code in this document may be copied
+;; and reused without restriction.</P>
+
+(defconstant sm_a 6378137.0)
+(defconstant sm_b 6356752.314)
+(defconstant sm_EccSquared 6.69437999013e-03)
+
+(defconstant UTMScaleFactor 0.9996)
+
+(defun DegToRad (deg)
+  (* (/ deg 180.0) pi))
+
+(defun RadToDeg (rad)
+  (/ rad (* pi 180.0)))
+
+(defun ArcLengthOfMeridian (phi)
+  (let* ((n (/ (- sm_a sm_b) (+ sm_a sm_b)))
+	 (alpha (* (/ (+ sm_a sm_b) 2)
+		   (+ 1
+		      (/ (expt n 2) 4)
+		      (/ (expt n 4) 64))))
+	 (beta (+ (/ (* -3 n) 2)
+		  (/ (* 9 (expt n 3)) 16)
+		  (/ (* -3 (expt n 5)) 32)))
+	 (gamma (+ (/ (* 15 (expt n 2)) 16)
+		   (/ (* -15 (expt n 4)) 32)))
+	 (delta (+ (/ (* -35 (expt n 3)) 48)
+		   (/ (* 105 (expt n 5)) 256)))
+	 (epsilon (/ (* 315 (expt n 4)) 512)))
+    (* alpha
+       (+ phi
+	  (* beta (sin (* 2 phi)))
+	  (* gamma (sin (* 4 phi)))
+	  (* delta (sin (* 6 phi)))
+	  (* epsilon (sin (* 8 phi)))))))
+
+(defun UTMCentralMeridian (zone)
+  (DegToRad (+ -183 (* zone 6))))
+
+(defun FootpointLatitude (y)
+  (let* ((n (/ (- sm_a sm_b) (+ sm_a sm_b)))
+	 (alpha (* (/ (+ sm_a sm_b) 2)
+		   (+ 1
+		      (/ (expt n 2) 4)
+		      (/ (expt n 4) 64))))
+	 (y (/ y alpha))
+	 (beta (+ (/ (* 3 n) 2)
+		  (/ (* -27 (expt n 3)) 32)
+		  (/ (* 269 (expt n 5) 512))))
+	 (gamma (+ (/ (* 21 (expt n 2)) 16)
+		   (/ (* -55 (expt n 4)) 32)))
+	 (delta (+ (/ (* 151 (expt n 3)) 96)
+		   (/ (* -417 (expt n 5)) 128)))
+	 (epsilon (/ (* 1097 (expt n 4)) 512)))
+    (+ y
+       (* beta (sin (* 2 y)))
+       (* gamma (sin (* 4 y)))
+       (* delta (sin (* 6 y)))
+       (* epsilon (sin (* 8 y))))))
+
+(defun MapLatLonToXY (phi lambda lambda0)
+  (let* ((ep2 (/ (- (expt sm_a 2) (expt sm_b 2)) (expt sm_b 2)))
+	 (nu2 (* ep2 (expt (cos phi) 2)))
+	 (N (/ (expt sm_a 2) (* sm_b (sqrt (+ 1 nu2)))))
+	 (t- (tan phi))
+	 (t2 (* t- t-))
+	 (l (- lambda lambda0))
+	 (l3coef (+ 1 (- t2) nu2))
+	 (l4coef (+ 5 (- t2) (* 9 nu2) (* 4 nu2 nu2)))
+	 (l5coef (+ 5 (- (* 18 t2)) (* t2 t2) (* 14 nu2) (- (* 58 t2 nu2))))
+	 (l6coef (+ 61 (- (* 58 t2)) (* t2 t2) (* 270 nu2) (- (* 330 t2 nu2))))
+	 (l7coef (+ 61 (- (* 479 t2)) (* 179 t2 t2) (- (* t2 t2 t2))))
+	 (l8coef (+ 1385 (- (* 3111 t2)) (* 543 t2 t2) (- (* t2 t2 t2))))
+	 (easting (+ (* N (cos phi) l)
+		     (/ N (* 6 (expt (cos phi) 3) l3coef (expt l 3)))
+		     (/ N (* 120 (expt (cos phi) 5) l5coef (expt l 5)))
+		     (/ N (* 5040 (expt (cos phi) 7) l7coef (expt l 7)))))
+	 (northing (+ (ArcLengthOfMeridian phi)
+		      (/ t- (* 2 N (expt (cos phi) 2) (expt l 2)))
+		      (/ t- (* 24 N (expt (cos phi) 4) l4coef (expt l 4)))
+		      (/ t- (* 720 N (expt (cos phi) 6) l6coef (expt l 6)))
+		      (/ t- (* 40320 N (expt (cos phi) 8) l8coef (expt l 8))))))
+    (list easting northing)))
+
+(defun MapXYToLatLon (x y lambda0)
+  (let* ((phif (FootpointLatitude y))
+	 (ep2 (/ (- (expt sm_a 2) (expt sm_b 2)) (expt sm_b 2)))
+	 (cf (cos phif))
+	 (nuf2 (* ep2 (expt cf 2)))
+	 (Nf (/ (expt sm_a 2) (* sm_b (sqrt (+ 1 nuf2)))))
+	 (tf (tan phif))
+	 (tf2 (* tf tf))
+	 (tf4 (* tf2 tf2))
+	 (x1frac (/ 1.0 (* Nf cf)))
+	 (x2frac (/ tf (* 2 (expt Nf 2))))
+	 (x3frac (/ 1.0 (* 6 (expt Nf 3) cf)))
+	 (x4frac (/ tf (* 24 (expt Nf 4))))
+	 (x5frac (/ 1.0 (* 120 (expt Nf 5) cf)))
+	 (x6frac (/ tf (* 720 (expt Nf 6))))
+	 (x7frac (/ 1.0 (* 5040 (expt Nf 7) cf)))
+	 (x8frac (/ tf (* 40320 (expt Nf 8))))
+	 (x2poly (- -1 nuf2))
+	 (x3poly (- -1 (* 2 tf2) nuf2))
+	 (x4poly (+ 5
+		    (* 3 tf2)
+		    (* 6 nuf2)
+		    (- (* 6 tf2 nuf2))
+		    (- (* 3 nuf2 nuf2))
+		    (- (* 9 tf2 nuf2 nuf2))))
+	 (x5poly (+ 5 (* 28 tf2) (* 24 tf4) (* 6 nuf2) (* 8 tf2 nuf2)))
+	 (x6poly (+ -61 (- (* 90 tf2)) (- (* 45 tf4)) (- (* 107 nuf2)) (* 162 tf2 nuf2)))
+	 (x7poly (- 6 (* 662 tf2) (* 1320 tf4) (* 720 tf4 tf2)))
+	 (x8poly (+ 1385 (* 3633 tf2) (* 4095 tf4) (* 1575 tf4 tf2)))
+	 (latitude (+ phif
+		      (* x2frac x2poly x x)
+		      (* x4frac x4poly (expt x 4))
+		      (* x6frac x6poly (expt x 6))
+		      (* x8frac x8poly (expt x 8))))
+	 (longitude (+ lambda0
+		       (* x1frac x)
+		       (* x3frac x3poly (expt x 3))
+		       (* x5frac x5poly (expt x 5))
+		       (* x7frac x7poly (expt x 7)))))
+    (list latitude longitude)))
+
+(defun lat-lon-to-utm-x-y (lat lon &optional (zone 50))
+  ;; The Javascript version claims that the zone is calculated if not
+  ;; provided, but I could not find the code that does it.  This
+  ;; should be added.  This version of the code requires that the
+  ;; correct zone is provided.
+  (destructuring-bind (easting northing)
+      (MapLatLonToXY (DegToRad lat) (DegToRad lon) (UTMCentralMeridian zone))
+    (list (+ (* easting UTMScaleFactor) 500000)
+	  (+ (* northing UTMScaleFactor)
+	     (if (minusp northing) 10000000 0))
+	  zone)))
+
+(defun utm-x-y-to-lat-lon (x y zone southhemi)
+  (destructuring-bind (lat lon)
+      (MapXYToLatLon (/ (- x 500000) UTMScaleFactor)
+		     (/ (- y (if southhemi 10000000 0)) UTMScaleFactor)
+		     (UTMCentralMeridian zone))
+    (list (RadToDeg lat) (RadToDeg lon))))
\ No newline at end of file

Added: branches/bos/projects/bos/m2/geometry.lisp
==============================================================================
--- (empty file)
+++ branches/bos/projects/bos/m2/geometry.lisp	Mon Jan 21 06:10:27 2008
@@ -0,0 +1,19 @@
+
+(in-package :geometry)
+
+(defun point-in-polygon-p (x y polygon)
+  (let (result
+	(py y))
+    (loop with (pjx . pjy) = (aref polygon (1- (length polygon)))
+	  for (pix . piy) across polygon
+	  when (and (or (and (<= piy py) (< py pjy))
+			(and (<= pjy py) (< py piy)))
+		    (< x
+		       (+ (/ (* (- pjx pix) (- py piy))
+			     (- pjy piy))
+			  pix)))
+	  do (setf result (not result))
+	  do (setf pjx pix
+		   pjy piy))
+    result))
+

Modified: branches/bos/projects/bos/m2/map.lisp
==============================================================================
--- branches/bos/projects/bos/m2/map.lisp	(original)
+++ branches/bos/projects/bos/m2/map.lisp	Mon Jan 21 06:10:27 2008
@@ -60,22 +60,6 @@
 ;;;; kann, sie ist jedoch makrobasiert und nicht flexibel - Sie
 ;;;; erlaubt ausschließlich das Iterieren über ein Image.
 
-(defun point-in-polygon-p (x y polygon)
-  (let (result
-	(py y))
-    (loop with (pjx . pjy) = (aref polygon (1- (length polygon)))
-	  for (pix . piy) across polygon
-	  when (and (or (and (<= piy py) (< py pjy))
-			(and (<= pjy py) (< py piy)))
-		    (< x
-		       (+ (/ (* (- pjx pix) (- py piy))
-			     (- pjy piy))
-			  pix)))
-	  do (setf result (not result))
-	  do (setf pjx pix
-		   pjy piy))
-    result))
-
 (defun colorize-pixel (pixel-rgb-value color-red color-green color-blue)
   "Colorize the given PIXEL-RGB-VALUE in the COLOR given.
 PIXEL-RGB-VALUE is a raw truecolor pixel with RGB components.  COLOR

Modified: branches/bos/projects/bos/m2/packages.lisp
==============================================================================
--- branches/bos/projects/bos/m2/packages.lisp	(original)
+++ branches/bos/projects/bos/m2/packages.lisp	Mon Jan 21 06:10:27 2008
@@ -1,5 +1,14 @@
 (in-package :cl-user)
 
+(defpackage :geometry
+  (:use :cl)
+  (:export #:point-in-polygon-p))
+
+(defpackage :geo-utm
+  (:use :cl)
+  (:export #:lat-lon-to-utm-x-y
+	   #:utm-x-y-to-lat-lon))
+
 (defpackage :bos.m2.config
   (:export #:+width+
 	   #:+nw-utm-x+
@@ -27,6 +36,7 @@
   (:use :cl
 	:cl-ppcre
 	:cl-interpol
+	:geometry
 	:bknr.utils
 	:bknr.indices
 	:bknr.datastore
@@ -211,12 +221,9 @@
   (:shadowing-import-from :cl-interpol #:quote-meta-chars)
   (:export #:cert-daemon))
 
-;;; maybe there is a nicer way to do this
-;;; if you want to test this run ./build.sh at least twice !
-(intern "POINT-IN-POLYGON-P" :bos.m2) 
-
 (defpackage :bos.m2.allocation-cache
-  (:use :cl			
+  (:use :cl
+	:geometry
 	:bknr.indices
 	:bknr.datastore
 	:bknr.user       
@@ -227,10 +234,10 @@
 	:bos.m2.config	
 	:iterate
 	:arnesi)
-  (:import-from :bos.m2 bos.m2::point-in-polygon-p)
   (:export #:find-exact-match
 	   #:add-area
 	   #:free-regions-count
 	   #:free-regions-pprint
 	   #:rebuild-cache
 	   #:allocation-cache-subsystem))
+



More information about the Bknr-cvs mailing list