3

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)))
user706838
  • 5,132
  • 14
  • 54
  • 78

1 Answers1

4

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
from shapely.geometry import Point
from shapely.ops import transform
from functools import partial

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

def latlonbuffer(lat, lon, radius_m):
    proj4str = '+proj=aeqd +lat_0=%s +lon_0=%s +x_0=0 +y_0=0' % (lat, lon)
    AEQD = pyproj.Proj(proj4str)
    project = partial(pyproj.transform, AEQD, WGS84)
    return transform(project, Point(0, 0).buffer(radius_m))

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

Mike T
  • 41,085
  • 18
  • 152
  • 203
  • cool. haven't tested yet, however, I had a look on how to implement your #1 suggestion. please have a look on my EDIT and let me know what you think. :) – user706838 Dec 12 '14 at 11:50
  • @user706838 cool, it looks like that would work too, but you need to make sure you are comparing two geometries from the same UTM zone. Places like Calgary this would not always work, since there is a divide between UTM zones right through the city. – Mike T Dec 13 '14 at 00:06
  • hmmm so what you are saying is that YOUR solution is more robust than the other one and hence I have to go with yours? Not problem with that, but please clarify :) – user706838 Dec 22 '14 at 09:11
  • It depends. If you have a small region where the points always coincide in the same UTM zone, then your technique will work well. (It would be nice if the `utm` module could force a zone). They should give the same answer. If you need really high precision distances, my first technique should be good for this. UTM zones tend to warp near the edges. – Mike T Dec 22 '14 at 19:02
  • +1 I am wondering why it is transform(project, Point(0, 0).buffer(radius_m)) instead of transform(project, Point(lon,lat).buffer(radius_m)) – Eduardo Jan 30 '17 at 17:51