I have a single coordinate and a GeoPandas Dataframe representing US counties whose geometry is polygons that form the boundary of counties. Each of these polygons is also given through lat/lon pairs.
I need to find the distance, in miles, between my single coordinate (in lat/lon) and a county. When I use the distance function like this: point.geometry.distance(county.geometry), the distance is calculated through use of Euclidean distance. Since the points are given in terms of lat/lon, the Euclidean distance does not give the real distance between points.
How do I find the distance between a coordinate and a county using GeoPandas?
Answer
One of the packages that Geopandas depends on is Shapely. Shapely uses 2D Cartesian coordinates (x, y) and calculations are performed in 2D plane.
If you use "projected coordinate system", no problem. The distance that you get is the distance on the map (not on the spherical earth). When using "geographic coordinate system - GCS", the distance that you get will be the shortest distance in 3D space.
So you should use a formula to calculate distance on the sphere, and that is Haversine formula. But it calculates great-circle distance between two points on a sphere given their longitudes and latitudes. However, you have a point and a polygon. So you need a point on the polygon which has nearest distance to the point.
I assume that the coordinates of county polygon are on GCS, not projected. If projected, you need to transform to GCS using .to_crs()
method. Then use that code to find nearest point coordinates on the polygon. (Code Reference: Find Coordinate of Closest Point on Polygon Shapely)
from shapely.geometry import LinearRing
pol_ext = LinearRing(county.geometry.exterior.coords)
d = pol_ext.project(point.geometry)
p = pol_ext.interpolate(d)
closest_point_coords = list(p.coords)[0]
And then, use the function below to calculate spherical distance that you need. (Code Reference: Haversine Formula in Python)
from math import radians, cos, sin, asin, sqrt
def haversine(lon1, lat1, lon2, lat2):
"""
Calculate the great circle distance between two points
on the earth (specified in decimal degrees)
"""
# convert decimal degrees to radians
lon1, lat1, lon2, lat2 = map(radians, [lon1, lat1, lon2, lat2])
# haversine formula
dlon = lon2 - lon1
dlat = lat2 - lat1
a = sin(dlat/2)**2 + cos(lat1) * cos(lat2) * sin(dlon/2)**2
c = 2 * asin(sqrt(a))
r = 3956 # Radius of earth in miles. Use 6371 for kilometers
return c * r
No comments:
Post a Comment