Extracting MultiPolygons with Holes using Shapely and Fiona

SpatialThoughts.com recently posted a challenge on LinkedIn to extract only building footprints with holes from a city-wide dataset.  

What made the challenge interesting is the building footprints are MultiPolygons instead of simple Polygons which stumped a few people's attempts. But this challenge is very straight-forward Shapely to select the polygons and Fiona to read and write the data.

The correct feature count if you extract all of the building footprints with holes or polygon inner rings is 45.  The challenge author later demonstrated how to do it in QGIS using a series of expressions.  I think I had the only Python solution. One guy posted a successful answer in R.

The challenge data is a city extracted from OSM: https://lnkd.in/g7dVqzMx

Here is my solution:

	"""Extract Polygons with holes"""
	#https://github.com/GeospatialPython/Learn/raw/master/BuildingFootprintsOSM.gpkg
	import fiona
	from shapely.geometry import shape, MultiPolygon
	with fiona.open("BuildingFootprintsOSM.gpkg", "r") as city:
	# Collection of footprints with holes
	hole_footprints = []
	# Loop through footprints and check
	# for interior rings (holes). If at
	# least one is present, copy it
	# to our list
	for footprint in city:
	shapes = MultiPolygon(shape(footprint["geometry"]))
	for s in shapes.geoms:
	if len(s.interiors):
	hole_footprints.append(footprint)
	# Create a new geopackage and write out the hole features
	with fiona.open("holes.gpkg", "w", **city.meta) as holes:
	for hole_footprint in hole_footprints:
	holes.write(hole_footprint)

view raw holes.py hosted with ❤ by GitHub