Shapely deprecation warning with combination multipolygon and polygon how to solve

I am trying to plot some data on a map, where counties are colored in based on certain values. This worked fine before (first picture), but now leaves the plots empty (second picture) with many Shapely deprecation warnings.

Reading up on the issue with the deprecation warning both here and here, I have tried using .geoms, but the data seems to be a combination of MultiPolygons and Polygons. This prevents me from using .geoms effectively.

A slice of the problematic data:

[<shapely.geometry.multipolygon.MultiPolygon object at 0x70ac02a82980>]
Length: 1, dtype: geometry
<GeometryArray>
[<shapely.geometry.multipolygon.MultiPolygon object at 0x70ac02a829b0>]
Length: 1, dtype: geometry
<GeometryArray>
[<shapely.geometry.multipolygon.MultiPolygon object at 0x70ac02a82a10>]
Length: 1, dtype: geometry
<GeometryArray>
[<shapely.geometry.polygon.Polygon object at 0x70ac02a82a70>]
Length: 1, dtype: geometry
<GeometryArray>
[<shapely.geometry.polygon.Polygon object at 0x70ac02a82ad0>]
Length: 1, dtype: geometry
<GeometryArray>
[<shapely.geometry.polygon.Polygon object at 0x70ac02a82b30>]
Length: 1, dtype: geometry

How can I edit my code such that I get the plots that I need with Shapely 1.8.0 and up?

My code originally:

poly = [df.loc[df['NAME_0'] == 'Sweden']['geometry'].values[0]]
pad1 = 6  #padding, degrees unit
exts = [poly[0].bounds[0] - pad1, poly[0].bounds[2] + pad1, poly[0].bounds[1] - pad1, poly[0].bounds[3] + pad1];
msk = Polygon(rect_from_bound(*exts)).difference( poly[0].simplify(0.01) )
msk_stm  = st_proj.project_geometry (msk, projection)  # project geometry to the projection used by stamen

cs1 = axs[0].contourf(lons, lats, FWI_sum, 60, transform=projection, levels=levels, cmap=cmap, extend='both')
cs11 = axs[0].contourf(lons2, lats2, forest_output, 60, transform=projection, levels=levels2, cmap=cmap2, extend='both')
axs[0].coastlines()
axs[0].add_feature(feature.OCEAN, facecolor = facecolor_oceanmask, edgecolor=edgecolor_oceanmask, zorder=zorder_oceanmask)
axs[0].add_feature(feature.BORDERS)
axs[0].set_extent(extent, projection)
axs[0].axis('off')
axs[0].add_geometries(poly, crs=projection, facecolor=facecolor_polygon, edgecolor=edgecolor_polygon)
axs[0].add_geometries(msk_stm, st_proj, facecolor=facecolor_landmask, edgecolor=edgecolor_landmask, zorder=zorder_landmask)
for county in counties:
    poly2 = df2.loc[df2['name'] == county]['geometry'].values
    axs[0].add_geometries(poly2, crs=projection, facecolor='None', edgecolor='black', zorder=1)

cs12 = axs[1].contourf(lons2, lats2, forest_output, 60, transform=projection, levels=levels2, cmap=cmap2, extend='both')
axs[1].set_extent(extent, projection)
axs[1].add_feature(feature.BORDERS)
axs[1].add_feature(feature.OCEAN, facecolor = facecolor_oceanmask, edgecolor=edgecolor_oceanmask, zorder=zorder_oceanmask)
axs[1].coastlines(resolution='10m')
axs[1].axis('off')
axs[1].add_geometries(poly, crs=projection, facecolor=facecolor_polygon, edgecolor=edgecolor_polygon)
axs[1].add_geometries(msk_stm, st_proj, facecolor=facecolor_landmask, edgecolor=edgecolor_landmask, zorder=zorder_landmask)

for county, fireextentcounty_norm in zip(counties, fireextentcounties_norm):
    poly3 = df2.loc[df2['name'] == county]['geometry'].values
    rgba2 = cmap3(fireextentcounty_norm)
    axs[1].add_geometries(poly3, crs=projection, facecolor=rgba2, edgecolor='black', zorder=1)

cs13 = axs[2].contourf(lons2, lats2, forest_output, 60, transform=projection, levels=levels2, cmap=cmap2, extend='both')
axs[2].set_extent(extent, projection)
axs[2].add_feature(feature.BORDERS)
axs[2].add_feature(feature.OCEAN, facecolor = facecolor_oceanmask, edgecolor=edgecolor_oceanmask, zorder=zorder_oceanmask)
axs[2].coastlines(resolution='10m')
axs[2].axis('off')
axs[2].add_geometries(poly, crs=projection, facecolor=facecolor_polygon, edgecolor=edgecolor_polygon)
axs[2].add_geometries(msk_stm, st_proj, facecolor=facecolor_landmask, edgecolor=edgecolor_landmask, zorder=zorder_landmask)

for county, firenumbercounty_norm in zip(counties, firenumbercounties_norm):
    poly4 = df2.loc[df2['name'] == county]['geometry'].values
    rgba3 = cmap4(firenumbercounty_norm)
    axs[2].add_geometries(poly4, crs=projection, facecolor=rgba3, edgecolor='black', zorder=1)

Which gave these warnings:

Warning (from warnings module):
  File "/usr/lib/python3/dist-packages/cartopy/feature/__init__.py", line 217
    self._geoms = tuple(geometries)
ShapelyDeprecationWarning: Iteration over multi-part geometries is deprecated and will be removed in Shapely 2.0. Use the `geoms` property to access the constituent parts of a multi-part geometry.

Warning (from warnings module):
  File "/usr/lib/python3/dist-packages/cartopy/feature/__init__.py", line 217
    self._geoms = tuple(geometries)
ShapelyDeprecationWarning: __len__ for multi-part geometries is deprecated and will be removed in Shapely 2.0. Check the length of the `geoms` property instead to get the  number of parts of a multi-part geometry.

Warning (from warnings module):
  File "/usr/lib/python3/dist-packages/cartopy/feature/__init__.py", line 217
    self._geoms = tuple(geometries)
ShapelyDeprecationWarning: Iteration over multi-part geometries is deprecated and will be removed in Shapely 2.0. Use the `geoms` property to access the constituent parts of a multi-part geometry.

Warning (from warnings module):
  File "/usr/lib/python3/dist-packages/cartopy/feature/__init__.py", line 217
    self._geoms = tuple(geometries)
ShapelyDeprecationWarning: __len__ for multi-part geometries is deprecated and will be removed in Shapely 2.0. Check the length of the `geoms` property instead to get the  number of parts of a multi-part geometry.

Warning (from warnings module):
  File "/usr/lib/python3/dist-packages/cartopy/feature/__init__.py", line 217
    self._geoms = tuple(geometries)
ShapelyDeprecationWarning: Iteration over multi-part geometries is deprecated and will be removed in Shapely 2.0. Use the `geoms` property to access the constituent parts of a multi-part geometry.

Warning (from warnings module):
  File "/usr/lib/python3/dist-packages/cartopy/feature/__init__.py", line 217
    self._geoms = tuple(geometries)
ShapelyDeprecationWarning: __len__ for multi-part geometries is deprecated and will be removed in Shapely 2.0. Check the length of the `geoms` property instead to get the  number of parts of a multi-part geometry.

My code now:

poly = [df.loc[df['NAME_0'] == 'Sweden']

['geometry'].values[0]]
pad1 = 6  #padding, degrees unit
exts = [poly[0].bounds[0] - pad1, poly[0].bounds[2] + pad1, poly[0].bounds[1] - pad1, poly[0].bounds[3] + pad1];
msk = Polygon(rect_from_bound(*exts)).difference( poly[0].simplify(0.01) )
msk_stm  = st_proj.project_geometry (msk, projection)  # project geometry to the projection used by stamen

cs1 = axs[0].contourf(lons, lats, FWI_sum, 60, transform=projection, levels=levels, cmap=cmap, extend='both')
cs11 = axs[0].contourf(lons2, lats2, forest_output, 60, transform=projection, levels=levels2, cmap=cmap2, extend='both')
axs[0].coastlines()
axs[0].add_feature(feature.OCEAN, facecolor = facecolor_oceanmask, edgecolor=edgecolor_oceanmask, zorder=zorder_oceanmask)
axs[0].add_feature(feature.BORDERS)
axs[0].set_extent(extent, projection)
axs[0].axis('off')
axs[0].add_geometries(poly[0].geoms, crs=projection, facecolor=facecolor_polygon, edgecolor=edgecolor_polygon)
axs[0].add_geometries(msk_stm, st_proj, facecolor=facecolor_landmask, edgecolor=edgecolor_landmask, zorder=zorder_landmask)
for county in counties:
    poly2 = df2.loc[df2['name'] == county]['geometry'].values
    print(poly2[0].geoms)
    axs[0].add_geometries(poly2[0].geoms, crs=projection, facecolor='None', edgecolor='black', zorder=1)

cs12 = axs[1].contourf(lons2, lats2, forest_output, 60, transform=projection, levels=levels2, cmap=cmap2, extend='both')
axs[1].set_extent(extent, projection)
axs[1].add_feature(feature.BORDERS)
axs[1].add_feature(feature.OCEAN, facecolor = facecolor_oceanmask, edgecolor=edgecolor_oceanmask, zorder=zorder_oceanmask)
axs[1].coastlines(resolution='10m')
axs[1].axis('off')
axs[1].add_geometries(poly[0].geoms, crs=projection, facecolor=facecolor_polygon, edgecolor=edgecolor_polygon)
axs[1].add_geometries(msk_stm, st_proj, facecolor=facecolor_landmask, edgecolor=edgecolor_landmask, zorder=zorder_landmask)

for county, fireextentcounty_norm in zip(counties, fireextentcounties_norm):
    poly3 = df2.loc[df2['name'] == county]['geometry'].values
    rgba2 = cmap3(fireextentcounty_norm)
    axs[1].add_geometries(poly3[0].geoms, crs=projection, facecolor=rgba2, edgecolor='black', zorder=1)

cs13 = axs[2].contourf(lons2, lats2, forest_output, 60, transform=projection, levels=levels2, cmap=cmap2, extend='both')
axs[2].set_extent(extent, projection)
axs[2].add_feature(feature.BORDERS)
axs[2].add_feature(feature.OCEAN, facecolor = facecolor_oceanmask, edgecolor=edgecolor_oceanmask, zorder=zorder_oceanmask)
axs[2].coastlines(resolution='10m')
axs[2].axis('off')
axs[2].add_geometries(poly[0].geoms, crs=projection, facecolor=facecolor_polygon, edgecolor=edgecolor_polygon)
axs[2].add_geometries(msk_stm, st_proj, facecolor=facecolor_landmask, edgecolor=edgecolor_landmask, zorder=zorder_landmask)

for county, firenumbercounty_norm in zip(counties, firenumbercounties_norm):
    poly4 = df2.loc[df2['name'] == county]['geometry'].values
    rgba3 = cmap4(firenumbercounty_norm)
    axs[2].add_geometries(poly4[0].geoms, crs=projection, facecolor=rgba3, edgecolor='black', zorder=1)

Which gives these errors:

Warning (from warnings module):
  File "/usr/lib/python3/dist-packages/cartopy/feature/__init__.py", line 217
    self._geoms = tuple(geometries)
ShapelyDeprecationWarning: Iteration over multi-part geometries is deprecated and will be removed in Shapely 2.0. Use the `geoms` property to access the constituent parts of a multi-part geometry.

Warning (from warnings module):
  File "/usr/lib/python3/dist-packages/cartopy/feature/__init__.py", line 217
    self._geoms = tuple(geometries)
ShapelyDeprecationWarning: __len__ for multi-part geometries is deprecated and will be removed in Shapely 2.0. Check the length of the `geoms` property instead to get the  number of parts of a multi-part geometry.
<shapely.geometry.base.GeometrySequence object at 0x74aa047dbf40>
<shapely.geometry.base.GeometrySequence object at 0x74aa047dbf10>
<shapely.geometry.base.GeometrySequence object at 0x74aa047fcee0>
<shapely.geometry.base.GeometrySequence object at 0x74aa047fdd50>
<shapely.geometry.base.GeometrySequence object at 0x74aa047fe890>
Traceback (most recent call last):
  File "/usr/lib/python3.10/idlelib/run.py", line 578, in runcode
    exec(code, self.locals)
  File "/home/els/Nextcloud/Documents/PhD/damage_index/FWImodel/Plotting/PlotYearSum.py", line 182, in <module>
    print(poly2[0].geoms)
AttributeError: 'Polygon' object has no attribute 'geoms'. Did you mean: '_geom'?

There are plans to make a change in shapely so single Polygons also have a .geoms property, but this isn't ready yet (link).

Until this change lands, the easiest workaround is to use shapely.get_parts. E.g. use shapely.get_parts(poly[0]) to replace poly[0].geoms.