How to use the cartopy.feature.ShapelyFeature function in Cartopy

def draw(self):
            ax = self.ax
            if len(self.df.geometry) == 0:
                return ax

            edgecolor = kwargs.pop('edgecolor', 'black')
            facecolor = kwargs.pop('facecolor', 'None')
            zorder = kwargs.pop('zorder', -1)

            if self.projection:
                for geom in self.df.geometry:
                    features = ShapelyFeature([geom], ccrs.PlateCarree())
                    ax.add_feature(
                        features, facecolor=facecolor, edgecolor=edgecolor, zorder=zorder,
                        **kwargs
                    )
            else:
                for geom in df.geometry:
                    try:  # Duck test for MultiPolygon.
                        for subgeom in geom:
                            feature = descartes.PolygonPatch(
                                subgeom, facecolor=facecolor, edgecolor=edgecolor, zorder=zorder,
                                **kwargs
                            )
                            ax.add_patch(feature)
                    except (TypeError, AssertionError):  # Shapely Polygon.
                        feature = descartes.PolygonPatch(
                            geom, facecolor=facecolor, edgecolor=edgecolor, zorder=zorder,

#ax.set_aspect("equal")
        self.Radar.add_product_CR_lonlat(XLON, YLAT)
        radar_data = self.Radar.product["CR_geo"].values
        lon, lat = np.meshgrid(XLON, YLAT, indexing="ij")
        cmaps = plt.get_cmap(cmap)
        levels = MaxNLocator(nbins=cmap_bins).tick_values(vmin, vmax)
        norm = BoundaryNorm(levels, ncolors=cmaps.N, clip=True)
        pm = ax.pcolormesh(lon, lat, radar_data, transform=self.transform, cmap=cmap, norm=norm, zorder=4, **kwargs)
        #ax.set_extent([min_lon, max_lon, min_lat, max_lat], crs=self.transform)
        ax.add_feature(cfeature.OCEAN.with_scale('50m'), zorder=0)
        ax.add_feature(cfeature.NaturalEarthFeature('physical', 'land', '50m', \
                                                    edgecolor='none', facecolor="white"), zorder=1)
        ax.add_feature(cfeature.LAKES.with_scale('50m'), zorder=2)
        ax.add_feature(cfeature.RIVERS.with_scale('50m'), zorder=3)

        ax.add_feature(cfeature.ShapelyFeature(CN_shp_info.geometries(), self.transform, \
                                               edgecolor='k', facecolor='none'), linewidth=0.5, \
                       linestyle='-', zorder=5, alpha=0.8)
        parallels = np.arange(int(min_lat), np.ceil(max_lat) + 1, 1)
        meridians = np.arange(int(min_lon), np.ceil(max_lon) + 1, 1)
        ax.set_xticks(meridians, crs=self.transform)
        ax.set_yticks(parallels, crs=self.transform)
        lon_formatter = LongitudeFormatter()
        lat_formatter = LatitudeFormatter()
        ax.xaxis.set_major_formatter(lon_formatter)
        ax.yaxis.set_major_formatter(lat_formatter)
        if cbar:
            cb = plt.colorbar(mappable=pm, ax=ax, orientation=orientation)
            if cbar_ticks is None:
                ticks = levels
            else:
                ticks = cbar_ticks

colors = ['None']*len(df)

    # Manipulate trace_kwargs.
    if trace:
        if trace_kwargs is None:
            trace_kwargs = dict()
        if 'edgecolor' not in trace_kwargs.keys():
            trace_kwargs['edgecolor'] = 'lightgray'
        if 'facecolor' not in trace_kwargs.keys():
            trace_kwargs['facecolor'] = 'None'

    # Draw traces first, if appropriate.
    if trace:
        if projection:
            for polygon in df.geometry:
                features = ShapelyFeature([polygon], ccrs.PlateCarree())
                ax.add_feature(features, **trace_kwargs)
        else:
            for polygon in df.geometry:
                try:  # Duck test for MultiPolygon.
                    for subgeom in polygon:
                        feature = descartes.PolygonPatch(subgeom, **trace_kwargs)
                        ax.add_patch(feature)
                except (TypeError, AssertionError):  # Shapely Polygon.
                    feature = descartes.PolygonPatch(polygon, **trace_kwargs)
                    ax.add_patch(feature)

    # Finally, draw the scaled geometries.
    for value, color, polygon in zip(values, colors, df.geometry):
        scale_factor = dscale(value)
        scaled_polygon = shapely.affinity.scale(polygon, xfact=scale_factor, yfact=scale_factor)
        if projection:

# By often creates overlapping polygons, to keep smaller polygons from being hidden by possibly overlapping
        # larger ones we have to bring the smaller ones in front in the plotting order. This bit of code does that.
        sorted_indices = np.array(sorted(enumerate(gpd.GeoSeries(sectors).area.values),
                                         key=lambda tup: tup[1])[::-1])[:, 0].astype(int)
        sectors = np.array(sectors)[sorted_indices]
        values = np.array(values)[sorted_indices]

        # Generate a colormap.
        cmap = _continuous_colormap(values, cmap, vmin, vmax)
        colors = [cmap.to_rgba(value) for value in values]

        #  Draw.
        for sector, color in zip(sectors, colors):
            if projection:
                features = ShapelyFeature([sector], ccrs.PlateCarree())
                ax.add_feature(features, facecolor=color, **kwargs)
            else:
                try:  # Duck test for MultiPolygon.
                    for subgeom in sector:
                        feature = descartes.PolygonPatch(subgeom, facecolor=color, **kwargs)
                        ax.add_patch(feature)
                except (TypeError, AssertionError):  # Shapely Polygon.
                    feature = descartes.PolygonPatch(sector, facecolor=color, **kwargs)
                    ax.add_patch(feature)

        # Set extent.
        extrema = (bxmin, bxmax, bymin, bymax)
        _set_extent(ax, projection, extent, extrema)

    else:
        # Set reasonable defaults for the n-params if appropriate.

geoms = build_voronoi_polygons(self.df)
            # Must take the .values of the output GeoDataFrame because assign is index-aligned.
            # So self.df can have any index. But set_clip constructs a new GeoSeries with a fresh
            # descending (1..N) index. If self.df doesn't also have a 1..N index, the join will
            # be misaligned and/or nan values will be inserted. The easiest way to assign in an
            # index-naive (e.g. index-based) manner is to provide a numpy array instead of a
            # GeoSeries by taking .values.
            self.df = self.df.assign(
                geometry=self.set_clip(gpd.GeoDataFrame(geometry=geoms)).values
            )
            for color, geom in zip(self.colors, self.df.geometry):
                if geom.is_empty:  # do not plot data points that return empty due to clipping
                    continue

                if self.projection:
                    feature = ShapelyFeature([geom], ccrs.PlateCarree())
                    ax.add_feature(
                        feature, facecolor=color, edgecolor=edgecolor, **self.kwargs
                    )
                else:
                    feature = descartes.PolygonPatch(
                        geom, facecolor=color, edgecolor=edgecolor, **self.kwargs
                    )
                    ax.add_patch(feature)

            return ax

def draw(self):
            ax = self.ax

            if len(df.geometry) == 0:
                return ax

            if self.projection:
                for color, geom in zip(self.colors, df.geometry):
                    features = ShapelyFeature([geom], ccrs.PlateCarree())
                    ax.add_feature(features, facecolor=color, **self.kwargs)
            else:
                for color, geom in zip(self.colors, df.geometry):
                    try:  # Duck test for MultiPolygon.
                        for subgeom in geom:
                            feature = descartes.PolygonPatch(
                                subgeom, facecolor=color, **self.kwargs
                            )
                            ax.add_patch(feature)
                    except (TypeError, AssertionError):  # Shapely Polygon.
                        feature = descartes.PolygonPatch(
                            geom, facecolor=color, **self.kwargs
                        )
                        ax.add_patch(feature)

            return ax

def draw(self):
            ax = self.ax
            if len(self.df.geometry) == 0:
                return ax

            for value, color, polygon in zip(self.sizes, self.colors, self.df.geometry):
                scale_factor = value
                scaled_polygon = shapely.affinity.scale(
                    polygon, xfact=scale_factor, yfact=scale_factor
                )
                if self.projection is not None:
                    features = ShapelyFeature([scaled_polygon], ccrs.PlateCarree())
                    ax.add_feature(features, facecolor=color, **kwargs)
                else:
                    try:  # Duck test for MultiPolygon.
                        for subgeom in scaled_polygon:
                            feature = descartes.PolygonPatch(
                                subgeom, facecolor=color, **self.kwargs
                            )
                            ax.add_patch(feature)
                    except (TypeError, AssertionError):  # Shapely Polygon.
                        feature = descartes.PolygonPatch(
                            scaled_polygon, facecolor=color, **self.kwargs
                        )
                        ax.add_patch(feature)

            return ax

if len(df.geometry) == 0:
        return ax

    # Set extent.
    extrema = np.min(xs), np.max(xs), np.min(ys), np.max(ys)
    _set_extent(ax, projection, extent, extrema)

    if projection:
        if clip is None:
            sns.kdeplot(pd.Series([p.x for p in df.geometry]), pd.Series([p.y for p in df.geometry]),
                        transform=ccrs.PlateCarree(), ax=ax, **kwargs)
        else:
            sns.kdeplot(pd.Series([p.x for p in df.geometry]), pd.Series([p.y for p in df.geometry]),
                        transform=ccrs.PlateCarree(), ax=ax, **kwargs)
            clip_geom = _get_clip(ax.get_extent(crs=ccrs.PlateCarree()), clip)
            feature = ShapelyFeature([clip_geom], ccrs.PlateCarree())
            ax.add_feature(feature, facecolor=(1,1,1), linewidth=0, zorder=100)
    else:
        if clip is None:
            sns.kdeplot(pd.Series([p.x for p in df.geometry]), pd.Series([p.y for p in df.geometry]), ax=ax, **kwargs)
        else:
            clip_geom = _get_clip(ax.get_xlim() + ax.get_ylim(), clip)
            polyplot(gpd.GeoSeries(clip_geom),
                     facecolor='white', linewidth=0, zorder=100, extent=ax.get_xlim() + ax.get_ylim(), ax=ax)
            sns.kdeplot(pd.Series([p.x for p in df.geometry]), pd.Series([p.y for p in df.geometry]),
                        ax=ax, **kwargs)
    return ax

fig, ax = plt.subplots(figsize=(figsizex, figsizey),
                           subplot_kw=subplot_kw)

    # set a margin around the data
    ax.set_xmargin(0.05)
    ax.set_ymargin(0.10)

    # add a background image for rendering
    ax.stock_img()

    # show the data from the geotiff RGB image
    img = ax.imshow(rgbdata[:3, :, :].transpose((1, 2, 0)),
                    extent=extent, origin='upper')

    # read shapefile and plot it onto the tiff image map
    shape_feature = ShapelyFeature(Reader(shapefile).geometries(),
                                   crs=shapeproj, edgecolor='yellow',
                                   facecolor='none')
    ax.add_feature(shape_feature)

    # add a title
    plt.title(plottitle)

    # set map extent
    ax.set_extent(mapextent, tifproj)

    # add coastlines
    ax.coastlines(resolution='10m', color='navy', linewidth=1)

    # add lakes and rivers
    ax.add_feature(cartopy.feature.LAKES, alpha=0.5)
    ax.add_feature(cartopy.feature.RIVERS)

Returns
        -------
        A :class:`cartopy.mpl.feature_artist.FeatureArtist` instance
            The instance responsible for drawing the feature.

        Note
        ----
            Matplotlib keyword arguments can be used when drawing the feature.
            This allows standard Matplotlib control over aspects such as
            'facecolor', 'alpha', etc.


        """
        styler = kwargs.pop('styler', None)
        feature = cartopy.feature.ShapelyFeature(geoms, crs, **kwargs)
        return self.add_feature(feature, styler=styler)