I have a set of polygons specified by geographic (WGS84) coordinates: they live on a sphere.
I have a point specified by a latitude-longitude pair.
I would like to (efficiently) find the minimum great circle distance between the point and the polygon.
My current stack includes fiona, shapely, gdal, and proj.
Similar questions on StackOverflow mostly seem to project features onto a plane and find distances there, or (disturbingly) omit mention of projections or lack thereof entirely.
This is not especially efficient since so much of the manipulation takes place in Python, rather than within a compiled library, but it does get the job done:
import shapely
import numpy as np
import math
def Pairwise(iterable):
"""
Iterate through an itertable returning adjacent pairs
:param iterable An iterable
:returns: Pairs of sequential, adjacent entries from the iterable
"""
it = iter(iterable)
a = next(it, None)
for b in it:
yield (a, b)
a = b
def LatLonToXYZ(lat,lon,radius):
"""
Convert a latitude-longitude pair to 3D-Cartesian coordinates
:param lat Latitude in degrees
:param lon Longitude in degrees
:param radius Radius in arbitrary units
:returns: x,y,z coordinates in arbitrary units
"""
lat = np.radians(lat)
lon = np.radians(lon)
x = radius * np.cos(lon) * np.cos(lat)
y = radius * np.sin(lon) * np.cos(lat)
z = radius * np.sin(lat)
return x,y,z
def XYZtoLatLon(x,y,z):
"""
Convert 3D-Cartesian coordinates to a latitude-longitude pair
:param x x-coordinate in arbitrary units
:param y y-coordinate in arbitrary units
:param z z-coordinate in arbitrary units
:returns A (lat,lon) pair in degrees
"""
radius = np.sqrt(x*x+y*y+z*z)
lat = np.degrees(np.arcsin(z/radius))
lon = np.degrees(np.arctan2(y, x))
return lat,lon
def Haversine(lat1, lon1, lat2, lon2):
"""
Calculate the Great Circle distance on Earth between two latitude-longitude
points
:param lat1 Latitude of Point 1 in degrees
:param lon1 Longtiude of Point 1 in degrees
:param lat2 Latitude of Point 2 in degrees
:param lon2 Longtiude of Point 2 in degrees
:returns Distance between the two points in kilometres
"""
Rearth = 6371
lat1 = np.radians(lat1)
lon1 = np.radians(lon1)
lat2 = np.radians(lat2)
lon2 = np.radians(lon2)
#Haversine formula
dlon = lon2 - lon1
dlat = lat2 - lat1
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 Rearth*c
#http://stackoverflow.com/a/1302268/752843
def NearestPointOnGC(alat1,alon1,alat2,alon2,plat,plon):
"""
Calculate the location of the nearest point on a Great Circle to a query point
:param lat1 Latitude of start of arc in degrees
:param lon1 Longtiude of start of arc in degrees
:param lat2 Latitude of end of arc in degrees
:param lon2 Longtiude of end of arc in degrees
:param plat Latitude of query point in degrees
:param plon Longitude of query point in degrees
:returns: A (lat,lon) pair in degrees of the closest point
"""
Rearth = 6371 #km
#Convert everything to Cartesian coordinates
a1 = np.array(LatLonToXYZ(alat1,alon1,Rearth))
a2 = np.array(LatLonToXYZ(alat2,alon2,Rearth))
p = np.array(LatLonToXYZ(plat, plon, Rearth))
G = np.cross(a1,a2) #Plane of the Great Circle containing A and B
F = np.cross(p,G) #Plane perpendicular to G that passes through query pt
T = np.cross(G,F) #Vector marking the intersection of these planes
T = Rearth*T/LA.norm(T) #Normalize to lie on the Great Circle
tlat,tlon = XYZtoLatLon(*T)
return tlat,tlon
def DistanceToGCArc(alat,alon,blat,blon,plat,plon):
"""
Calculate the distance from a query point to the nearest point on a
Great Circle Arc
:param lat1 Latitude of start of arc in degrees
:param lon1 Longtiude of start of arc in degrees
:param lat2 Latitude of end of arc in degrees
:param lon2 Longtiude of end of arc in degrees
:param plat Latitude of query point in degrees
:param plon Longitude of query point in degrees
:returns: The distance in kilometres from the query point to the great circle
arc
"""
tlat,tlon = NearestPointOnGC(alat,alon,blat,blon,plat,plon) #Nearest pt on GC
abdist = Haversine(alat,alon,blat,blon) #Length of arc
atdist = Haversine(alat,alon,tlat,tlon) #Distance arc start to nearest pt
tbdist = Haversine(tlat,tlon,blat,blon) #Distance arc end to nearest pt
#If the nearest point T on the Great Circle lies within the arc, then the
#length of the arc is approximately equal to the distance from T to each end
#of the arc, accounting for floating-point errors
PRECISION = 1e-3 #km
#We set the precision to a relatively high value because super-accuracy is not
#to needed here and a failure to catch this can lead to vast under-estimates
#of distance
if np.abs(abdist-atdist-tbdist)<PRECISION: #Nearest point was on the arc
return Haversine(tlat,tlon,plat,plon)
#Okay, the nearest point wasn't on the arc, so the nearest point is one of the
#ends points of the arc
apdist = Haversine(alat,alon,plat,plon)
bpdist = Haversine(blat,blon,plat,plon)
return min(apdist,bpdist)
def Distance3dPointTo3dPolygon(lat,lon,geom):
"""
Calculate the closest distance between a latitude-longitude query point
and a `shapely` polygon defined by latitude-longitude points using only
spherical mathematics
:param lat Latitude of query point in degrees
:param lon Longitude of query point in degrees
:param geom A `shapely` geometry whose points are in latitude-longitude space
:returns: The minimum distance in kilometres between the polygon and the
query point
"""
if geom.type == 'Polygon':
dist = math.inf
xy = features[0]['geometry'][0].exterior.xy
#Polygons are closed rings, so the first-last pair is automagically delt with
for p1, p2 in Pairwise(zip(*xy)):
dist = min(dist,DistanceToGCArc(p1[1],p1[0],p2[1],p2[0],lat,lon))
elif geom.type == 'MultiPolygon':
dist = min(*[Distance3dPointTo3dPolygon(lat,lon,part) for part in geom])
return dist
Of course, you could speed things up by only considering points and ignoring great circle arcs, as would be appropriate for polygons with suitable dense boundary specification:
def Distance3dPointTo3dPolygonQuick(lat,lon,geom):
"""
Calculate the closest distance between a polygon and a latitude-longitude
point, using only spherical considerations. Ignore edges.
:param lat Latitude of query point in degrees
:param lon Longitude of query point in degrees
:param geom A `shapely` geometry whose points are in latitude-longitude space
:returns: The minimum distance in kilometres between the polygon and the
query point
"""
if geom.type == 'Polygon':
dist = math.inf
x,y = geom.exterior.xy
#Polygons are closed rings, so the first-last pair is automagically delt with
dist = np.min(Haversine(x,y,lat,lon))
elif geom.type == 'MultiPolygon':
dist = min(*[Distance3dPointTo3dPolygonQuick(lat,lon,part) for part in geom])
return dist