Commits

Marco Yuen committed e213ed3

Change layout of index. Add implementations for calculating bounding boxes.

Comments (0)

Files changed (4)

research/geodist.py

+#!/usr/bin/python
+
+import math
+
+# http://www.platoscave.net/blog/2009/oct/5/calculate-distance-latitude-longitude-python/
+def haversine(origin, destination):
+    lat1, lon1 = origin
+    lat2, lon2 = destination
+    radius = 6371.01 # km
+
+    dlat = math.radians(lat2-lat1)
+    dlon = math.radians(lon2-lon1)
+    a = math.sin(dlat/2) * math.sin(dlat/2) + math.cos(math.radians(lat1)) \
+        * math.cos(math.radians(lat2)) * math.sin(dlon/2) * math.sin(dlon/2)
+    c = 2 * math.atan2(math.sqrt(a), math.sqrt(1-a))
+    d = radius * c
+
+    return d
+
+home    = [40.7599573, -73.9921683]
+ny_penn = [40.750046, -73.992358]
+pton    = [40.343398, -74.659872]
+
+
+dist = haversine(home, ny_penn)
+print(dist)
+
+dist = haversine(home, pton)
+print(dist)
+
+# http://www.movable-type.co.uk/scripts/latlong.html
+def lawOfCos(origin, destination):
+    lat1, lon1 = origin
+    lat2, lon2 = destination
+    rlat1, rlon1 = math.radians(lat1), math.radians(lon1)
+    rlat2, rlon2 = math.radians(lat2), math.radians(lon2)
+    radius = 6371.01 # km
+
+    a = math.sin(rlat1) * math.sin(rlat2)
+    b = math.cos(rlat1) * math.cos(rlat2) * math.cos(rlon2-rlon1)
+    d = math.acos(a + b) * radius
+
+    return d
+
+dist = lawOfCos(home, ny_penn)
+print(dist)
+
+dist = lawOfCos(home, pton)
+print(dist)
+
+# http://janmatuschek.de/LatitudeLongitudeBoundingCoordinates
+MIN_LAT = math.radians(-90.0)
+MAX_LAT = math.radians(90.0)
+MIN_LON = math.radians(-180.0)
+MAX_LON = math.radians(180.0)
+def boundingCoords(coords, distance):
+    radius = 6371.01 # km
+    lat, lon   = coords
+    rlat, rlon = math.radians(lat), math.radians(lon)
+    rdist = distance / radius
+
+    minLat = rlat - rdist
+    maxLat = rlat + rdist
+    minLon = 0.0
+    maxLon = 0.0
+
+    if minLat > MIN_LAT and maxLat < MAX_LAT:
+        deltaLon = math.asin(math.sin(rdist) / math.cos(rlat))
+        minLon   = rlon - deltaLon
+        if minLon < MIN_LON:
+            minLon += 2.0 * math.pi
+        maxLon = rlon + deltaLon
+        if maxLon > MAX_LON:
+            maxLon -= 2.0 * math.pi
+    else:
+        minLat = max(minLat, MIN_LAT)
+        maxLat = min(maxLat, MAX_LAT)
+        minLon = MIN_LON
+        maxLon = MAX_LON
+
+    return ([math.degrees(minLat), math.degrees(minLon)],
+            [math.degrees(maxLat), math.degrees(maxLon)])
+
+print("----- 5Km box from %r" % home)
+box = boundingCoords(home, 5) # 5Km
+print(box)
+
+print("----- 2Km box from %r" % home)
+box = boundingCoords(home, 2) # 5Km
+print(box)
+
+# http://stackoverflow.com/questions/238260/how-to-calculate-the-bounding-box-for-a-given-lat-lng-location
+
+# Semi-axes of WGS-84 geoidal reference
+WGS84_a = 6378137.0  # Major semiaxis [m]
+WGS84_b = 6356752.3  # Minor semiaxis [m]
+
+# Earth radius at a given latitude, according to the WGS-84 ellipsoid [m]
+def WGS84EarthRadius(lat):
+    # http://en.wikipedia.org/wiki/Earth_radius
+    An = WGS84_a*WGS84_a * math.cos(lat)
+    Bn = WGS84_b*WGS84_b * math.sin(lat)
+    Ad = WGS84_a * math.cos(lat)
+    Bd = WGS84_b * math.sin(lat)
+    return math.sqrt( (An*An + Bn*Bn)/(Ad*Ad + Bd*Bd) )
+
+# Bounding box surrounding the point at given coordinates,
+# assuming local approximation of Earth surface as a sphere
+# of radius given by WGS84
+def boundingBox(origin, halfSideInKm):
+    latitudeInDegrees, longitudeInDegrees = origin
+    lat = math.radians(latitudeInDegrees)
+    lon = math.radians(longitudeInDegrees)
+    halfSide = 1000*halfSideInKm
+
+    # Radius of Earth at given latitude
+    radius = WGS84EarthRadius(lat)
+    # Radius of the parallel at given latitude
+    pradius = radius*math.cos(lat)
+
+    latMin = lat - halfSide/radius
+    latMax = lat + halfSide/radius
+    lonMin = lon - halfSide/pradius
+    lonMax = lon + halfSide/pradius
+
+    rad2deg = math.degrees
+
+    return ([rad2deg(latMin), rad2deg(lonMin)], [rad2deg(latMax), rad2deg(lonMax)])
+
+print("----- 5Km box (Simplified) from %r" % home)
+box = boundingBox(home, 5) # 5Km
+print(box)
+
+print("----- 2Km box (Simplified) from %r" % home)
+box = boundingBox(home, 2) # 2Km
+print(box)
+

research/haversine.py

-#!/usr/bin/python
-
-import math
-
-# http://www.platoscave.net/blog/2009/oct/5/calculate-distance-latitude-longitude-python/
-def haversine(origin, destination):
-    lat1, lon1 = origin
-    lat2, lon2 = destination
-    radius = 6371 # km
-
-    dlat = math.radians(lat2-lat1)
-    dlon = math.radians(lon2-lon1)
-    a = math.sin(dlat/2) * math.sin(dlat/2) + math.cos(math.radians(lat1)) \
-        * math.cos(math.radians(lat2)) * math.sin(dlon/2) * math.sin(dlon/2)
-    c = 2 * math.atan2(math.sqrt(a), math.sqrt(1-a))
-    d = radius * c
-
-    return d
-
-home    = [40.7599573, -73.9921683]
-ny_penn = [40.750046, -73.992358]
-pton    = [40.343398, -74.659872]
-
-
-dist = haversine(home, ny_penn)
-print(dist)
-
-dist = haversine(home, pton)
-print(dist)
-

src/nextrain/pages.clj

                   "http://code.jquery.com/mobile/1.1.0/jquery.mobile-1.1.0.min.js")]
     [:body
       [:div {:data-rote "page"}
-        [:div#lng "longitude"]
-        [:div#lat "latitude"]]]))
+        [:div#lat "latitude"]
+        [:div#lng "longitude"]]]))
 

src/nextrain/utils.clj

 ;(def ny_penn [40.750046, -73.992358])
 ;(def pton    [40.343398, -74.659872])
 
-(def earth-radius 6371)
+(def earth-radius 6371.01)
 (def cos    #(Math/cos %))
 (def sin    #(Math/sin %))
 (def radian #(Math/toRadians %))