Actually, I found difficulties to plot polygons with some holes. The code modified as below grants at least one hole, as I tested it.

As this should be the main post on stackoverflow for compatibility between Shapely 2.x and descartes, the contribution below includes EXTERIOR and INTERIOR coordinates fixes, in order to have a full integration.

Modify this:

if hasattr(polygon, 'geom_type'):  # Shapely
    ptype = polygon.geom_type
    if ptype == 'Polygon':
        polygon = [Polygon(polygon)]
    elif ptype == 'MultiPolygon':
        polygon = [Polygon(p) for p in polygon]

with this:

if hasattr(polygon, 'geom_type'):
    ptype = polygon.geom_type
    if ptype == 'Polygon':
        polygon = [polygon]
    elif ptype == 'MultiPolygon':
        polygon = list(polygon.geoms)

And this:

 if ptype == 'Polygon':
    polygon = [Polygon(polygon)]
 elif ptype == 'MultiPolygon':
    polygon = [Polygon(p) for p in polygon['coordinates']]
   

with this:

  if ptype == 'Polygon':
      polygon = [Polygon(polygon['coordinates'])]
  elif ptype == 'MultiPolygon':
      polygon = [Polygon(p) for p in polygon['coordinates']]

And finally this:

vertices = concatenate([
        concatenate([asarray(t.exterior)[:, :2]] +
                    [asarray(r)[:, :2] for r in t.interiors])
        for t in polygon])
codes = concatenate([
    concatenate([coding(t.exterior)] +
                [coding(r) for r in t.interiors]) for t in polygon])

With this:

vertices = concatenate([
        concatenate([asarray(t.exterior.coords)[:, :2]] +
                    [asarray(r.coords)[:, :2] for r in t.interiors if len(r.coords) > 0])
        for t in polygon
    ])
codes = concatenate([
    concatenate([coding(t.exterior)] +
                [coding(r) for r in t.interiors if len(r.coords) > 0])
    for t in polygon
    ])

It seems to work for me.