-
-
Save kennwhite/b8a43c3b29601ad53693d71c2847c8d1 to your computer and use it in GitHub Desktop.
Testing geo bounds http://janmatuschek.de/LatitudeLongitudeBoundingCoordinates and http://gis.stackexchange.com/questions/19221/find-tangent-point-on-circle-furthest-east-or-west ; Detailed diagram for the same here http://gis.stackexchange.com/a/195148/61719
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
from __future__ import division | |
import numpy as np | |
RADIUS_OF_EARTH_IN_KM = 6371.01 | |
def haversine(lat1, lon1, lat2, lon2): | |
""" | |
Utility to calcutlate distance between two pointtodo explain regarding height | |
coverting from geodisc co-ordinate to cartestian gives errors when distances are further apart | |
""" | |
dLat = np.radians(lat2 - lat1) | |
dLon = np.radians(lon2 - lon1) | |
lat1 = np.radians(lat1) | |
lat2 = np.radians(lat2) | |
a = np.sin(dLat/2)**2 + np.cos(lat1)*np.cos(lat2)*np.sin(dLon/2)**2 | |
c = 2*np.arcsin(np.sqrt(a)) | |
return RADIUS_OF_EARTH_IN_KM * c | |
from __future__ import division | |
import numpy as np | |
def get_bounding_box(distance, latittude, longitude): | |
""" | |
Given a Geographic co-ordintate in degress, aim is to get a bounding box of a given point, | |
basically a maximum latitude, minmun latitude and max longitude and min longitude based on explanation | |
and formulae by these two sites | |
http://janmatuschek.de/LatitudeLongitudeBoundingCoordinates | |
http://gis.stackexchange.com/questions/19221/find-tangent-point-on-circle-furthest-east-or-west | |
:param distance: bounding box distance in kilometer | |
:param latittude: of the point of the center of the box | |
:param longitude: of the point of the center of the box | |
:return: top left co-ordinates, and top-right co-ordinates of the bounding box | |
""" | |
# Convert to radians | |
latittude = np.radians(latittude) | |
longitude = np.radians(longitude) | |
if not MIN_LATITUDE <= latittude <= MAX_LATITUDE: | |
raise ValueError("Invalid latitude passed in") | |
if not MIN_LONGITUDE <= longitude <= MAX_LONGITUDE: | |
raise ValueError("Invalid longitude passed in") | |
distance_from_point_km = distance | |
angular_distance = distance_from_point_km / RADIUS_OF_EARTH_IN_KM | |
# Moving along a latitude north or south, keeping the longitude fixed, is travelling along | |
# a great circle; which means basically that doing so,the Radius from center of the earth , | |
# to the circle traced by such a path is constant; So to derive the top latitude and bottom latitude | |
# without any error you just need to add or suntract the angular distance | |
lat_min = latittude - angular_distance | |
lat_max = latittude + angular_distance | |
# Handle poles here | |
if lat_min < MIN_LATITUDE or lat_max > MAX_LATITUDE: | |
lat_min = max(lat_min, MIN_LATITUDE) | |
lat_max = min(lat_max, MAX_LATITUDE) | |
lon_max = MAX_LONGITUDE | |
lon_min = MIN_LONGITUDE | |
lat_max = np.degrees(lat_max) | |
lat_min = np.degrees(lat_min) | |
lon_max = np.degrees(lon_max) | |
lon_min = np.degrees(lon_min) | |
return (lat_min, lon_min), (lat_max, lon_max), lat_max | |
# the latitude which intersects T1 and T2 in [1] and [2] is based on speherical trignomerty | |
# this is just calculated for verifying this method | |
latitude_of_intersection = np.arcsin(np.sin(latittude) / np.cos(angular_distance)) | |
latitude_of_intersection = np.degrees(latitude_of_intersection) | |
# if we draw a circle with distance around the lat and longitude, then the point where the | |
# logitude touches the circle will be at a latitude a bit apart from the original latitiude as the | |
# circle is on the sphere , the earths sphere | |
delta_longitude = np.arcsin(np.sin(angular_distance) / np.cos(latittude)) | |
lon_min = longitude - delta_longitude | |
lon_max = longitude + delta_longitude | |
lon_min = np.degrees(lon_min) | |
lat_max = np.degrees(lat_max) | |
lon_max = np.degrees(lon_max) | |
lat_min = np.degrees(lat_min) | |
return (lat_min, lon_min), (lat_max, lon_max), latitude_of_intersection |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment