Finding shortest distance from point(s) to number of polygons at varying angles using shapely?

Question

My aim is to calculate the distance a point is from land in Python, trying to replicate the functionality of fetchR.

Given an origin and distance, I create shapely LineStrings which radiate out from this origin with 1 degree difference in angles. I then create a geopandas GeoSeries from this and run this GeoSeries of LineStrings against a GeoSeries of land_polygons first running intersects to establish if there is an intersection and then intersection to grab all of them and find the closest one.

I am really only interested in the closest intersection for that angle ( in terms of calculating fetch distances ) is there a way of constructing an intersects_first() query?

I wish to carry out this calculation for thousands of points and would be interested in ways to make the query take advantage of multiple cores

My rough implementation using Shapefile polygons and no shortening, it takes a couple of minutes in the function which uses geopandas.intersects to compute all 360 degrees on a i7 laptop.

Answer 1

Make concentric circles

Instead of making 360 different line-strings, why not make <360 different concentric circles around the point of interest, until you find the one that intersects land within your desired accuracy.

For example, start with a concentric circle around the point at a distance of $n$. If you have a shapefile in unit distance coordinates, like UTM (as opposed to lon/lat coords), you can take a point (x, y) and do:

distance = 100000 # represents 100 km if in UTM meters
point = shapely.geometry.Point(x, y)
circle = point.buffer(distance)

First, determine what you want to be the maximum range at which the function will work. Given the limits of UTM projections, this method won't be too accurate at over 1000km. So, start at 1000 km. If the actual distance is 356 km; then a search to 1 km accuracy would look like 1000->500->250->375->313->344->360->352->356->354->355. There will always be 11 circles (log₂ 1000 ~= 10 + 1 initial point).

This method will likely generate and call the intersect function on fewer shapes, though the shapes will be more complex. I'm pretty sure this will save you computationally.

Haversine everything

What is the accuracy of your shapes representing land? How many points make up those shapes? What if you just measure the distance of every shape to the point of origin using a haversine formula?

There is a lot of computational complexity building shapes and calculating intersects. haversine is just brute force math, and that tends to run really fast in numpy, in my experience. If your shapes are not that accurate (i.e. shapes with 1km accuracy instead of 20m, or something) then his approach may actually be faster.

While the speed is debatable, it is certainly

1. simpler,
2. more accurate...you get the absolute closest point, and
3. works between any two points on Earth, unaffected by projection issues.