Tell me more ×
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.
up vote 4 down vote favorite

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.

share|improve this question

1 Answer

up vote 6 down vote accepted

I've no idea if this works or if it's fast enough, but I'd try this:

  • compute distance (Z) between poly A and point B
  • build a buffer geometry (C) around point B of "radius" Z
  • compute the intersection between poly C and poly A as geometry D
  • compute centroid of geometry D
  • compute heading between D and A

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:

  • the buffer size or
  • the buffer's accuracy or
  • both

until you get an intersection.

Here's a diagram of a buffer slightly bigger then the distance to illustrate what I mean:

enter image description here

share|improve this answer
(+1) I love this idea--I've been using it for a long time. It does have problems, though. You can overcome the problem of using approximate circles by expanding them just enough to guarantee intersection. There may be more than one point of intersection, so you have to select just one. The centroid can be a poor approximation when the intersection region is long and slender (which often happens), so some post-processing of the intersection sliver to search for a reasonable closest point may be needed. Finally, all this has to be done in projected coordinates or using spherical geometry. – whuber♦ Sep 30 '11 at 14:58
The other thing I had in mind was to compute points at N, S, E, W at 10% distance away from B. Then compute distances between N, S, E, W points and polygon to figure out heading quadrant. And then keep diving the quadrant until it's enough for the required precision. But this seems less deterministic compared to the method above. – diciu Sep 30 '11 at 15:40
I don't think that second idea will work in general. It's easy to construct examples where the headings from those initial four points are all quite wrong. – whuber♦ Sep 30 '11 at 16:01
@whuber - you're right. I didn't take into account the fact that the four points might be closer to a different point on the poly in which case the whole thing fails. – diciu Sep 30 '11 at 16:12
1  
Shapely is derived from GEOS which is derived from JTS. They all share the buffer function which approximates a circle (when built around a point geometry). See here: gispython.org/shapely/docs/1.0/manual.html#buffer – diciu Oct 1 '11 at 2:46
show 1 more comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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