5

I need to calculate the minimum distance (in meters) of two polygons which are defined in lat/long coordinates (EPSG:4326) using Python.

shapely geometries have distance() method which almost does what I need but as I understand first I need to reproject my polygons to some other coordinate reference system (maybe using pyproj module) to get the distance in meters.

As I understand the coordinate reference system I should use depends on which region on earth my polygons are located, but how should I find the right projection I need given the polygon coordinates.

My two polygons might be located anywhere on the earth but they will be within 10km distance from each other and 0.5m accuracy is good enough.

I have found this code on Internet, which I suppose does what I need (I'm not sure), but I understand 'EPSG:26944' is meant to be used only in California, so wouldn't work let's say in Australia, right?

from functools import partial
from shapely import ops
import pyproj

def reproject(geom, from_proj=None, to_proj=None):
    tfm = partial(pyproj.transform, pyproj.Proj(init=from_proj), pyproj.Proj(init=to_proj))
    return ops.transform(tfm, geom)

polygon1_m = reproject(shapely_polygon1, 'EPSG:4326', 'EPSG:26944')
polygon2_m = reproject(shapely_polygon2, 'EPSG:4326', 'EPSG:26944')
distance_in_meters = polygon1_m.distance(polygon2_m)
Improve this question
8
  • is your question about the reprojecting, the distance calc, or both?
    – Paul H
    Apr 17, 2019 at 16:23
  • @PaulH my goal is to calculate the minimum distance between polygons. So my question is how to do it. I just guessed for that I need to do reprojection, but I don't know how as I don't understand coordinate reference systems very well.
    – Arshak
    Apr 17, 2019 at 17:09
  • you need to read up the projected/geometric coordinate systems then. it's a large topic that is fundamental to geospatial analysis. The US has the State Plane system, developed so that engineers can do geometric calculation with acceptable levels of error over very long distances. Australia might have a similar system.
    – Paul H
    Apr 17, 2019 at 17:20
  • @PaulH Isn't there any simple way in python to calculate metric distance between two polygons without having to research large topic of coordinate systems?
    – Arshak
    Apr 17, 2019 at 18:28
  • There is. But to do so you need to be sure to pick the best coordinate system based on the geographic location and scale of your geometries. Based on the information you've provided in this question, only you are in a position to do that.
    – Paul H
    Apr 17, 2019 at 20:04

1 Answer 1

Reset to default
2

You are doing it right, for most of it. since your polygons could be anywhere on the earth you should use a coordinates system that uses meters as unit and covers the whole earth, the best fit would be the pseudo mercator 'epsg:3857'

just replace the projection you choose with this one and it should be ok

EDIT 1:

You should also try with Haversine formula on lon lat coordiantes (epsg:4326) https://pypi.org/project/haversine/

EDIT 2:

another way is to query this API and try to find the best CRS, an exemple of the query is http://www.epsg-registry.org/query.htm?name=**&geometry=bbox&north=12.2&west=10.2&south=1.2&east=2.3&validOnly=true&pagesize=10&random=66644433249

where you can get the north, west, south, east from the bounding box in lat / lon It returns many CRS depending on the area of the bounding box

Improve this answer
5
  • Hicham Zouarhi, thanks for your answer! I tried both EPSG:3857 and EPSG:26944 projections on 2 small polygons( less than 5 m diameter each) in California, for EPSG:26944 it gave a distance of 191 m between polygons, while for EPSG:3857 it was 249m. This seems too much difference to me :( Any ideas why and how to get something more accurate (which would work anywhere in the world)? Btw, I tried to measure the distance between the centroids of those polygons using great_circle distance from geopy module, it was around 193 m, so it matches with the results from EPSG:26944 projection.
    – Arshak
    Apr 19, 2019 at 22:11
  • The polygons I was trying are: polygon1 = [(-121.483071, 39.702249), (-121.483071, 39.702231), (-121.483035, 39.702231)] polygon2 = [(-121.483361, 39.700525), (-121.483361, 39.700507), (-121.483380, 39.700507)]
    – Arshak
    Apr 19, 2019 at 22:13
  • @Arshak I've edited my answer, you should look for the haversine distance either from haversine package in pypi or from numpy, it's based on geodesic coordinates so no need to reproject
    – Hicham Zouarhi
    Apr 21, 2019 at 11:56
  • Hicham Zouarhi, haversine distance from haversine package gives distance between 2 points, but how can I get the minimum distance between 2 polygons? As I mentioned earlier, my polygons are always quite small and close to each other, so any algorithm assuming flat surface should still work
    – Arshak
    Apr 23, 2019 at 12:58
  • @Arshak check this website epsg-registry.org you can get the CRS from lat/lon bounding box, I've updated the answer with an exemple of a query to its API
    – Hicham Zouarhi
    Apr 23, 2019 at 14:10

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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