diff --git a/.gitignore b/.gitignore index da240e2..464a34b 100644 --- a/.gitignore +++ b/.gitignore @@ -13,3 +13,7 @@ pom.xml.asc .lein-failures .nrepl-port .cpcache/ + +\.idea/ + +yyy-data\.iml diff --git a/project.clj b/project.clj index 0af7757..dd898a1 100644 --- a/project.clj +++ b/project.clj @@ -7,4 +7,7 @@ [generateme/fastmath "1.0.1"] [org.clojure/clojure "1.8.0"] [org.clojure/data.json "0.2.6"] - [net.mikera/core.matrix "0.62.0"]]) + [org.clojure/math.numeric-tower "0.0.4"] + [net.mikera/core.matrix "0.62.0"]] + :plugins [[lein-gorilla "0.4.0"]] + ) diff --git a/src/yyy_data/core.clj b/src/yyy_data/core.clj index d09c62a..0e16326 100644 --- a/src/yyy_data/core.clj +++ b/src/yyy_data/core.clj @@ -4,11 +4,41 @@ [fastmath.core :refer [radians degrees sin cos tan sqrt atan2]] [clojure.core.matrix :as mx] [clojure.data.json :as json] + [clojure.math.numeric-tower :refer [expt]] [clojure.string :as s])) +;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; +;;;; +;;;; yyy-data.core +;;;; +;;;; This program is free software; you can redistribute it and/or +;;;; modify it under the terms of the GNU General Public License +;;;; as published by the Free Software Foundation; either version 2 +;;;; of the License, or (at your option) any later version. +;;;; +;;;; This program is distributed in the hope that it will be useful, +;;;; but WITHOUT ANY WARRANTY; without even the implied warranty of +;;;; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +;;;; GNU General Public License for more details. +;;;; +;;;; You should have received a copy of the GNU General Public License +;;;; along with this program; if not, write to the Free Software +;;;; Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, +;;;; USA. +;;;; +;;;; Copyright (C) 2018 Simon Brooke +;;;; +;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; + +;;; There are a number of bad design decisions in my implementation of this. +;;; Using protocols and records was probably a mistake. +;;; +;;; The decision not only to adopt but to extend the untypable variable names +;;; of the original was DEFINITELY a mistake. + + (declare geopoint-to-osgrid vector3d-to-geopoint geopoint-to-vector3d osgrid-to-geopoint) -;;; Coordinate system conversion cribbed from https://www.movable-type.co.uk/scripts/latlong-os-gridref.html (def ellipsoids "Ellipsoid parameters; major axis (a), minor axis (b), and flattening (f) for each ellipsoid." @@ -65,10 +95,12 @@ (defprotocol Transformable + "Something which can be transformed with a matrix transform." (apply-transform [this transform])) (defrecord Vector3d + ;; "A vector from the centre of the earth, which intercepts the surface at a point." [x y z] Location ;; datum is a bit meaningless for a Vector3d; get the default. @@ -77,11 +109,12 @@ (ellipsoid [x] (:ellipsoid (datum x))) ;; I already am vector3d; return myself (vector3d [x] x) - (geopoint [this datum] (vector3d-to-geopoint)) + (geopoint [this datum] (vector3d-to-geopoint this datum)) Transformable (apply-transform [this transform] + (println (str "(apply-transform " this " " transform ")")) (let - [s1 (+ (/ (:s transform) 1e6) 1) ;; scale + [s1 (+ (/ (:s transform) 1e6) 1) ;; scale rx (radians (/ (:rx transform) 3600)) ;; x-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 @@ -117,7 +150,9 @@ #(if (number? (t %)) (hash-map % (- 0 (t %)))) ;; (hash-map % (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 "Round the number `n` to `p` places after the decimal point." @@ -126,19 +161,13 @@ (double (/ (int (*' n m)) m)))) -(defn to-datum - ;; TODO: this obviously doesn't work but I'm trying to debug a compilation problem - [geopoint datum] - geopoint) - - ;; (to-fixed 1234.56789 4) + (defrecord GeoPoint [latitude longitude datum] - ;; "A point with an `x` co-ordinate, a `y` co-ordinate, and a datum `d`. We're - ;; agnostic as to whether `d` is passed as a keyword or a map, but it must be - ;; taken from `datums`, q.v." + ;; "A point with an `x` co-ordinate, a `y` co-ordinate, and a datum `d`. `d` + ;; must be a key taken from `datums`, q.v." Location (datum [location] (cond @@ -156,38 +185,37 @@ (:ellipsoid (datum [location]))) (vector3d [this] (geopoint-to-vector3d this)) (geopoint [location new-datum] - (let [od (datum location) - nd (if (keyword? new-datum) (datums new-datum) new-datum) - c (vector3d location)] - (cond - (= od nd) location - (= (:key od) :WGS84) (geopoint - (apply-transform c (:transform nd)) - (:datum location)) - (= (:key nd) :WGS84) (geopoint - (apply-transform - c - (inverse-transform (datums (:datum location)))) - (:datum location)) - true - (to-datum (to-datum location :WGS84) nd)))) + (println (str "(geopoint " location " " new-datum ")")) + (if + (= (:datum location) new-datum) + location + (let [od (datum location) + nd (datums new-datum) + c (vector3d location)] + (cond + (= od nd) location + (= (:key od) :WGS84) (geopoint + (apply-transform c (:transform nd)) + (:datum location)) + (= (:key nd) :WGS84) (geopoint + (apply-transform + c + (inverse-transform + (:transform od))) + (:datum location)) + true + (geopoint (geopoint location :WGS84) new-datum))))) (latitude [location] - (:latitude (to-datum location :WGS84))) + (:latitude (geopoint location :WGS84))) (longitude [location] - (:longitude (to-datum location :WGS84))) + (:longitude (geopoint location :WGS84))) (grid-x [location] (:e (osgrid location (:datum location)))) (grid-y [location] (:n (osgrid location (:datum location)))) (osgrid [location datum] - (geopoint-to-osgrid location (:datum location)))) + (geopoint-to-osgrid location datum))) -(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 [e n] @@ -207,11 +235,13 @@ (defn vector3d-to-geopoint + ;; 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 [this datum] (let - [a (:a (:ellipsoid datum)) - b (:b (:ellipsoid datum)) - f (:f (:ellipsoid datum)) + [a (:a (:ellipsoid (datums datum))) + b (:b (:ellipsoid (datums datum))) + f (:f (:ellipsoid (datums datum))) e² (- (* 2 f) (* f f)) ;; first eccentricity squared ε² (/ e² (- 1 e²)) ;; second eccentricity squared p (sqrt (+ (* (.x this) (.x this)) (* (.y this) (.y this)))) @@ -231,9 +261,9 @@ (defn geopoint-to-osgrid - [geopoint datum] + [gp datum] ;; for bizarrely over-complicated trigonometry, look no further. - (let [point (.geopoint geopoint (or datum :OSGB36)) + (let [point (geopoint gp datum) φ (radians (latitude point)) λ (radians (longitude point)) a 6377563.396 ;; Airy 1830 major & minor semi-axes @@ -259,7 +289,7 @@ cosφ (cos φ) v (* a (/ F0 (sqrt (- 1 (* e² sin²φ))))) ;; nu = transverse radius of curvature - ρ (/ (* a F0 (- 1 e²)) (pow (- 1 (* e² sin²φ)) 1.5)) + ρ (/ (* a F0 (- 1 e²)) (expt (- 1 (* e² sin²φ)) 1.5)) ;; rho = meridional radius of curvature η2 (/ v (- ρ 1)) ;; beware confusing η2 with n² Ma (* (+ 1 n (* (/ 5 4) n²) (* (/ 5 4) n³)) Δφ) @@ -299,11 +329,8 @@ (defn geopoint-to-vector3d [geopoint] - (let [datum (cond - (keyword? (:datum geopoint)) - (datums (:datum geopoint)) - (map? (:datum geopoint)) - (:datum geopoint)) + (println (str "(geopoint-to-vector3d " geopoint ")")) + (let [datum (datum geopoint) φ (radians (latitude geopoint)) λ (radians (longitude geopoint)) h 0 @@ -324,7 +351,7 @@ (defn osgrid-to-geopoint [osgrid datum] (let - [d (or datum :WGS84) + [d datum E (.e osgrid) N (.n osgrid) a 6377563.396 ;; Airy 1830 major semi-axis @@ -362,7 +389,7 @@ v⁵ (reduce * 1 (repeat 5 v)) v⁷ (reduce * 1 (repeat 7 v)) ;; nu = transverse radius of curvature - ρ (/ (* a F0 (- 1 e²)) (pow (- 1 (* e² sin²φ)) 1.5)) + ρ (/ (* a F0 (- 1 e²)) (expt (- 1 (* e² sin²φ)) 1.5)) ;; rho = meridional radius of curvature η2 (/ v (- ρ 1)) ;; beware confusing η2 with n² tanφ (tan φ) @@ -391,7 +418,7 @@ Δe⁶ (reduce * 1 (repeat 6 Δe)) Δe⁷ (reduce * 1 (repeat 6 Δ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))) @@ -416,11 +443,44 @@ :key-fn #(keyword (.toLowerCase %)))))))) -;; (def home (GeoPoint. 54.8240911 -3.9170342 :WGS84)) +;;(def home (GeoPoint. 54.8240911 -3.9170342 :WGS84)) -;; (.osgrid home :WGS84) +;; -;; (.datum home) +;; (:datum home) ;; (datums (:datum home)) ;; (:ellipsoid (datums (:datum home))) +;; @@ +;; => +;;; {"type":"html","content":"#'yyy-data.core/home","value":"#'yyy-data.core/home"} +;; <= +;; @@ +;;(.osgrid home :WGS84) +;; @@ + +;; @@ +;;(:datum home) + +;; @@ +;; => +;;; {"type":"html","content":":WGS84","value":":WGS84"} +;; <= + +;; @@ +;; (datums (:datum home)) +;; @@ +;; => +;;; {"type":"list-like","open":"{","close":"}","separator":", ","items":[{"type":"list-like","open":"","close":"","separator":" ","items":[{"type":"html","content":":key","value":":key"},{"type":"html","content":":WGS84","value":":WGS84"}],"value":"[:key :WGS84]"},{"type":"list-like","open":"","close":"","separator":" ","items":[{"type":"html","content":":ellipsoid","value":":ellipsoid"},{"type":"list-like","open":"{","close":"}","separator":", ","items":[{"type":"list-like","open":"","close":"","separator":" ","items":[{"type":"html","content":":a","value":":a"},{"type":"html","content":"6378137","value":"6378137"}],"value":"[:a 6378137]"},{"type":"list-like","open":"","close":"","separator":" ","items":[{"type":"html","content":":b","value":":b"},{"type":"html","content":"6356752.314245","value":"6356752.314245"}],"value":"[:b 6356752.314245]"},{"type":"list-like","open":"","close":"","separator":" ","items":[{"type":"html","content":":f","value":":f"},{"type":"html","content":"0.0033528106647474805","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":":transform","value":":transform"},{"type":"list-like","open":"{","close":"}","separator":", ","items":[{"type":"list-like","open":"","close":"","separator":" ","items":[{"type":"html","content":":tx","value":":tx"},{"type":"html","content":"0.0","value":"0.0"}],"value":"[:tx 0.0]"},{"type":"list-like","open":"","close":"","separator":" ","items":[{"type":"html","content":":ty","value":":ty"},{"type":"html","content":"0.0","value":"0.0"}],"value":"[:ty 0.0]"},{"type":"list-like","open":"","close":"","separator":" ","items":[{"type":"html","content":":tz","value":":tz"},{"type":"html","content":"0.0","value":"0.0"}],"value":"[:tz 0.0]"},{"type":"list-like","open":"","close":"","separator":" ","items":[{"type":"html","content":":s","value":":s"},{"type":"html","content":"0.0","value":"0.0"}],"value":"[:s 0.0]"},{"type":"list-like","open":"","close":"","separator":" ","items":[{"type":"html","content":":rx","value":":rx"},{"type":"html","content":"0.0","value":"0.0"}],"value":"[:rx 0.0]"},{"type":"list-like","open":"","close":"","separator":" ","items":[{"type":"html","content":":ry","value":":ry"},{"type":"html","content":"0.0","value":"0.0"}],"value":"[:ry 0.0]"},{"type":"list-like","open":"","close":"","separator":" ","items":[{"type":"html","content":":rz","value":":rz"},{"type":"html","content":"0.0","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":"#namespace[yyy-data.core]","value":"#namespace[yyy-data.core]"} +;; <= + +;; @@ + +;; @@