Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
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
@alexcpn

This comment has been minimized.

Copy link
Owner Author

@alexcpn alexcpn commented Feb 20, 2016

Bounding Box
lat,long,angular_distance 13.08955 78.0335055556 0.000784584484057
delta_longitude 0.000805513923644 0.0461525481638
lat_min,lat_max 13.0445966204 13.1345033796
lon_min,lon_max 77.9873530074 78.0796581038
Distance from Point to T1 4.99859559377
Distnace from Point to T2 4.99859559377
Total time = 0.00556182861328

@alexcpn

This comment has been minimized.

Copy link
Owner Author

@alexcpn alexcpn commented Feb 20, 2016

Poles are not handled currently

@alexcpn

This comment has been minimized.

Copy link
Owner Author

@alexcpn alexcpn commented Apr 6, 2016

poles also handled

@jonathan-golorry

This comment has been minimized.

Copy link

@jonathan-golorry jonathan-golorry commented May 23, 2016

Looks like you forgot to convert lat_min to degrees.

@alexcpn

This comment has been minimized.

Copy link
Owner Author

@alexcpn alexcpn commented May 27, 2016

good catch; will change; However the (latitude_of_intersection, min_longitude) is what defines the top left corner of the bounding box

   latitude1, longitude1 = 13.04738626, 77.61946793
    bounding_dist_in_meter = .500
    min_bound, max_bound, latitude_of_intersection = get_bounding_box(bounding_dist_in_meter, latitude1,
                                                                               longitude1)
    min_latitude, min_longitude = min_bound
    max_latitude, max_longitude = max_bound
    print min_latitude
    dist_to_t1 = haversine(latitude_of_intersection, min_longitude, latitude1, longitude1)
    dist_to_t2 = haversine(latitude_of_intersection, max_longitude, latitude1, longitude1)
    self.assertAlmostEqual(bounding_dist_in_meter, dist_to_t1, places=2)
    self.assertAlmostEqual(bounding_dist_in_meter, dist_to_t2, places=2)
@geospatialology

This comment has been minimized.

Copy link

@geospatialology geospatialology commented Sep 5, 2020

Nice! Appreciate you sharing this. Thought you might want to know latitude is misspelled in the get_bounding_box function in several places stemming from the typo in the parameters for the function. Cheers

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment