4

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?

3

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.

enter image description here

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
2

You need to reproject your geodataframes into a CRS that adequately preserves distances; UTM is a common choice, and the results would be in meters. You must determine which UTM zone to use, find the appropriate EPSG, then use the .to_crs(epsg=XXXX) method on your geodataframes to convert them.

Choosing the correct CRS in many cases is nontrivial; I have just provided a starting point that may work for you. UTM projections do not preserve distances, but over small areas, the distortion is often acceptable. The Plate Caree apparently preserves distances (epsg 32662).

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.