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.

  • Your code is 150+ lines of non-trivial code without docstrings or comments, please rather list and explain the major steps you take in plain text right here on GIS.SA. A random suggestion in case you are not using spatial indexing yet, put your land geometries into a STRTree for faster intersection checks. And maybe do it all without GeoPandas? No idea what impact on performance it might have.
    – inc42
    Commented Nov 30, 2019 at 9:01

1 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 (log2 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.
  • Thanks, that is an interesting approach. I will attempt an implementation & share the results. I was only considering a fetch length of up to 500km ( I think conversion to UTM29N would be best for my case ), the accuracy isn't too critical the intention was to develop a broad method of quantifying the exposure of a location by comparison with other locations. Commented Jun 27, 2019 at 8:39

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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