Still not working, but very close now.

This commit is contained in:
Simon Brooke 2018-06-25 08:29:59 +01:00
parent ef6538da9c
commit 379115a864

View file

@ -33,7 +33,7 @@
accurate to better than ±1 metre. No transformation should be assumed to be accurate to better accurate to better than ±1 metre. No transformation should be assumed to be accurate to better
than a meter; for many datums somewhat less. than a meter; for many datums somewhat less.
Yes, I know that the plural of datums is data..." Yes, I know that the plural of datum is data..."
{ {
;; transforms: t in metres, s in ppm, r in arcseconds tx ty tz s rx ry rz ;; transforms: t in metres, s in ppm, r in arcseconds tx ty tz s rx ry rz
:ED50 { :key :ED50 :ellipsoid (ellipsoids :Intl1924) :transform { :tx 89.5, :ty 93.8 :tz 123.1 :s -1.2 :rx 0.0 :ry 0.0 :rz 0.156 }} :ED50 { :key :ED50 :ellipsoid (ellipsoids :Intl1924) :transform { :tx 89.5, :ty 93.8 :tz 123.1 :s -1.2 :rx 0.0 :ry 0.0 :rz 0.156 }}
@ -101,7 +101,7 @@
(- 0 (* (.x this) ry)) (- 0 (* (.x this) ry))
(* (.y this) rx) (* (.y this) rx)
(* (.z this) s1))] (* (.z this) s1))]
(.Vector3d x2 y2 z2)))) (Vector3d. x2 y2 z2))))
(defn inverse-transform (defn inverse-transform
@ -114,10 +114,10 @@
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))))
;; (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 })
(defn to-fixed (defn to-fixed
"Round the number `n` to `p` places after the decimal point." "Round the number `n` to `p` places after the decimal point."
@ -142,16 +142,16 @@
Location Location
(datum [location] (datum [location]
(cond (cond
(keyword? (:d location)) (keyword? (:datum location))
(datums (:d location)) (datums (:datum location))
(map? (:d location)) (map? (:datum location))
(:d location))) (:datum location)))
(datum-key [location] (datum-key [location]
(cond (cond
(keyword? (:d location)) (keyword? (:datum location))
(:d location) (:datum location)
(map? (:d location)) (map? (:datum location))
(:key (:d location)))) (:key (:datum location))))
(ellipsoid [location] (ellipsoid [location]
(:ellipsoid (datum [location]))) (:ellipsoid (datum [location])))
(vector3d [this] (geopoint-to-vector3d this)) (vector3d [this] (geopoint-to-vector3d this))
@ -167,7 +167,7 @@
(= (:key nd) :WGS84) (geopoint (= (:key nd) :WGS84) (geopoint
(apply-transform (apply-transform
c c
(inverse-transform (:datum location))) (inverse-transform (datums (:datum location))))
(:datum location)) (:datum location))
true true
(to-datum (to-datum location :WGS84) nd)))) (to-datum (to-datum location :WGS84) nd))))
@ -182,7 +182,13 @@
(osgrid [location datum] (osgrid [location datum]
(geopoint-to-osgrid location (:datum location)))) (geopoint-to-osgrid location (:datum location))))
;; clojure.lang.Compiler$CompilerException: java.lang.ClassCastException: clojure.lang.PersistentList cannot be cast to clojure.lang.IPersistentVector, compiling:(/home/simon/workspace/yyy-data/src/yyy_data/core.clj:129:1) (def location (GeoPoint. 54.8240911 -3.9170342 :WGS84))
(def od (datum location))
(def nd (datums :NTF))
(def c (vector3d location))
(apply-transform c (:transform nd))
(defrecord OsGrid (defrecord OsGrid
[e n] [e n]
Location Location
@ -227,7 +233,7 @@
(defn geopoint-to-osgrid (defn geopoint-to-osgrid
[geopoint datum] [geopoint datum]
;; for bizarrely over-complicated trigonometry, look no further. ;; for bizarrely over-complicated trigonometry, look no further.
(let [point (geopoint geopoint (or datum :OSGB36)) (let [point (.geopoint geopoint (or datum :OSGB36))
φ (radians (latitude point)) φ (radians (latitude point))
λ (radians (longitude point)) λ (radians (longitude point))
a 6377563.396 ;; Airy 1830 major & minor semi-axes a 6377563.396 ;; Airy 1830 major & minor semi-axes
@ -293,11 +299,16 @@
(defn geopoint-to-vector3d (defn geopoint-to-vector3d
[geopoint] [geopoint]
(let [φ (radians (latitude geopoint)) (let [datum (cond
(keyword? (:datum geopoint))
(datums (:datum geopoint))
(map? (:datum geopoint))
(:datum geopoint))
φ (radians (latitude geopoint))
λ (radians (longitude geopoint)) λ (radians (longitude geopoint))
h 0 h 0
a (:a (:ellipsoid (:datum geopoint))) a (:a (:ellipsoid datum))
f (:f (:ellipsoid (:datum geopoint))) f (:f (:ellipsoid datum))
sinφ (sin φ) sinφ (sin φ)
cosφ (cos φ) cosφ (cos φ)
sinλ (sin λ) sinλ (sin λ)
@ -403,3 +414,13 @@
(json/read-str (json/read-str
(slurp filename) (slurp filename)
:key-fn #(keyword (.toLowerCase %)))))))) :key-fn #(keyword (.toLowerCase %))))))))
;; (def home (GeoPoint. 54.8240911 -3.9170342 :WGS84))
;; (.osgrid home :WGS84)
;; (.datum home)
;; (datums (:datum home))
;; (:ellipsoid (datums (:datum home)))