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