Geographic Information Systems Asked by Arshak on February 2, 2021
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)
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
Answered by Hicham Zouarhi on February 2, 2021
Get help from others!
Recent Answers
Recent Questions
© 2024 TransWikia.com. All rights reserved. Sites we Love: PCI Database, UKBizDB, Menu Kuliner, Sharing RPP