pythonshapely
I have counties mapped using polygons in Shapely with latitude/longitude for coordinates. Given a complex polygon x and single point y, how would I find the shortest distance between x and y as well as the heading needed?
Using the distance function in Shapely, I can get the closest distance, but I'm unable to come up with a fast solution to get the heading. I tried determining the closest point, but this solution was unreasonably slow to run on 10 polygons.
The C++ library geographiclib has been implemented in Python. Ray tracing, or more accurately - the Geodesic Direct problem, is as easy as:
geodesic.Geodesic.WGS84.Direct(lat1, lon1, azimuth_degrees, distance_meters)
Let's test the accuracy:
from geographiclib import geodesic
from geopy.distance import vincenty
lat1, lon1 = 32.074322, 34.792081 # Azrieli Centre, Tel Aviv
distances=[]
for degree in range(360):
result=geodesic.Geodesic.WGS84.Direct(lat1, lon1, degree, 15000)
lat2, lon2 = result["lat2"], result["lon2"]
distances.append(vincenty( (lat1, lon1), (lat2, lon2)).meters)
print min(distances), max(distances)
Gives:
14999.9999876 14999.9999999
The accuracy is between 5 and 7 digits, which is great for my needs.
Links:
Instead of Haversine distance formula, you can use directly pyproj python module. For instance, for producing a 1000 m buffer (EPSG:32612) around this point (-112.171129818, 40.1731054115), EPSG:4326, next code does the job.
from shapely.geometry import Point, Polygon
import pyproj
srcProj = pyproj.Proj(init='EPSG:4326')
dstProj = pyproj.Proj(init='EPSG:32612')
pt = (-112.171129818, 40.1731054115)
x, y = pyproj.transform(srcProj, dstProj, pt[0], pt[1])
pt = Point(x,y)
buffer = pt.buffer(1000)
buffer_points = zip(buffer.exterior.coords.xy[0], buffer.exterior.coords.xy[1])
proj_buffer_points = []
for point in buffer_points:
x = point[0]
y = point[1]
x, y = pyproj.transform(dstProj, srcProj, x, y)
proj_buffer_points.append((x, y))
print Polygon(proj_buffer_points)
After running the code, at console is printed buffer (WKT format) in long/lat coordinates (EPSG:4326). With the help of QuickWKT plugin, it is visualized (EPSG:4326) at next image with QGIS.
When 'on the fly' CRS transformation to EPSG:32612 is enabled, it can be observed circular buffer shape.
Best Answer
I've no idea if this works or if it's fast enough, but I'd try this:
What makes or breaks this approach is if the buffer intersects A or not (because buffer is a circle's approximation, it might not).
If you don't get an intersection, you could increase either:
until you get an intersection.
Here's a diagram of a buffer slightly bigger then the distance to illustrate what I mean: