Maximizing Code Performance for Shapely?

Question

I've written this code in order to calculate the percent cover of grassland within a national park, within 20km^2 of points of recorded occurrences for a species. It is designed to only consider the areas within the park, and not those outside. It works, except for one huge issue... it's *incredibly* slow. I think I've tracked it down to the `avail_lc_type = pt_buffer.intersection(gl).area`. The gl polygon has > 24,000 polygons in it. It is a raster to polygon transformation (it has been dissolved), so it has a lot of little polygons within. Right now it's not a huge problem since I'm running it on ~300 points (still takes > 1 hr), but I plan to run it on a few million later, so I need to improve it. Any ideas? <!-- language: lang-py --> import numpy as np import pprint import shapely from shapely.geometry import* import fiona from fiona import collection import math traps = fiona.open('some_points.shp', 'r') #point file: focal points around which statistics are being derived study_area = fiona.open('available_areas.shp', 'r') #polygon file: represents area available for analysis for i in study_area: #for every record in 'study_area' sa = shape(i['geometry']) #make a variable called 'sa' that is a polygon grassland = fiona.open('land_cover_type_of_interest.shp', 'r') #polygon file: want to calculate percent cover of this lc type within study_area, and within areaKM2 (next variable) of each focal point pol = grassland.next() gl = MultiPolygon([shape(pol['geometry']) for pol in grassland]) areaKM2 = 20 #hyp home range size of specie of interest with traps as input: #calculate initial area in meters, set radius areaM2 = areaKM2 * 1000000 r = (math.sqrt(areaM2/math.pi)) #begin buffering and calculating available area (i.e. within study area) for each point for point in input: pt_buffer = shape(point['geometry']).buffer(r) avail_area = pt_buffer.intersection(sa).area #check and adjust radius of buffer until it covers desired available area within study area while avail_area < areaM2: r += 300 pt_buffer = shape(point['geometry']).buffer(r) avail_area = pt_buffer.intersection(sa).area #then, calulate percent cover of land cover type of interest within adjusted buffer area #print to check avail_lc_type = pt_buffer.intersection(gl).area perc_cov = (avail_lc_type/areaM2) * 100 print perc_cov @MWrenn - Thanks for the suggestion. I was able to profile my code and came up with this: 55.3555590078 415 function calls (365 primitive calls) in 48.633 seconds Ordered by: standard name ncalls tottime percall cumtime percall filename:lineno(function) 1 0.001 0.001 48.632 48.632 <module1>:15(neighb_func) 1 0.000 0.000 48.633 48.633 <string>:1(<module>) 2 0.000 0.000 0.001 0.000 <string>:531(write) 1 0.000 0.000 0.000 0.000 __init__.py:1118(debug) 1 0.000 0.000 0.000 0.000 __init__.py:1318(getEffectiveLevel) 1 0.000 0.000 0.000 0.000 __init__.py:1332(isEnabledFor) 1 0.000 0.000 0.000 0.000 _abcoll.py:483(update) 3 0.000 0.000 0.000 0.000 _weakrefset.py:68(__contains__) 1 0.000 0.000 0.000 0.000 abc.py:128(__instancecheck__) 1 0.000 0.000 0.000 0.000 abc.py:148(__subclasscheck__) 11 0.000 0.000 0.000 0.000 base.py:195(_is_empty) 11 0.003 0.000 0.003 0.000 base.py:202(empty) 7 0.000 0.000 0.003 0.000 base.py:212(__del__) 25 0.000 0.000 0.000 0.000 base.py:231(_geom) 2 0.000 0.000 0.000 0.000 base.py:235(_geom) 3 0.000 0.000 0.000 0.000 base.py:313(geometryType) 3 0.000 0.000 0.000 0.000 base.py:316(type) 3 0.000 0.000 0.001 0.000 base.py:383(area) 2 0.000 0.000 0.001 0.000 base.py:443(buffer) 8 0.000 0.000 0.000 0.000 base.py:46(geometry_type_name) 5 0.000 0.000 0.000 0.000 base.py:52(geom_factory) 3 0.000 0.000 48.626 16.209 base.py:529(intersection) 10 0.000 0.000 0.000 0.000 brine.py:106(_dump_int) 2 0.000 0.000 0.000 0.000 brine.py:150(_dump_str) 12/2 0.000 0.000 0.000 0.000 brine.py:179(_dump_tuple) 24/2 0.000 0.000 0.000 0.000 brine.py:202(_dump) 2 0.000 0.000 0.000 0.000 brine.py:332(dump) 10/2 0.000 0.000 0.000 0.000 brine.py:360(dumpable) 14/8 0.000 0.000 0.000 0.000 brine.py:369(<genexpr>) 2 0.000 0.000 0.000 0.000 channel.py:56(send) 1 0.000 0.000 0.000 0.000 collection.py:186(filter) 1 0.000 0.000 0.000 0.000 collection.py:274(__iter__) 1 0.000 0.000 0.000 0.000 collection.py:364(closed) 1 0.000 0.000 0.000 0.000 collections.py:37(__init__) 25 0.000 0.000 0.000 0.000 collections.py:53(__setitem__) 4 0.000 0.000 0.000 0.000 compat.py:17(BYTES_LITERAL) 2 0.000 0.000 0.000 0.000 coords.py:20(required) 2 0.000 0.000 0.000 0.000 geo.py:20(shape) 5 0.000 0.000 0.000 0.000 geos.py:484(errcheck_predicate) 8 0.000 0.000 0.000 0.000 impl.py:43(__getitem__) 2 0.000 0.000 0.000 0.000 point.py:124(_set_coords) 2 0.000 0.000 0.000 0.000 point.py:188(geos_point_from_py) 2 0.000 0.000 0.000 0.000 point.py:37(__init__) 2 0.000 0.000 0.001 0.000 protocol.py:220(_send) 2 0.000 0.000 0.001 0.000 protocol.py:227(_send_request) 2 0.000 0.000 0.000 0.000 protocol.py:241(_box) 2 0.000 0.000 0.001 0.000 protocol.py:438(_async_request) 2 0.000 0.000 0.000 0.000 stream.py:173(write) 11 0.000 0.000 0.000 0.000 topology.py:14(_validate) 3 0.000 0.000 0.000 0.000 topology.py:33(__call__) 3 48.626 16.209 48.626 16.209 topology.py:40(__call__) 2 0.001 0.000 0.001 0.000 topology.py:57(__call__) 5 0.000 0.000 0.000 0.000 utf_8.py:15(decode) 5 0.000 0.000 0.000 0.000 {__import__} 5 0.000 0.000 0.000 0.000 {_codecs.utf_8_decode} 3 0.000 0.000 0.000 0.000 {_ctypes.byref} 6/2 0.000 0.000 0.000 0.000 {all} 6 0.000 0.000 0.000 0.000 {getattr} 5 0.000 0.000 0.000 0.000 {globals} 10 0.000 0.000 0.000 0.000 {hasattr} 5 0.000 0.000 0.000 0.000 {isinstance} 31 0.000 0.000 0.000 0.000 {len} 5 0.000 0.000 0.000 0.000 {locals} 1 0.000 0.000 0.000 0.000 {math.sqrt} 2 0.000 0.000 0.000 0.000 {method 'acquire' of 'thread.lock' objects} 24 0.000 0.000 0.000 0.000 {method 'append' of 'list' objects} 1 0.000 0.000 0.000 0.000 {method 'disable' of '_lsprof.Profiler' objects} 28 0.000 0.000 0.000 0.000 {method 'get' of 'dict' objects} 1 0.000 0.000 0.000 0.000 {method 'items' of 'dict' objects} 2 0.000 0.000 0.000 0.000 {method 'join' of 'str' objects} 2 0.000 0.000 0.000 0.000 {method 'lower' of 'str' objects} 5 0.000 0.000 0.000 0.000 {method 'pack' of 'Struct' objects} 2 0.000 0.000 0.000 0.000 {method 'release' of 'thread.lock' objects} 2 0.000 0.000 0.000 0.000 {method 'send' of '_socket.socket' objects} 2 0.000 0.000 0.000 0.000 {next} I'm still trying to figure out exactly what it means, but I think it's telling me that the `intersection` function is taking ~16 seconds each time it's called. So, per @gene 's suggestion, I'm working on using a spatial index. I just have to figure out how to get libspatialindex going so that I can use Rtrees.

Answer 1

from shapely.strtree import STRtree

# Assuming you have your polygons as a list

grassland_polygons = [shape(pol['geometry']) for pol in grassland]

# Build the spatial index

spatial_index = STRtree(grassland_polygons)

for point in input:

pt_geom = shape(point['geometry'])

pt_buffer = pt_geom.buffer(r)

# Query spatial index for polygons that intersect the buffer

candidate_polygons = spatial_index.query(pt_buffer)

# Now intersect only these polygons and sum their areas

intersected_area = sum(poly.intersection(pt_buffer).area for poly in candidate_polygons)

perc_cov = (intersected_area / areaM2) * 100

print(perc_cov)