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)) +