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

bottom = 0.1
width = 0.7
height = 0.8
'''
left = 0.25
bottom = 0.01
width = 0.65
height = 0.98
rect = [left, bottom, width, height]

ax = plt.axes(rect, projection=ccrs.PlateCarree(), )
ax.set_extent((150, 155, -30, -23))

ax.coastlines(resolution='10m', zorder=2)

LAND_10m = cartopy.feature.NaturalEarthFeature('physical', 'land', '10m',
                                               edgecolor='face',
                                               facecolor=cartopy.feature.COLORS['land'])
RIVERS_10m = cartopy.feature.NaturalEarthFeature('physical', 'rivers_lake_centerlines', '10m',
                                                 edgecolor=cartopy.feature.COLORS['water'],
                                                 facecolor='none')
BORDERS2_10m = cartopy.feature.NaturalEarthFeature('cultural', 'admin_1_states_provinces',
                                                   '10m', edgecolor='black', facecolor='none')
ax.add_feature(LAND_10m)
ax.add_feature(RIVERS_10m)
ax.add_feature(BORDERS2_10m, edgecolor='grey')
ax.stock_img()
# stock image is good enough for example, but OCEAN_10m could be used, but very slow
#       ax.add_feature(OCEAN_10m)

ax.gridlines(draw_labels=True, xlocs=[150, 152, 154, 155])

# Convert the fraction value into a code of 0-8, which can be used to pull out
# the appropriate symbol
data['cloud_coverage'] = (8 * data_arr['cloud_fraction']).fillna(10).values.astype(int)

# Map weather strings to WMO codes, which we can use to convert to symbols
# Only use the first symbol if there are multiple
wx_text = data_arr['weather'].fillna('')
data['present_weather'] = [wx_code_map[s.split()[0] if ' ' in s else s] for s in wx_text]

###########################################
# All the data wrangling is finished, just need to set up plotting and go:
# Set up the map projection and set up a cartopy feature for state borders
proj = ccrs.LambertConformal(central_longitude=-95, central_latitude=35,
                             standard_parallels=[35])
state_boundaries = feat.NaturalEarthFeature(category='cultural',
                                            name='admin_1_states_provinces_lines',
                                            scale='110m', facecolor='none')

###########################################
# The payoff
# ----------

# Change the DPI of the resulting figure. Higher DPI drastically improves the
# look of the text rendering
plt.rcParams['savefig.dpi'] = 255

# Create the figure and an axes set to the projection
fig = plt.figure(figsize=(20, 10))
add_metpy_logo(fig, 1080, 290, size='large')
ax = fig.add_subplot(1, 1, 1, projection=proj)

def main():
    # Define the two coordinate systems with different ellipses.
    wgs84 = ccrs.PlateCarree(globe=ccrs.Globe(datum='WGS84',
                                              ellipse='WGS84'))
    sphere = ccrs.PlateCarree(globe=ccrs.Globe(datum='WGS84',
                                               ellipse='sphere'))

    # Define the coordinate system of the data we have from Natural Earth and
    # acquire the 1:10m physical coastline shapefile.
    geodetic = ccrs.Geodetic(globe=ccrs.Globe(datum='WGS84'))
    dataset = cfeature.NaturalEarthFeature(category='physical',
                                           name='coastline',
                                           scale='10m')

    # Create a Stamen map tiler instance, and use its CRS for the GeoAxes.
    tiler = Stamen('terrain-background')
    fig = plt.figure()
    ax = fig.add_subplot(1, 1, 1, projection=tiler.crs)
    ax.set_title('The effect of incorrectly referencing the Solomon Islands')

    # Pick the area of interest. In our case, roughly the Solomon Islands, and
    # get hold of the coastlines for that area.
    extent = [155, 163, -11.5, -6]
    ax.set_extent(extent, geodetic)
    geoms = list(dataset.intersecting_geometries(extent))

    # Add the Stamen aerial imagery at zoom level 7.

###########################################
# The payoff
# ----------

# Change the DPI of the resulting figure. Higher DPI drastically improves the
# look of the text rendering.
plt.rcParams['savefig.dpi'] = 255

# Create the figure and an axes set to the projection.
fig = plt.figure(figsize=(20, 10))
add_metpy_logo(fig, 1080, 290, size='large')
ax = fig.add_subplot(1, 1, 1, projection=proj)

# Set up a cartopy feature for state borders.
state_boundaries = feat.NaturalEarthFeature(category='cultural',
                                            name='admin_1_states_provinces_lines',
                                            scale='110m', facecolor='none')

# Add some various map elements to the plot to make it recognizable.
ax.add_feature(feat.LAND, zorder=-1)
ax.add_feature(feat.OCEAN, zorder=-1)
ax.add_feature(feat.LAKES, zorder=-1)
ax.coastlines(resolution='110m', zorder=2, color='black')
ax.add_feature(state_boundaries, edgecolor='black')
ax.add_feature(feat.BORDERS, linewidth=2, edgecolor='black')

# Set plot bounds
ax.set_extent((-118, -73, 23, 50))

#
# Here's the actual station plot

"but the cartopy package is not installed")
            
    if isinstance(subplot, gridspec.SubplotSpec):
        ax = plt.subplot(subplot, projection=crs)
    else:
        ax = plt.subplot(subplot[0], subplot[1], subplot[2], projection=crs)

    ax.add_feature(cfeature.NaturalEarthFeature("physical", "ocean", scale = "50m" if scale is "10m" else scale,
        edgecolor="none", facecolor=np.array([0.59375, 0.71484375, 0.8828125])), zorder=0)
    ax.add_feature(cfeature.NaturalEarthFeature("physical", "land",
       scale=scale, edgecolor="none", facecolor=np.array([0.9375, 0.9375, 0.859375])), zorder=0)
    ax.add_feature(cfeature.NaturalEarthFeature("physical", "coastline", scale=scale,
        edgecolor="black", facecolor="none", linewidth=lw), zorder=2)
    ax.add_feature(cfeature.NaturalEarthFeature("physical", "lakes", scale=scale,
        edgecolor="none", facecolor=np.array([0.59375, 0.71484375, 0.8828125])), zorder=0)
    ax.add_feature(cfeature.NaturalEarthFeature("physical", "rivers_lake_centerlines",
        scale=scale, edgecolor=np.array([ 0.59375, 0.71484375, 0.8828125]),
        facecolor="none"), zorder=0)
    ax.add_feature(cfeature.NaturalEarthFeature("cultural", "admin_0_boundary_lines_land",
        scale=scale, edgecolor="black", facecolor="none", linewidth=lw), zorder=2)

    if drawlonlatlines:
        ax.gridlines(crs=ccrs.PlateCarree())
    
    ax.set_extent(extent, crs)

    return ax

land color and ``zorder=0.9``.

    Returns
    -------
    ax : axes handle

    Note
    ----
    plot internally calls :py:func:`Regions.plot_regions`.

    """
    import cartopy.crs as ccrs
    import cartopy.feature as cfeature
    import matplotlib.pyplot as plt

    NEF = cfeature.NaturalEarthFeature

    if proj is None:
        proj = ccrs.PlateCarree()

    if ax is None:
        ax = plt.axes(projection=proj)

    if ocean_kws is None:
        ocean_kws = dict(color=cfeature.COLORS["water"], zorder=0.9)

    if land_kws is None:
        land_kws = dict(color=cfeature.COLORS["land"], zorder=0.9)

    if coastline_kws is None:
        coastline_kws = dict(color="0.4", lw=0.5)

def cartopy_coastlines(self):
        """ Get coastlines using cartopy. """
        return cartopy.feature.NaturalEarthFeature(
            category='physical', name='coastline', scale='10m',
            facecolor='none')

max_lat = np.max(self.Radar.fields[0].lat)
        else:
            min_lon, max_lon, min_lat, max_lat = extend
        XLON = np.arange(min_lon, max_lon, 0.01)
        YLAT = np.arange(min_lat, max_lat, 0.01)
        #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:

Function for creating himawari-8 image plots, called by create_plot()
        when cloud fraction is the selected plot type.
        '''

        simim_cube = self.datasets['simim']['data'][self.current_type][self.current_time]
        lats = simim_cube.coords('grid_latitude')[0].points
        lons = simim_cube.coords('grid_longitude')[0].points
        self.main_plot = self.current_axes.pcolormesh(lons,
                                                      lats,
                                                      simim_cube.data,
                                                      cmap=self.plot_options[self.current_type]['cmap'],
                                                      norm=self.plot_options[self.current_type]['norm']
                                                      )

        # Add coastlines to the map created by contourf
        coastline_50m = cartopy.feature.NaturalEarthFeature('physical',
                                                            'coastline',
                                                            '50m',
                                                            edgecolor='g',
                                                            alpha=0.5,
                                                            facecolor='none')

        self.current_axes.add_feature(coastline_50m)
        self.current_axes.set_extent((90, 154, -18, 30))
        self.update_title()

max_lat = np.max(self.Radar.fields[0].lat)
        else:
            min_lon, max_lon, min_lat, max_lat = extend
        XLON = np.arange(min_lon, max_lon, 0.01)
        YLAT = np.arange(min_lat, max_lat, 0.01)
        # ax.set_aspect("equal")
        self.Radar.add_product_CAPPI_lonlat(XLON, YLAT, level_height)
        radar_data = self.Radar.product["CAPPI_geo_%d" % level_height].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: