I am trying to determine a rough distance between points in my shapefile.

The shapefile I am working with has 28 million points in the structure x y z. I am hoping to find the most common distance between each pair of points in metres.

I have come across this post How to measure distance using shapely which uses the code structure

import fiona
import shapely.geometry

# capital is a shapely geometry object

with fiona.open(path_to_data) as src:
    for feature in src:
        geom = shapely.geometry.shape(feature["geometry"])
        distance_between_pts = capital.distance(geom)

Whilst this example is comparing the distance between cities and capitals (2 different things), I'm wondering is there a way I can use fiona and shapely to give me the distance between points within 1 shapefile?

Yes I know 28 million calculations is a large magnitude. But if I could even limit the calculations to a much smaller number of points and just take a rough distance as my calculation

  • 1
    (28M*28M)/2 is a lot of distance calculations
    – Ian Turton
    Apr 28, 2022 at 15:52
  • 1
    Between each pair of point ? Do you have that powerful machine ? Wow
    – Taras
    Apr 28, 2022 at 15:53
  • @IanTurton I know ... Is there even a way to limit the calculation to the first 100 or so and I can take a rough calculation? Apr 28, 2022 at 15:57
  • I guess that you want to get that information to support the bathymetry challenge of rasterizing a depth point data set that appears in your other questions. But maybe an effective method for estimating the distances and also for getting a visual image of distribution could be to take another raster route and create a kernel density map like here desktop.arcgis.com/en/arcmap/10.3/tools/spatial-analyst-toolbox/….
    – user30184
    Apr 28, 2022 at 16:21

1 Answer 1


You can do it with itertools.combinations(). In your case:

import fiona
import shapely.geometry
import itertools

with fiona.open(path_to_data) as src:
    for feat_1, feat_2 in itertools.combinations(src, 2):
        geom_1 = shapely.geometry.shape(feat_1["geometry"])
        geom_2 = shapely.geometry.shape(feat_2["geometry"])
        distance_between_pts = geom_1.distance(geom_2)

What itertools.combinations(src, 2) does is that, given some iterator (in this case src), it returns all combinations of the features in the iterator that have length 2.

  • Fantastic! Thank you!! Quick question - is there any way to limit the number of geometries it searches for?? (to save time) Apr 28, 2022 at 19:55
  • for i, (feat_1, feat_2) in enumerate(itertools.combinations(src, 2)): ... if i > 999: break. But unless you also randomize the results may be influenced by the order they appear in the data source
    – mikewatt
    Apr 28, 2022 at 23:08
  • The way I would do it is first sample randomly a number of points from your data. And then compute the distances between them live above. This way, the results won't be influenced by the order they appear in the file
    – Jan Pisl
    Apr 29, 2022 at 7:04
  • Also I would try sample different number of points and compare the results. If for example the results you get when sampling 1k, 10k and 50k points are very similar, then it's probably a pretty good estimate of the real value you'd get if using all points
    – Jan Pisl
    Apr 29, 2022 at 7:07

Your Answer

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

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