I have a shapefile with thousands of features (polygons). I need to randomly sample one polygon at every 1 sq. km. area if at least one polygon is available. I wrote the following code which works. However, it is very slow (took few days to complete).

import numpy as np, pandas as pd, geopandas as gpd
from shapely.geometry import Polygon

fi = 'test.shp'
gdf = gpd.read_file(fi)[['geometry']]
print (gdf.shape)        
crs = gdf.crs     

ULX, LRY, LRX, ULY = gdf.total_bounds
ts = 0.008983157 #1km
lats, lons = np.mgrid[ULY:LRY-ts:-ts, ULX:LRX+ts:ts]

ulx_s = lons[:-1,:-1].ravel()
uly_s = lats[:-1,:-1].ravel()                   
lrx_s = lons[1:,1:].ravel()
lry_s = lats[1:,1:].ravel()
tiles = [[(ulx,uly),(ulx,lry),(lrx,lry),(lrx,uly)] for ulx, uly, lrx, lry in zip(ulx_s, uly_s, lrx_s, lry_s)]

#### this is the slowest part ###########    
outputs = []          
for tile in tiles:                
    poly = Polygon(tile)
    ok = gdf[gdf.geometry.intersects(poly)]
    if ok.shape[0] >= 1:
        out = ok.sample(1)

out_df = pd.concat(outputs, axis=0)
out_gdf = gpd.GeoDataFrame(out_df, geometry='geometry')
out_gdf.crs = crs

out_gdf.to_file(fo, layer='test', driver='ESRI Shapefile')

How can I do it faster?

  • What is result of len(tiles) for your data? – Kadir Şahbaz Dec 24 '19 at 19:28
  • How many tiles do you have? – Kadir Şahbaz Dec 25 '19 at 18:37
  • You can use generators for tiles. I've editted the answer. I suggest you to review the link. – Kadir Şahbaz Dec 25 '19 at 20:21

You can use spatial index by sindex method in geopandas. I've tested on three datasets include 100, 1000, 10000 points (instead of polygons), respectively. I've used different number of tiles.

# without spatial index (for loop in the question)
outputs = []          
for tile in tiles:                
    poly = Polygon(tile)
    ok = gdf[gdf.geometry.intersects(poly)]
    if ok.shape[0] >= 1:
        out = ok.sample(1)

# with spatial index
sindex = gdf.sindex
outputs = []          
for tile in tiles:
    poly = Polygon(tile)
    candidates_index = list(sindex.intersection(poly.bounds))
    candidates = gdf.iloc[candidates_index]
    matches = candidates[candidates.intersects(poly)]
    if matches.shape[0] >= 1:
        out = matches.sample(1)

RESULTS: (times for for loop in seconds)

   Number Of        No Index   Index
Tiles   Points      (sec)     (sec)
        100         0.10       0.10
40      1000        0.50       0.12
        10000       3.50       0.23
        100         1.4        1.6
560     1000        5.6        1.6
        10000       50         1.6
        100         3.5        4.5
1420    1000        15         4.5
        10000       132        4.0
        100         8          10
3096    1000        34         10
        10000       392        10

As you can see, increase in number of points increases times extremely when not using index, but no changing when using index. When using index, in that case, number of tiles is important.

EDIT: If you have memory problem with tiles list, then you can use generator.

# Just change outer [] into (). tiles is not a list anymore, but a generator.
# convert tiles = [ ... ] to tiles = ( ... )
tiles = ([(ulx, uly), (ulx, lry), (lrx, lry), (lrx, uly)] for ulx, uly, lrx, lry in zip(ulx_s, uly_s, lrx_s, lry_s))
# remove print line. because a generator has no len function
  • Maybe, I have no experience in multiprocessing/multithreading. – Kadir Şahbaz Dec 25 '19 at 19:25

If there are (far) more polygons than grid cells, you should invert your computation, making the outer loop over the polygons. Something like:

for poly in  polygons:
  bb = boundingBox(poly)
  compute list of grid cells intersecting/containing the bb. #Note this is NOT a polygon    
      #intersection, it's a simple comparison of bounds
  for each overlapping grid cell, add poly to the list of overlapping boxes

for each cell in grid_cells:
  sample one overlapping box from list
  test to see if the polygon actually intersects the grid cell
  if false, delete the box from the list and sample again
  else add poly to your output

I also note that you are say you want 1km grid cells, but you're working in lat/lon coordinates and using a conversion of 0.008983157 degrees = 1km. That's correct for longitudes at the equator, but gets increasingly bad as you move away from the equator. You really should work in a projected coordinate system, like UTM, where the coordinates are in distance units.

  • How many polygons are in your shapefile? what is the extent (in kilometers)? – Llaves Dec 25 '19 at 18:23

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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