Connecting MultiPolygon Edges to Lines

This post is another Spatial Thoughts Academy Weekly Challenge solutions. The challenge is to find the edge of the polygon in a set of building footprints whose center point is closest to a street and then connect it with the nearest point on the street.  Shapely and Fiona were of course up to the challenge.  The basic algorithm is to grab the exterior-most ring of the building footprints, break out the edges, determine the midpoint of each edge, measure the distance to each road until the shortest distance is found, then grab the nearest point on the road, and finally draw a line between the edge midpoint and the nearest road point.  Shapely's nearest_points() method is really versatile here because you can throw any two geometries at it even if they are different types.  Fiona is used to read and write geopackages which I'm really liking lately as more an more people are using them finally. Here's the code:

	import fiona
	from shapely.geometry import mapping, shape, LineString, MultiLineString, MultiPolygon
	from shapely.ops import nearest_points
	# Did a quick manual check to see what the layer names are.
	# You could do this programmatically by checking the geometry
	# type.
	#>>>fiona.listlayers("city.gpkg")
	# ['buildings', 'roads', 'layer_styles']
	# Sample data: https://github.com/GeospatialPython/Learn/raw/master/city.gpkg
	# Now we open the buildings layer and then the roads layer
	with fiona.open("city.gpkg", layer="buildings") as buildings:
	with fiona.open("city.gpkg", layer="roads") as roads:
	# The list of lines containing our result
	face_lines = []
	# Grab the roads metadata file as a template for our results output
	meta = roads.meta
	# Now we loop through each building to begin our algorithm
	for building in buildings:
	# Covert the building geometry to a Shapely feature
	building_geometry = MultiPolygon(shape(building["geometry"])).geoms
	# Create the variables that monitor our loop checks until
	# we have the right answer.
	# The road closest to the building edge facing it.
	closest_road = None
	# The distance to the closest road to the building.
	closest_road_distance = None
	# The midpoint of the closest building edge to the road.
	closest_edge_midpoint = None
	# Loop through each road which is made up of segments
	for road in roads:
	# Convert the road to a Shapely Feature
	road_geometries = MultiLineString(shape(road["geometry"])).geoms
	# The buildings are multipolygons so we have to extract their
	# exterior rings.
	building_exterior_coords = []
	for poly in building_geometry:
	building_exterior_coords.extend(poly.exterior.coords)
	building_edges = LineString(building_exterior_coords)
	# Break the exterior ring feature back into points for
	# edge detection
	edge_points = list(building_edges.coords)
	# Loop through the edges
	for i,j in zip(edge_points, edge_points[1:]):
	# Compare edge midpoints to each road to see which is
	# the closest. This takes awhile.
	current_edge = LineString((i,j))
	current_edge_midpoint = current_edge.interpolate(0.5, normalized = True)
	for road_geometry in road_geometries:
	current_road_distance = current_edge_midpoint.distance(road_geometry)
	if closest_road:
	if current_road_distance < closest_road_distance:
	closest_road_distance = current_road_distance
	closest_road = road_geometry
	closest_edge_midpoint = current_edge_midpoint
	else:
	closest_road = road_geometry
	closest_road_distance = current_road_distance
	closest_edge_midpoint = current_edge_midpoint
	# Now that we found the closest edge midpoint to each road,
	# draw a line between the midpoint and the closest line
	# on the road.
	midpoint_to_road = nearest_points(closest_edge_midpoint, closest_road)
	face_lines.append(LineString(midpoint_to_road))
	# Now write out our result to a new geopackage!
	meta["schema"]["geometry"] = "LineString"
	meta["schema"]["properties"] = {"id": "int"}
	with fiona.open("faces.gpkg", "w", **meta) as faces:
	i = 0
	for face_line in face_lines:
	faces.write({
	"geometry": mapping(face_line),
	"properties": {"id": i}
	})
	i += 1

view raw face.py hosted with ❤ by GitHub

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