Plotting multiple maps using GeoPandas

I usually plot multiple maps as following.

import pandas as pd
import geopandas as gpd
import matplotlib.pyplot as plt


fig, axs = plt.subplots(5, 5, figsize=(25, 25))
axs = axs.flatten()

for idx in range(len(regions)):
    regions.iloc[[idx]].plot(ax=axs[idx], color='#ffffff', edgecolor='#6a6a6a', linewidth=2)
    urbfab.iloc[[idx]].plot(ax=axs[idx], color='#000000', edgecolor='#000000')
    axs[idx].set_title(f"{regions.iloc[idx]['NAMEREG']} #{idx}")
    axs[idx].axis('off')
plot.show()

Spatial join the region names to the second df. Groupby region and plot

import geopandas as gpd
import matplotlib.pyplot as plt

#Read two shapefiles as geodataframes
regions = gpd.read_file(r'D:/gisdata/regions.shp')
buildings = gpd.read_file(r'D:/gisdata/buildings.shp')

region_name_column = "KOMMUNNAMN" #The column name identifying each region in the region df

#Spatial join region attribute to the buildings
buildings = buildings.sjoin(regions[[region_name_column, "geometry"]], how="inner", predicate="intersects")


#Create a plot with 3 rows and 3 columns
fig, axes = plt.subplots(nrows=3, ncols=3, figsize=(20,20))
axes = axes.flatten() #From 2D array with shape 3,3 to 1D array with shape 9. Easier to iterate over
[ax.set_axis_off() for ax in axes] #Hide all axes


#For each region, plot it and the buildings in it
for ax, (region_name, region_subframe) in zip(axes, regions.groupby(region_name_column)):
    region_subframe.plot(ax=ax, zorder=1, facecolor="None", edgecolor="red", linewidth=4)
    
    #Select buildings by the joined region attribute column1
    buildings_subframe = buildings.loc[buildings[region_name_column]==region_name]
    buildings_subframe.plot(ax=ax, zorder=2, facecolor="black", edgecolor="black")
    
    ax.set_title(region_name, fontsize=20)