I have two geojson files, each with a list of polygons, that I can import easily with shapely:
import json
from shapely.geometry import shape, GeometryCollection, Point
with open("background.json") as f:
features = json.load(f)["features"]
background = GeometryCollection([shape(feature["geometry"]).buffer(0) for feature in features])
with open("foreground.json") as f:
features = json.load(f)["features"]
foreground = GeometryCollection([shape(feature["geometry"]).buffer(0) for feature in features])
The two sets of polygons represent two "layers". All polygons of the foreground
are inside a bigger polygon of the background
. Each polygon of the background
may contain zero, one, two, or any number of polygon of the foreground
.
I want to know in which polygon of the background every polygon of the foreground is contained. A trivial algorithm to do that would consist to loop on all polygons of the foreground f
, and test if f
is within every polygon b
of the background layer, and stops when one is found.
However this would be terrible in terms of performance, with a O(Nf*Nb)
scaling. Would there be an easy way to do this faster with the algorithms provided by shapely
?