Commits

Marco Yuen committed 98d3502

Add an implementation to calculate a bounding box for a point.

Comments (0)

Files changed (1)

src/nextrain/utils.clj

 (ns nextrain.utils)
 
 (def origin  [40.7599573, -73.9921683])
-;(def ny_penn [40.750046, -73.992358])
-;(def pton    [40.343398, -74.659872])
+(def ny_penn [40.750046, -73.992358])
+(def pton    [40.343398, -74.659872])
 
-(def earth-radius 6371.01)
+(def ^Double earth-radius 6371.01)
 (def cos    #(Math/cos %))
+(def acos   #(Math/acos %))
 (def sin    #(Math/sin %))
 (def radian #(Math/toRadians %))
+(def degree #(Math/toDegrees %))
 (def sqrt   #(Math/sqrt %))
 (def atan2  #(Math/atan2 % %2))
 (def pow    #(Math/pow % %2))
                     (sqrt (- 1 a))))]
     (* earth-radius c)))
 
+(defn distance
+  "A simplified and less accurate way of calculating distance between two
+   points. http://www.movable-type.co.uk/scripts/latlong.html"
+  [[src-lat src-lon] [dest-lat dest-lon]]
+  (let [a (* (sin (radian src-lat))
+             (sin (radian dest-lat)))
+        b (* (cos (radian src-lat))
+             (cos (radian dest-lat))
+             (cos (- (radian dest-lon)
+                     (radian src-lon))))
+        d (acos (+ a b))]
+    (* d earth-radius)))
+
+(defn- wgs84-earth-radius
+  [lat]
+  (let [wgs84-a 6378137.0
+        wgs84-b 6356752.3
+        an (* (pow wgs84-a 2) (cos lat))
+        bn (* (pow wgs84-b 2) (sin lat))
+        ad (* wgs84-a (cos lat))
+        bd (* wgs84-b (sin lat))
+        n  (+ (pow an 2)
+              (pow bn 2))
+        d  (+ (pow ad 2)
+              (pow bd 2))]
+    (sqrt (/ n d))))
+
+(defn bounding-box
+  "A less accurate way to find a bounding box.
+  http://janmatuschek.de/LatitudeLongitudeBoundingCoordinates
+  http://stackoverflow.com/questions/238260/how-to-calculate-the-bounding-box-for-a-given-lat-lng-location"
+  [[dlat dlon] halfside]
+  (let [lat (radian dlat)
+        lon (radian dlon)
+        radius  (wgs84-earth-radius lat)
+        pradius (* radius
+                   (cos lat))
+        lat-min (- lat (/ halfside radius))
+        lat-max (+ lat (/ halfside radius))
+        lon-min (- lon (/ halfside pradius))
+        lon-max (+ lon (/ halfside pradius))]
+    [[(degree lat-min) (degree lon-min)]
+     [(degree lat-max) (degree lon-max)]]))
+