Longitude is close. Latitude is way out. Progress.
This commit is contained in:
parent
d8f28ef3cc
commit
fd886daa49
2
.gitignore
vendored
2
.gitignore
vendored
|
|
@ -17,3 +17,5 @@ pom.xml.asc
|
||||||
\.idea/
|
\.idea/
|
||||||
|
|
||||||
yyy-data\.iml
|
yyy-data\.iml
|
||||||
|
|
||||||
|
resources/js/
|
||||||
|
|
|
||||||
|
|
@ -8,6 +8,8 @@
|
||||||
[org.clojure/clojure "1.8.0"]
|
[org.clojure/clojure "1.8.0"]
|
||||||
[org.clojure/data.json "0.2.6"]
|
[org.clojure/data.json "0.2.6"]
|
||||||
[org.clojure/math.numeric-tower "0.0.4"]
|
[org.clojure/math.numeric-tower "0.0.4"]
|
||||||
|
[lein-light-nrepl "0.0.4"]
|
||||||
[net.mikera/core.matrix "0.62.0"]]
|
[net.mikera/core.matrix "0.62.0"]]
|
||||||
:plugins [[lein-gorilla "0.4.0"]]
|
:plugins [[lein-gorilla "0.4.0"]]
|
||||||
|
:repl-options {:nrepl-middleware [lighttable.nrepl.handler/lighttable-ops]}
|
||||||
)
|
)
|
||||||
|
|
|
||||||
|
|
@ -30,7 +30,9 @@
|
||||||
;;;;
|
;;;;
|
||||||
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
|
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
|
||||||
|
|
||||||
;;; There are a number of bad design decisions in my implementation of this.
|
;;; Coordinate system conversion cribbed from https://www.movable-type.co.uk/scripts/latlong-os-gridref.html
|
||||||
|
|
||||||
|
;;; There are a number of bad design decisions in my re-implementation of this.
|
||||||
;;; Using protocols and records was probably a mistake.
|
;;; Using protocols and records was probably a mistake.
|
||||||
;;;
|
;;;
|
||||||
;;; The decision not only to adopt but to extend the untypable variable names
|
;;; The decision not only to adopt but to extend the untypable variable names
|
||||||
|
|
@ -95,8 +97,7 @@
|
||||||
(latitude [location])
|
(latitude [location])
|
||||||
(longitude [location])
|
(longitude [location])
|
||||||
(osgrid [location datum])
|
(osgrid [location datum])
|
||||||
(vector3d [location])
|
(vector3d [location]))
|
||||||
)
|
|
||||||
|
|
||||||
|
|
||||||
(defprotocol Transformable
|
(defprotocol Transformable
|
||||||
|
|
@ -119,26 +120,29 @@
|
||||||
(apply-transform [this transform]
|
(apply-transform [this transform]
|
||||||
(println (str "(apply-transform " this " " transform ")"))
|
(println (str "(apply-transform " this " " transform ")"))
|
||||||
(let
|
(let
|
||||||
[s1 (+ (/ (:s transform) 1e6) 1) ;; scale
|
[s1 (+' (/ (:s transform) 1e6) 1) ;; scale
|
||||||
rx (radians (/ (:rx transform) 3600)) ;; x-rotation: normalise arcseconds to radians
|
rx (radians (/ (:rx transform) 3600)) ;; x-rotation: normalise arcseconds to radians
|
||||||
ry (radians (/ (:ry transform) 3600)) ;; y-rotation: normalise arcseconds to radians
|
ry (radians (/ (:ry transform) 3600)) ;; y-rotation: normalise arcseconds to radians
|
||||||
rz (radians (/ (:rz transform) 3600)) ;; z-rotation: normalise arcseconds to radians
|
rz (radians (/ (:rz transform) 3600)) ;; z-rotation: normalise arcseconds to radians
|
||||||
x2 (-
|
;; var x2 = tx + x1*s1 - y1*rz + z1*ry;
|
||||||
(+ (:tx transform)
|
x2 (-'
|
||||||
(* (.x this) s1))
|
(+' (:tx transform)
|
||||||
(+
|
(*' (:x this) s1))
|
||||||
(* (.y this) rz)
|
(+'
|
||||||
(* (.z this) ry)))
|
(*' (:y this) rz)
|
||||||
y2 (+
|
(*' (:z this) ry)))
|
||||||
|
;; var y2 = ty + x1*rz + y1*s1 - z1*rx;
|
||||||
|
y2 (+'
|
||||||
(:ty transform)
|
(:ty transform)
|
||||||
(* (.x this) rz)
|
(*' (:x this) rz)
|
||||||
(* (.y this) s1)
|
(*' (:y this) s1)
|
||||||
(- 0 (* (.z this) rx)))
|
(-' 0 (*' (:z this) rx)))
|
||||||
z2 (+
|
;; var z2 = tz - x1*ry + y1*rx + z1*s1;
|
||||||
(.z this)
|
z2 (+'
|
||||||
(- 0 (* (.x this) ry))
|
(:tz transform)
|
||||||
(* (.y this) rx)
|
(-' 0 (*' (:x this) ry))
|
||||||
(* (.z this) s1))]
|
(*' (:y this) rx)
|
||||||
|
(*' (:z this) s1))]
|
||||||
(Vector3d. x2 y2 z2))))
|
(Vector3d. x2 y2 z2))))
|
||||||
|
|
||||||
|
|
||||||
|
|
@ -152,7 +156,7 @@
|
||||||
merge
|
merge
|
||||||
{}
|
{}
|
||||||
(map
|
(map
|
||||||
#(if (number? (t %)) (hash-map % (- 0 (t %)))) ;; (hash-map % (t %)))
|
#(if (number? (t %)) (hash-map % (-' 0 (t %)))) ;; (hash-map % (t %)))
|
||||||
(keys t))))
|
(keys t))))
|
||||||
|
|
||||||
|
|
||||||
|
|
@ -163,7 +167,7 @@
|
||||||
"Round the number `n` to `p` places after the decimal point."
|
"Round the number `n` to `p` places after the decimal point."
|
||||||
[n p]
|
[n p]
|
||||||
(let [m (reduce *' 1 (repeat p 10))]
|
(let [m (reduce *' 1 (repeat p 10))]
|
||||||
(double (/ (int (*' n m)) m))))
|
(double (/ (bigint (*' n m)) m))))
|
||||||
|
|
||||||
|
|
||||||
;; (to-fixed 1234.56789 4)
|
;; (to-fixed 1234.56789 4)
|
||||||
|
|
@ -196,10 +200,10 @@
|
||||||
location
|
location
|
||||||
(let [od (datum location)
|
(let [od (datum location)
|
||||||
nd (datums new-datum)
|
nd (datums new-datum)
|
||||||
c (vector3d location)]
|
c (geopoint-to-vector3d location)]
|
||||||
(cond
|
(cond
|
||||||
(= od nd) location
|
(= od nd) location
|
||||||
(= (:key od) :WGS84) (geopoint
|
(= (:key od) :WGS84) (vector3d-to-geopoint
|
||||||
(apply-transform c (:transform nd))
|
(apply-transform c (:transform nd))
|
||||||
(:datum location))
|
(:datum location))
|
||||||
(= (:key nd) :WGS84) (geopoint
|
(= (:key nd) :WGS84) (geopoint
|
||||||
|
|
@ -239,45 +243,54 @@
|
||||||
(osgrid [location datum] location))
|
(osgrid [location datum] location))
|
||||||
|
|
||||||
|
|
||||||
(defn vector3d-to-geopoint
|
(defn raw-vector3d-to-geopoint
|
||||||
;; My initial design decision to allow datums to be passed as keywords or as maps
|
;; My initial design decision to allow datums to be passed as keywords or as maps
|
||||||
;; was WRONG: a datum SHALL BE passed as a map
|
;; was WRONG: a datum SHALL BE passed as a key
|
||||||
[this datum]
|
([v datum]
|
||||||
|
(vector3d-to-geopoint (:x v) (:y v) (:z v) datum))
|
||||||
|
([x y z d]
|
||||||
(let
|
(let
|
||||||
[a (:a (:ellipsoid (datums datum)))
|
[a (:a (:ellipsoid (datums d)))
|
||||||
b (:b (:ellipsoid (datums datum)))
|
b (:b (:ellipsoid (datums d)))
|
||||||
f (:f (:ellipsoid (datums datum)))
|
f (:f (:ellipsoid (datums d)))
|
||||||
e² (- (* 2 f) (* f f)) ;; first eccentricity squared
|
e² (-' (*' 2 f) (*' f f)) ;; first eccentricity squared
|
||||||
ε² (/ e² (- 1 e²)) ;; second eccentricity squared
|
ε² (/ e² (-' 1 e²)) ;; second eccentricity squared
|
||||||
p (sqrt (+ (* (.x this) (.x this)) (* (.y this) (.y this))))
|
p (sqrt (+' (*' x x) (*' y y)))
|
||||||
;; distance from minor radius
|
;; distance from minor radius
|
||||||
R (sqrt (+ (* p p) (* (:z this) (:z this))))
|
R (sqrt (+' (*' p p) (*' z z)))
|
||||||
;; polar radius
|
;; polar radius
|
||||||
tanβ (* (/ (* b (:z this))(* a p)) (/ (* b (+ 1 ε²)) R))
|
tanβ (*' (/ (*' b z)(*' a p)) (/ (*' b (+' 1 ε²)) R))
|
||||||
sinβ (/ tanβ (sqrt (+ 1 (* tanβ tanβ))))
|
sinβ (/ tanβ (sqrt (+' 1 (*' tanβ tanβ))))
|
||||||
cosβ (/ sinβ tanβ)
|
cosβ (/ sinβ tanβ)
|
||||||
φ (if
|
φ (if
|
||||||
(Double/isNaN cosβ) 0
|
(Double/isNaN cosβ) 0
|
||||||
(atan2 (+ (.z this) (* ε² b sinβ sinβ sinβ))
|
(atan2 (+' z (*' ε² b sinβ sinβ sinβ))
|
||||||
(- p (* e² a cosβ cosβ cosβ))))
|
(-' p (*' e² a cosβ cosβ cosβ))))
|
||||||
λ (atan2 (:y this) (:x this))
|
λ (atan2 y x)
|
||||||
v (/ a (sqrt (- 1 (* e² (sin φ) (sin φ)))))]
|
v (/ a (sqrt (-' 1 (*' e² (sin φ) (sin φ)))))]
|
||||||
(GeoPoint. (degrees φ) (degrees λ) datum)))
|
(GeoPoint. (degrees φ) (degrees λ) d))))
|
||||||
|
|
||||||
|
|
||||||
(defn geopoint-to-osgrid
|
;; memoisation here is used more to break mutual recursion than to speed things
|
||||||
[gp datum]
|
;; up, although of course it will also speed things up a bit.
|
||||||
|
(def vector3d-to-geopoint (memoize raw-vector3d-to-geopoint))
|
||||||
|
|
||||||
|
|
||||||
|
(defn raw-geopoint-to-osgrid
|
||||||
|
([gp datum]
|
||||||
|
(geopoint-to-osgrid (:latitude gp) (:longitude gp) (:datum gp) datum))
|
||||||
|
([latitude longitude from-datum to-datum]
|
||||||
;; for bizarrely over-complicated trigonometry, look no further.
|
;; for bizarrely over-complicated trigonometry, look no further.
|
||||||
(let [point (geopoint gp datum)
|
(let [point (GeoPoint. latitude longitude to-datum)
|
||||||
φ (radians (latitude point))
|
φ (radians latitude)
|
||||||
λ (radians (longitude point))
|
λ (radians longitude)
|
||||||
a 6377563.396 ;; Airy 1830 major & minor semi-axes
|
a 6377563.396 ;; Airy 1830 major & minor semi-axes
|
||||||
b 6356256.909
|
b 6356256.909
|
||||||
F0 0.9996012717 ;; OS Grid scale factor on central meridian
|
F0 0.9996012717 ;; OS Grid scale factor on central meridian
|
||||||
φ0 (radians 49) ;; OS Grid true origin latitude
|
φ0 (radians 49) ;; OS Grid true origin latitude
|
||||||
λ0 (radians -2) ;; OS Grid true origin longitude
|
λ0 (radians -2) ;; OS Grid true origin longitude
|
||||||
Δφ (- φ φ0) ;; latitude offset from origin
|
Δφ (-' φ φ0) ;; latitude offset from origin
|
||||||
Δλ (- λ λ0) ;; longitude offset from origin
|
Δλ (-' λ λ0) ;; longitude offset from origin
|
||||||
Δλ² (reduce * 1 (repeat 2 Δλ))
|
Δλ² (reduce * 1 (repeat 2 Δλ))
|
||||||
Δλ³ (reduce * 1 (repeat 3 Δλ))
|
Δλ³ (reduce * 1 (repeat 3 Δλ))
|
||||||
Δλ⁴ (reduce * 1 (repeat 4 Δλ))
|
Δλ⁴ (reduce * 1 (repeat 4 Δλ))
|
||||||
|
|
@ -285,59 +298,63 @@
|
||||||
Δλ⁶ (reduce * 1 (repeat 6 Δλ))
|
Δλ⁶ (reduce * 1 (repeat 6 Δλ))
|
||||||
N0 -100000 ;; northing of true origin, metres
|
N0 -100000 ;; northing of true origin, metres
|
||||||
E0 400000 ;; easting of true origin, metres
|
E0 400000 ;; easting of true origin, metres
|
||||||
e² (- 1 (/ (* b b) (* a a))) ;; eccentricity squared
|
e² (-' 1 (/ (*' b b) (*' a a))) ;; eccentricity squared
|
||||||
n (/ (- a b) (+ a b))
|
n (/ (-' a b) (+' a b))
|
||||||
n² (* n n)
|
n² (*' n n)
|
||||||
n³ (* n n n)
|
n³ (*' n n n)
|
||||||
sinφ (sin φ)
|
sinφ (sin φ)
|
||||||
sin²φ (* (sin φ) (sin φ))
|
sin²φ (*' (sin φ) (sin φ))
|
||||||
cosφ (cos φ)
|
cosφ (cos φ)
|
||||||
v (* a (/ F0 (sqrt (- 1 (* e² sin²φ)))))
|
v (*' a (/ F0 (sqrt (-' 1 (*' e² sin²φ)))))
|
||||||
;; nu = transverse radius of curvature
|
;; nu = transverse radius of curvature
|
||||||
ρ (/ (* a F0 (- 1 e²)) (expt (- 1 (* e² sin²φ)) 1.5))
|
ρ (/ (*' a F0 (-' 1 e²)) (expt (-' 1 (*' e² sin²φ)) 1.5))
|
||||||
;; rho = meridional radius of curvature
|
;; rho = meridional radius of curvature
|
||||||
η2 (/ v (- ρ 1)) ;; beware confusing η2 with n²
|
η2 (/ v (-' ρ 1)) ;; beware confusing η2 with n²
|
||||||
Ma (* (+ 1 n (* (/ 5 4) n²) (* (/ 5 4) n³)) Δφ)
|
Ma (*' (+' 1 n (*' (/ 5 4) n²) (*' (/ 5 4) n³)) Δφ)
|
||||||
Mb (*
|
Mb (*
|
||||||
(+ (* 3 n) (* 3 n²) (* (/ 21 8) n³))
|
(+' (*' 3 n) (*' 3 n²) (*' (/ 21 8) n³))
|
||||||
(sin Δφ) (cos (+ φ φ0)))
|
(sin Δφ) (cos (+' φ φ0)))
|
||||||
Mc (*
|
Mc (*
|
||||||
(+ (* (/ 15 8) n²) (* (/ 15 8) n³))
|
(+' (*' (/ 15 8) n²) (*' (/ 15 8) n³))
|
||||||
(sin (* 2 Δφ)) (cos (* 2 (+ φ φ0))))
|
(sin (*' 2 Δφ)) (cos (*' 2 (+' φ φ0))))
|
||||||
Md (* (/ 35 24) n³ (sin (* 3 Δφ)) (cos (* 3 (+ φ φ0))))
|
Md (*' (/ 35 24) n³ (sin (*' 3 Δφ)) (cos (*' 3 (+' φ φ0))))
|
||||||
M (* b F0 (+ (- Ma Mb) (- Mc Md)))
|
M (*' b F0 (+' (-' Ma Mb) (-' Mc Md)))
|
||||||
;; meridional arc
|
;; meridional arc
|
||||||
cos³φ (reduce * 1 (repeat 3 cosφ))
|
cos³φ (reduce * 1 (repeat 3 cosφ))
|
||||||
cos⁵φ (reduce * 1 (repeat 5 cosφ))
|
cos⁵φ (reduce * 1 (repeat 5 cosφ))
|
||||||
tan²φ (reduce * 1 (repeat 2 (tan φ)))
|
tan²φ (reduce * 1 (repeat 2 (tan φ)))
|
||||||
tan⁴φ (reduce * 1 (repeat 4 (tan φ)))
|
tan⁴φ (reduce * 1 (repeat 4 (tan φ)))
|
||||||
I (+ M N0)
|
I (+' M N0)
|
||||||
II (* (/ v 2) sinφ cosφ)
|
II (*' (/ v 2) sinφ cosφ)
|
||||||
III (* (/ v 24) sinφ cos³φ (- 5 (+ tan²φ (* 9 η2))))
|
III (*' (/ v 24) sinφ cos³φ (-' 5 (+' tan²φ (*' 9 η2))))
|
||||||
;; var IIIA = (ν/720)*sinφ*cos5φ*(61-58*tan2φ+tan4φ);
|
;; var IIIA = (ν/720)*sinφ*cos5φ*(61-58*tan2φ+tan4φ);
|
||||||
;; TODO: CHECK JAVASCRIPT OPERATOR PRECEDENCE!
|
;; TODO: CHECK JAVASCRIPT OPERATOR PRECEDENCE!
|
||||||
IIIA (* (/ v 720) sinφ cos⁵φ (- 61 (+ (* 58 tan²φ) tan⁴φ)))
|
IIIA (*' (/ v 720) sinφ cos⁵φ (-' 61 (+' (*' 58 tan²φ) tan⁴φ)))
|
||||||
IV (* v cosφ)
|
IV (*' v cosφ)
|
||||||
V (* (/ v 6) cos³φ (- (/ v ρ) tan²φ))
|
V (*' (/ v 6) cos³φ (-' (/ v ρ) tan²φ))
|
||||||
;; var VI = (ν/120) * cos5φ * (5 - 18*tan2φ + tan4φ + 14*η2 - 58*tan2φ*η2);
|
;; var VI = (ν/120) * cos5φ * (5 - 18*tan2φ + tan4φ + 14*η2 - 58*tan2φ*η2);
|
||||||
VI (*
|
VI (*
|
||||||
(/ v 120)
|
(/ v 120)
|
||||||
cos⁵φ
|
cos⁵φ
|
||||||
(- 5
|
(-' 5
|
||||||
(+ (* 18 tan²φ)
|
(+' (*' 18 tan²φ)
|
||||||
tan⁴φ
|
tan⁴φ
|
||||||
(* 14 η2) (- 0 (* 58 tan²φ η2)))))
|
(*' 14 η2) (-' 0 (*' 58 tan²φ η2)))))
|
||||||
N (to-fixed (+ I (* II Δλ²) (+ III Δλ⁴) (* IIIA Δλ⁶)) 3)
|
N (to-fixed (+' I (*' II Δλ²) (+' III Δλ⁴) (*' IIIA Δλ⁶)) 3)
|
||||||
E (to-fixed (+ E0 (* IV Δλ) (* V Δλ³) (* VI Δλ⁵)) 3)]
|
E (to-fixed (+' E0 (*' IV Δλ) (*' V Δλ³) (*' VI Δλ⁵)) 3)]
|
||||||
(OsGrid. E N)))
|
(OsGrid. E N))))
|
||||||
|
|
||||||
|
|
||||||
(defn geopoint-to-vector3d
|
(def geopoint-to-osgrid (memoize raw-geopoint-to-osgrid))
|
||||||
[geopoint]
|
|
||||||
|
|
||||||
|
(defn raw-geopoint-to-vector3d
|
||||||
|
([gp]
|
||||||
(println (str "(geopoint-to-vector3d " geopoint ")"))
|
(println (str "(geopoint-to-vector3d " geopoint ")"))
|
||||||
(let [datum (datum geopoint)
|
(geopoint-to-vector3d (:latitude gp) (:longitude gp) (datums (:datum gp))))
|
||||||
φ (radians (latitude geopoint))
|
([latitude longitude datum]
|
||||||
λ (radians (longitude geopoint))
|
(let [φ (radians latitude)
|
||||||
|
λ (radians longitude)
|
||||||
h 0
|
h 0
|
||||||
a (:a (:ellipsoid datum))
|
a (:a (:ellipsoid datum))
|
||||||
f (:f (:ellipsoid datum))
|
f (:f (:ellipsoid datum))
|
||||||
|
|
@ -345,86 +362,91 @@
|
||||||
cosφ (cos φ)
|
cosφ (cos φ)
|
||||||
sinλ (sin λ)
|
sinλ (sin λ)
|
||||||
cosλ (cos λ)
|
cosλ (cos λ)
|
||||||
e² (- (* 2 f) (* f f))
|
e² (-' (*' 2 f) (*' f f))
|
||||||
v (/ a (sqrt (- 1 (* e² sinφ sinφ))))]
|
v (/ a (sqrt (-' 1 (*' e² sinφ sinφ))))]
|
||||||
(Vector3d.
|
(Vector3d.
|
||||||
(* (+ v h) cosφ cosλ)
|
(*' (+' v h) cosφ cosλ)
|
||||||
(* (+ v h) cosφ sinλ)
|
(*' (+' v h) cosφ sinλ)
|
||||||
(* v (+ h (- 1 e²)) sinφ))))
|
(*' v (+' h (-' 1 e²)) sinφ)))))
|
||||||
|
|
||||||
|
|
||||||
(defn osgrid-to-geopoint
|
(def geopoint-to-vector3d (memoize raw-geopoint-to-vector3d))
|
||||||
[osgrid datum]
|
|
||||||
|
|
||||||
|
(defn raw-osgrid-to-geopoint
|
||||||
|
([osgrid datum]
|
||||||
|
(osgrid-to-geopoint (:e osgrid) (:n osgrid) datum))
|
||||||
|
([E N datum]
|
||||||
(let
|
(let
|
||||||
[d datum
|
[a 6377563.396 ;; Airy 1830 major semi-axis
|
||||||
E (.e osgrid)
|
|
||||||
N (.n osgrid)
|
|
||||||
a 6377563.396 ;; Airy 1830 major semi-axis
|
|
||||||
b 6356256.909 ;; Airy 1830 minor semi-axis
|
b 6356256.909 ;; Airy 1830 minor semi-axis
|
||||||
F0 0.9996012717 ;; national grid scale factor on central meridian
|
F0 0.9996012717 ;; national grid scale factor on central meridian
|
||||||
φ0 (radians 49) ;; national grid true origin latitude
|
φ0 (radians 49) ;; national grid true origin latitude
|
||||||
λ0 (radians -2) ;; national grid true origin longitude
|
λ0 (radians -2) ;; national grid true origin longitude
|
||||||
N0 -100000 ;; northing of true origin, metres
|
N0 -100000 ;; northing of true origin, metres
|
||||||
E0 400000 ;; easting of true origin, metres
|
E0 400000 ;; easting of true origin, metres
|
||||||
e² (- 1
|
e² (-' 1
|
||||||
(/
|
(/
|
||||||
(* b b)
|
(*' b b)
|
||||||
(* a a))) ;; eccentricity squared
|
(*' a a))) ;; eccentricity squared
|
||||||
n (/ (- a b) (+ a b))
|
n (/ (-' a b) (+' a b))
|
||||||
n² (* n n)
|
n² (*' n n)
|
||||||
n³ (* n n n)
|
n³ (*' n n n)
|
||||||
[M φ] (loop [φ φ0 M 0]
|
[M φ] (loop [φ φ0 M 0]
|
||||||
(let
|
(let
|
||||||
[φ₁ (+ φ (/ (- N N0 M) (* a F0)))
|
[φ₁ (+' φ (/ (-' N N0 M) (*' a F0)))
|
||||||
Δφ (- φ φ0)
|
Δφ (-' φ φ0)
|
||||||
Ma (* (+ 1 n (* 5/4 n²) (* 5/4 n³)) Δφ)
|
Ma (*' (+' 1 n (*' 5/4 n²) (*' 5/4 n³)) Δφ)
|
||||||
Mb (* (+ (* n 3) (* n² 3) (* n³ 21/8)) (sin Δφ) (cos Δφ))
|
Mb (*' (+' (*' n 3) (*' n² 3) (*' n³ 21/8)) (sin Δφ) (cos Δφ))
|
||||||
Mc (* (+ (* n² 15/8) (* n³ 15/8)) (sin (* 2 Δφ)) (cos (* 2 Δφ)))
|
Mc (*' (+' (*' n² 15/8) (*' n³ 15/8)) (sin (*' 2 Δφ)) (cos (*' 2 Δφ)))
|
||||||
Md (* 35/24 n³ (sin (* 3 Δφ)) (cos (* 3 Δφ)))
|
Md (*' 35/24 n³ (sin (*' 3 Δφ)) (cos (*' 3 Δφ)))
|
||||||
M₁ (* b F0 (+ (- 0 Ma) Mb Mc (- 0 Md)))]
|
M₁ (*' b F0 (+' (-' 0 Ma) Mb Mc (-' 0 Md)))]
|
||||||
(if
|
(if
|
||||||
(>= (- N N0 M₁) 0.00001)
|
(>= (-' N N0 M₁) 0.00001)
|
||||||
(recur M₁ φ₁)
|
(recur M₁ φ₁)
|
||||||
[M₁ φ₁])))
|
[M₁ φ₁])))
|
||||||
sinφ (sin φ)
|
sinφ (sin φ)
|
||||||
sin²φ (* (sin φ) (sin φ))
|
sin²φ (*' (sin φ) (sin φ))
|
||||||
cosφ (cos φ)
|
cosφ (cos φ)
|
||||||
v (* a (/ F0 (sqrt (- 1 (* e² sin²φ)))))
|
v (*' a (/ F0 (sqrt (-' 1 (*' e² sin²φ)))))
|
||||||
v³ (reduce * 1 (repeat 3 v))
|
v³ (reduce * 1 (repeat 3 v))
|
||||||
v⁵ (reduce * 1 (repeat 5 v))
|
v⁵ (reduce * 1 (repeat 5 v))
|
||||||
v⁷ (reduce * 1 (repeat 7 v))
|
v⁷ (reduce * 1 (repeat 7 v))
|
||||||
;; nu = transverse radius of curvature
|
;; nu = transverse radius of curvature
|
||||||
ρ (/ (* a F0 (- 1 e²)) (expt (- 1 (* e² sin²φ)) 1.5))
|
ρ (/ (*' a F0 (-' 1 e²)) (expt (-' 1 (*' e² sin²φ)) 1.5))
|
||||||
;; rho = meridional radius of curvature
|
;; rho = meridional radius of curvature
|
||||||
η2 (/ v (- ρ 1)) ;; beware confusing η2 with n²
|
η2 (/ v (-' ρ 1)) ;; beware confusing η2 with n²
|
||||||
tanφ (tan φ)
|
tanφ (tan φ)
|
||||||
tan²φ (reduce * 1 (repeat 2 tanφ))
|
tan²φ (reduce * 1 (repeat 2 tanφ))
|
||||||
tan⁴φ (reduce * 1 (repeat 4 tanφ))
|
tan⁴φ (reduce * 1 (repeat 4 tanφ))
|
||||||
tan⁶φ (reduce * 1 (repeat 6 tanφ))
|
tan⁶φ (reduce * 1 (repeat 6 tanφ))
|
||||||
secφ (/ 1 cosφ)
|
secφ (/ 1 cosφ)
|
||||||
VII (/ tanφ (* 2 ρ v))
|
VII (/ tanφ (*' 2 ρ v))
|
||||||
;; TODO: check Javascript operator precedence!
|
;; TODO: check Javascript operator precedence!
|
||||||
VIII (* (/ tanφ (* 24 ρ v³) (+ 5 (* 3 tan²φ) η2 (- 0 (* 9 tan²φ η2)))))
|
VIII (*' (/ tanφ (*' 24 ρ v³) (+' 5 (*' 3 tan²φ) η2 (-' 0 (*' 9 tan²φ η2)))))
|
||||||
;; var VIII = tanφ/(24*ρ*ν3)*(5+3*tan2φ+η2-9*tan2φ*η2);
|
;; var VIII = tanφ/(24*ρ*ν3)*(5+3*tan2φ+η2-9*tan2φ*η2);
|
||||||
IX (* (/ tanφ (* 720 ρ v⁵)) (+ 61 (* 90 tan²φ) (* 45 tan⁴φ)))
|
IX (*' (/ tanφ (*' 720 ρ v⁵)) (+' 61 (*' 90 tan²φ) (*' 45 tan⁴φ)))
|
||||||
;; var IX = tanφ/(720*ρ*ν5)*(61+90*tan2φ+45*tan4φ);
|
;; var IX = tanφ/(720*ρ*ν5)*(61+90*tan2φ+45*tan4φ);
|
||||||
X (/ secφ v)
|
X (/ secφ v)
|
||||||
XI (* (/ secφ (* 6 v³)) (+ (/ v ρ) (* 2 tan²φ)))
|
XI (*' (/ secφ (*' 6 v³)) (+' (/ v ρ) (*' 2 tan²φ)))
|
||||||
;; var XI = secφ/(6*ν3)*(ν/ρ+2*tan2φ);
|
;; var XI = secφ/(6*ν3)*(ν/ρ+2*tan2φ);
|
||||||
XII (* (/ secφ (* 120 v⁵)) (+ 5 (* 28 tan²φ) (* 24 tan⁴φ)))
|
XII (*' (/ secφ (*' 120 v⁵)) (+' 5 (*' 28 tan²φ) (*' 24 tan⁴φ)))
|
||||||
;; var XII = secφ/(120*ν5)*(5+28*tan2φ+24*tan4φ);
|
;; var XII = secφ/(120*ν5)*(5+28*tan2φ+24*tan4φ);
|
||||||
XIIA (* (/ secφ (* 5040 v⁷)) (+ 61 (* 622 tan²φ) (* 1322 tan⁴φ) (* 720 tan⁶φ)))
|
XIIA (*' (/ secφ (*' 5040 v⁷)) (+' 61 (*' 622 tan²φ) (*' 1322 tan⁴φ) (*' 720 tan⁶φ)))
|
||||||
;; var XIIA = secφ/(5040*ν7)*(61+662*tan2φ+1320*tan4φ+720*tan6φ);
|
;; var XIIA = secφ/(5040*ν7)*(61+662*tan2φ+1320*tan4φ+720*tan6φ);
|
||||||
Δe (- E E0)
|
Δe (-' E E0)
|
||||||
Δe² (reduce * 1 (repeat 2 Δe))
|
Δe² (reduce *' 1 (repeat 2 Δe))
|
||||||
Δe³ (reduce * 1 (repeat 3 Δe))
|
Δe³ (reduce *' 1 (repeat 3 Δe))
|
||||||
Δe⁴ (reduce * 1 (repeat 4 Δe))
|
Δe⁴ (reduce *' 1 (repeat 4 Δe))
|
||||||
Δe⁵ (reduce * 1 (repeat 5 Δe))
|
Δe⁵ (reduce *' 1 (repeat 5 Δe))
|
||||||
Δe⁶ (reduce * 1 (repeat 6 Δe))
|
Δe⁶ (reduce *' 1 (repeat 6 Δe))
|
||||||
Δe⁷ (reduce * 1 (repeat 6 Δe))
|
Δe⁷ (reduce *' 1 (repeat 6 Δe))
|
||||||
φ₁ (+ φ (- 0 (+ VII Δe²)) (* VIII Δe⁴) (- 0 (* IX Δe⁷)))
|
φ₁ (+' φ (-' 0 (+' VII Δe²)) (*' VIII Δe⁴) (-' 0 (*' IX Δe⁷)))
|
||||||
λ₁ (+ λ0 (* X Δe) (- 0 (* XI Δe³)) (* IX Δe⁶) (- 0 (* XIIA Δe⁷)))]
|
λ₁ (+' λ0 (*' X Δe) (-' 0 (*' XI Δe³)) (*' IX Δe⁶) (-' 0 (*' XIIA Δe⁷)))]
|
||||||
(.geopoint (GeoPoint. (degrees φ₁) (degrees λ₁) :OSGB36) datum)))
|
(GeoPoint. (degrees φ₁) (degrees λ₁) :WGS84))))
|
||||||
|
|
||||||
|
|
||||||
|
(def osgrid-to-geopoint (memoize raw-osgrid-to-geopoint))
|
||||||
|
|
||||||
|
|
||||||
(defn post-code-data-to-addresses
|
(defn post-code-data-to-addresses
|
||||||
|
|
@ -452,42 +474,26 @@
|
||||||
|
|
||||||
;; (def home (GeoPoint. 54.8240911 -3.9170342 :WGS84))
|
;; (def home (GeoPoint. 54.8240911 -3.9170342 :WGS84))
|
||||||
|
|
||||||
;;
|
;; ;;
|
||||||
|
|
||||||
;; (:datum home)
|
;; ;; (:datum home)
|
||||||
;; (datums (:datum home))
|
;; (datums (:datum home))
|
||||||
;; (:ellipsoid (datums (:datum home)))
|
;; (:ellipsoid (datums (:datum home)))
|
||||||
;; @@
|
|
||||||
;; =>
|
|
||||||
;;; {"type":"html","content":"<span class='clj-var'>#'yyy-data.core/home</span>","value":"#'yyy-data.core/home"}
|
|
||||||
;; <=
|
|
||||||
|
|
||||||
;; @@
|
;; ;; @@
|
||||||
;;(.osgrid home :WGS84)
|
;; (osgrid home :WGS84)
|
||||||
;; @@
|
|
||||||
|
|
||||||
;; @@
|
;; (def hazelfield (OsGrid. 277656 549165))
|
||||||
;;(:datum home)
|
|
||||||
|
|
||||||
;; @@
|
;; (apply-transform
|
||||||
;; =>
|
;; (apply-transform
|
||||||
;;; {"type":"html","content":"<span class='clj-keyword'>:WGS84</span>","value":":WGS84"}
|
;; (Vector3d. 1 1 1)
|
||||||
;; <=
|
;; { :tx 89.5, :ty 93.8 :tz 123.1 :s -1.2 :rx 0.0 :ry 0.0 :rz 0.156 })
|
||||||
|
;; (inverse-transform { :tx 89.5, :ty 93.8 :tz 123.1 :s -1.2 :rx 0.0 :ry 0.0 :rz 0.156 }))
|
||||||
|
|
||||||
;; @@
|
;; (inverse-transform { :tx 89.5, :ty 93.8 :tz 123.1 :s -1.2 :rx 0.0 :ry 0.0 :rz 0.156 })
|
||||||
;; (datums (:datum home))
|
|
||||||
;; @@
|
|
||||||
;; =>
|
|
||||||
;;; {"type":"list-like","open":"<span class='clj-map'>{</span>","close":"<span class='clj-map'>}</span>","separator":", ","items":[{"type":"list-like","open":"","close":"","separator":" ","items":[{"type":"html","content":"<span class='clj-keyword'>:key</span>","value":":key"},{"type":"html","content":"<span class='clj-keyword'>:WGS84</span>","value":":WGS84"}],"value":"[:key :WGS84]"},{"type":"list-like","open":"","close":"","separator":" ","items":[{"type":"html","content":"<span class='clj-keyword'>:ellipsoid</span>","value":":ellipsoid"},{"type":"list-like","open":"<span class='clj-map'>{</span>","close":"<span class='clj-map'>}</span>","separator":", ","items":[{"type":"list-like","open":"","close":"","separator":" ","items":[{"type":"html","content":"<span class='clj-keyword'>:a</span>","value":":a"},{"type":"html","content":"<span class='clj-long'>6378137</span>","value":"6378137"}],"value":"[:a 6378137]"},{"type":"list-like","open":"","close":"","separator":" ","items":[{"type":"html","content":"<span class='clj-keyword'>:b</span>","value":":b"},{"type":"html","content":"<span class='clj-double'>6356752.314245</span>","value":"6356752.314245"}],"value":"[:b 6356752.314245]"},{"type":"list-like","open":"","close":"","separator":" ","items":[{"type":"html","content":"<span class='clj-keyword'>:f</span>","value":":f"},{"type":"html","content":"<span class='clj-double'>0.0033528106647474805</span>","value":"0.0033528106647474805"}],"value":"[:f 0.0033528106647474805]"}],"value":"{:a 6378137, :b 6356752.314245, :f 0.0033528106647474805}"}],"value":"[:ellipsoid {:a 6378137, :b 6356752.314245, :f 0.0033528106647474805}]"},{"type":"list-like","open":"","close":"","separator":" ","items":[{"type":"html","content":"<span class='clj-keyword'>:transform</span>","value":":transform"},{"type":"list-like","open":"<span class='clj-map'>{</span>","close":"<span class='clj-map'>}</span>","separator":", ","items":[{"type":"list-like","open":"","close":"","separator":" ","items":[{"type":"html","content":"<span class='clj-keyword'>:tx</span>","value":":tx"},{"type":"html","content":"<span class='clj-double'>0.0</span>","value":"0.0"}],"value":"[:tx 0.0]"},{"type":"list-like","open":"","close":"","separator":" ","items":[{"type":"html","content":"<span class='clj-keyword'>:ty</span>","value":":ty"},{"type":"html","content":"<span class='clj-double'>0.0</span>","value":"0.0"}],"value":"[:ty 0.0]"},{"type":"list-like","open":"","close":"","separator":" ","items":[{"type":"html","content":"<span class='clj-keyword'>:tz</span>","value":":tz"},{"type":"html","content":"<span class='clj-double'>0.0</span>","value":"0.0"}],"value":"[:tz 0.0]"},{"type":"list-like","open":"","close":"","separator":" ","items":[{"type":"html","content":"<span class='clj-keyword'>:s</span>","value":":s"},{"type":"html","content":"<span class='clj-double'>0.0</span>","value":"0.0"}],"value":"[:s 0.0]"},{"type":"list-like","open":"","close":"","separator":" ","items":[{"type":"html","content":"<span class='clj-keyword'>:rx</span>","value":":rx"},{"type":"html","content":"<span class='clj-double'>0.0</span>","value":"0.0"}],"value":"[:rx 0.0]"},{"type":"list-like","open":"","close":"","separator":" ","items":[{"type":"html","content":"<span class='clj-keyword'>:ry</span>","value":":ry"},{"type":"html","content":"<span class='clj-double'>0.0</span>","value":"0.0"}],"value":"[:ry 0.0]"},{"type":"list-like","open":"","close":"","separator":" ","items":[{"type":"html","content":"<span class='clj-keyword'>:rz</span>","value":":rz"},{"type":"html","content":"<span class='clj-double'>0.0</span>","value":"0.0"}],"value":"[:rz 0.0]"}],"value":"{:tx 0.0, :ty 0.0, :tz 0.0, :s 0.0, :rx 0.0, :ry 0.0, :rz 0.0}"}],"value":"[:transform {:tx 0.0, :ty 0.0, :tz 0.0, :s 0.0, :rx 0.0, :ry 0.0, :rz 0.0}]"}],"value":"{:key :WGS84, :ellipsoid {:a 6378137, :b 6356752.314245, :f 0.0033528106647474805}, :transform {:tx 0.0, :ty 0.0, :tz 0.0, :s 0.0, :rx 0.0, :ry 0.0, :rz 0.0}}"}
|
|
||||||
;; <=
|
|
||||||
|
|
||||||
;; @@
|
|
||||||
*ns*
|
|
||||||
;; @@
|
|
||||||
;; =>
|
|
||||||
;;; {"type":"html","content":"<span class='clj-unkown'>#namespace[yyy-data.core]</span>","value":"#namespace[yyy-data.core]"}
|
|
||||||
;; <=
|
|
||||||
|
|
||||||
;; @@
|
;; (geopoint (osgrid-to-geopoint hazelfield :froboz) :WGS84)
|
||||||
|
|
||||||
|
|
||||||
;; @@
|
|
||||||
|
|
|
||||||
|
|
@ -1,7 +1,65 @@
|
||||||
(ns yyy-data.core-test
|
(ns yyy-data.core-test
|
||||||
(:require [clojure.test :refer :all]
|
(:require [clojure.test :refer :all]
|
||||||
[adl-data.core :refer :all]))
|
[yyy-data.core :refer :all]))
|
||||||
|
|
||||||
|
(defn abs [n]
|
||||||
|
(if (pos? n) n (- 0 n)))
|
||||||
|
|
||||||
|
(deftest vector-tests
|
||||||
|
(testing "vector transformation round-trip."
|
||||||
|
(let [unit-vector (yyy_data.core.Vector3d. 1 1 1)
|
||||||
|
transformed (apply-transform
|
||||||
|
(apply-transform unit-vector (:transform (:ED50 datums)))
|
||||||
|
(inverse-transform (:transform (:ED50 datums))))]
|
||||||
|
(is
|
||||||
|
(< (abs (- 1 (:x transformed))) 0.001)
|
||||||
|
(str "Components should be close to 1: x: " (:x transformed)))
|
||||||
|
(is
|
||||||
|
(< (abs (- 1 (:y transformed))) 0.001)
|
||||||
|
(str "Components should be close to 1: y: " (:y transformed)))
|
||||||
|
(is
|
||||||
|
(< (abs (- 1 (:z transformed))) 0.001)
|
||||||
|
(str "Components should be close to 1: z: " (:z transformed))))))
|
||||||
|
|
||||||
|
|
||||||
|
(deftest osgrid-to-geopoint-tests
|
||||||
|
(testing "osgrid to geopoint one way"
|
||||||
|
(let
|
||||||
|
[hazelfield (yyy_data.core.OsGrid. 277656 549165)
|
||||||
|
actual (geopoint hazelfield :WGS84)]
|
||||||
|
(is
|
||||||
|
(< (abs (- (:latitude actual) 54.822218)) 0.001)
|
||||||
|
(str "Latitude should be 54.822219: " (:latitude actual)))
|
||||||
|
(is
|
||||||
|
(< (abs (- (:longitude actual) -3.906009)) 0.001)
|
||||||
|
(str "Longitude should be -3.906009: " (:longitude actual))))))
|
||||||
|
|
||||||
|
|
||||||
|
(deftest geopoint-to-osgrid-tests
|
||||||
|
(testing "geopoint to osgrid one way"
|
||||||
|
(let
|
||||||
|
[hazelfield (yyy_data.core.GeoPoint. 54.822218 -3.906009 :WGS84)
|
||||||
|
actual (osgrid hazelfield :WGS84)]
|
||||||
|
(is
|
||||||
|
(< (abs (- (:n actual) 549165)) 0.001)
|
||||||
|
(str "Northing should be 549165: " (:n actual)))
|
||||||
|
(is
|
||||||
|
(< (abs (- (:e actual) 277656)) 0.001)
|
||||||
|
(str "Easting should be 277656: " (:e actual))))))
|
||||||
|
|
||||||
|
|
||||||
|
(deftest osgrid-to-geopoint-round-trip-tests
|
||||||
|
(testing "osgrid to geopoint round trip"
|
||||||
|
(let
|
||||||
|
[hazelfield (yyy_data.core.OsGrid. 277656 549165)
|
||||||
|
latlon (geopoint hazelfield :WGS84)
|
||||||
|
hazelfield' (osgrid latlon :WGS84)]
|
||||||
|
(is
|
||||||
|
(< (abs (- (:e hazelfield) (:e hazelfield'))) 0.001)
|
||||||
|
(str "Components should be close to one another: easting: "
|
||||||
|
(:e hazelfield) ", " (:e hazelfield')))
|
||||||
|
(is
|
||||||
|
(< (abs (- (:n hazelfield) (:n hazelfield'))) 0.001)
|
||||||
|
(str "Components should be close to one another: northing: "
|
||||||
|
(:n hazelfield) ", " (:n hazelfield'))))))
|
||||||
|
|
||||||
(deftest a-test
|
|
||||||
(testing "FIXME, I fail."
|
|
||||||
(is (= 0 1))))
|
|
||||||
|
|
|
||||||
Loading…
Reference in a new issue