Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

Find the intersection between two geographical data points

I have two pairs of lat/lon (expressed in decimal degrees) along with their radius (expressed in meters). What I am trying to achieve is to find if an intersect between these two points exits (of course, it is obvious that this doesn't hold here but the plan is to try this algorithm in many other data points). In order to check this I am using Shapely's intersects() function. My question however is how should I deal with the different units? Should I make some sort of transformation \ projection first (same units for both lat\lon and radius)?

48.180759,11.518950,19.0
47.180759,10.518950,10.0

EDIT:

I found this library here (https://pypi.python.org/pypi/utm) which seems helpfull. However, I am not 100% sure if I apply it correctly. Any ideas?

X = utm.from_latlon(38.636782, 21.414384)
A = geometry.Point(X[0], X[1]).buffer(30.777)
Y = utm.from_latlon(38.636800, 21.414488)
B = geometry.Point(Y[0], Y[1]).buffer(23.417)
A.intersects(B)

SOLUTION:

So, I finally managed to solve my problem. Here are two different implementations that both solve the same problem:

X = from_latlon(48.180759, 11.518950)
Y = from_latlon(47.180759, 10.518950)

print(latlonbuffer(48.180759, 11.518950, 19.0).intersects(latlonbuffer(47.180759, 10.518950, 19.0)))
print(latlonbuffer(48.180759, 11.518950, 100000.0).intersects(latlonbuffer(47.180759, 10.518950, 100000.0)))

X = from_latlon(48.180759, 11.518950)
Y = from_latlon(47.180759, 10.518950)

print(geometry.Point(X[0], X[1]).buffer(19.0).intersects(geometry.Point(Y[0], Y[1]).buffer(19.0)))
print(geometry.Point(X[0], X[1]).buffer(100000.0).intersects(geometry.Point(Y[0], Y[1]).buffer(100000.0)))
like image 210
user706838 Avatar asked Oct 28 '25 15:10

user706838


1 Answers

Shapely only uses the Cartesian coordinate system, so in order to make sense of metric distances, you would need to either:

  1. project the coordinates into a local projection system that uses distance units in metres, such as a UTM zone.
  2. buffer a point from (0,0), and use a dynamic azimuthal equidistant projection centered on the lat/lon point to project to geographic coords.

Here's how to do #2, using shapely.ops.transform and pyproj

import pyproj
import shapely
from shapely.geometry import Point

WGS84 = pyproj.Proj(init='epsg:4326')

def latlonbuffer(lat, lon, radius_m):
    proj4str = f'+proj=aeqd +lat_0={lat} +lon_0={lon} +x_0=0 +y_0=0 +units=m'
    AEQD = pyproj.CRS.from_proj4(proj4str)
    proj = pyproj.Transformer.from_crs(crs_from=AEQD, crs_to=WGS84).transform
    p0 = Point(0, 0).buffer(radius_m)
    return shapely.ops.transform(proj, p0)

A = latlonbuffer(48.180759, 11.518950, 19.0)
B = latlonbuffer(47.180759, 10.518950, 10.0)
print(A.intersects(B))  # False

Your two buffered points don't intersect. But these do:

A = latlonbuffer(48.180759, 11.518950, 100000.0)
B = latlonbuffer(47.180759, 10.518950, 100000.0)
print(A.intersects(B))  # True

As shown by plotting the lon/lat coords (which distorts the circles):

img

like image 145
Mike T Avatar answered Oct 31 '25 06:10

Mike T



Donate For Us

If you love us? You can donate to us via Paypal or buy me a coffee so we can maintain and grow! Thank you!