How to use the cartopy.feature function in Cartopy

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])

lon0, lon1, lat0, lat1 = ax.get_extent()

data_level['u'].metpy.unit_array[wind_slice, wind_slice].to('knots'),
         data_level['v'].metpy.unit_array[wind_slice, wind_slice].to('knots'),
         length=6)

# Plot heights and temperature as contours
h_contour = ax.contour(x, y, data_level['height'], colors='k', levels=range(5400, 6000, 60))
h_contour.clabel(fontsize=8, colors='k', inline=1, inline_spacing=8,
                 fmt='%i', rightside_up=True, use_clabeltext=True)
t_contour = ax.contour(x, y, data_level['temperature'], colors='xkcd:deep blue',
                       levels=range(-26, 4, 2), alpha=0.8, linestyles='--')
t_contour.clabel(fontsize=8, colors='xkcd:deep blue', inline=1, inline_spacing=8,
                 fmt='%i', rightside_up=True, use_clabeltext=True)

# Add geographic features
ax.add_feature(cfeature.LAND.with_scale('50m'), facecolor=cfeature.COLORS['land'])
ax.add_feature(cfeature.OCEAN.with_scale('50m'), facecolor=cfeature.COLORS['water'])
ax.add_feature(cfeature.STATES.with_scale('50m'), edgecolor='#c7c783', zorder=0)
ax.add_feature(cfeature.LAKES.with_scale('50m'), facecolor=cfeature.COLORS['water'],
               edgecolor='#c7c783', zorder=0)

# Set a title and show the plot
ax.set_title('500 hPa Heights (m), Temperature (\u00B0C), Humidity (%) at '
             + time[0].dt.strftime('%Y-%m-%d %H:%MZ').item())
plt.show()

hght_500 = data.variables['Geopotential_height_isobaric'][0, plev.index(500)]

# Calculate and smooth 1000-500 hPa thickness
thickness_1000_500 = gaussian_filter(hght_500 - hght_1000, sigma=3.0)

###############################
# Set map and data projections for use in mapping

# Set projection of map display
mapproj = ccrs.LambertConformal(central_latitude=45., central_longitude=-100.)

# Set projection of data
dataproj = ccrs.PlateCarree()

# Grab data for plotting state boundaries
states_provinces = cfeature.NaturalEarthFeature(
        category='cultural',
        name='admin_1_states_provinces_lakes',
        scale='50m',
        facecolor='none')

###############################
# Create figure and plot data
fig = plt.figure(1, figsize=(17., 11.))
ax = plt.subplot(111, projection=mapproj)

# Set extent and plot map lines
ax.set_extent([-145., -70, 20., 60.], ccrs.PlateCarree())
ax.coastlines('50m', edgecolor='black', linewidth=0.75)
ax.add_feature(states_provinces, edgecolor='black', linewidth=0.5)

# Plot thickness with multiple colors

def main():
    fig = plt.figure(figsize=[10, 5])
    ax1 = fig.add_subplot(1, 2, 1, projection=ccrs.SouthPolarStereo())
    ax2 = fig.add_subplot(1, 2, 2, projection=ccrs.SouthPolarStereo(),
                          sharex=ax1, sharey=ax1)
    fig.subplots_adjust(bottom=0.05, top=0.95,
                        left=0.04, right=0.95, wspace=0.02)

    # Limit the map to -60 degrees latitude and below.
    ax1.set_extent([-180, 180, -90, -60], ccrs.PlateCarree())

    ax1.add_feature(cfeature.LAND)
    ax1.add_feature(cfeature.OCEAN)

    ax1.gridlines()
    ax2.gridlines()

    ax2.add_feature(cfeature.LAND)
    ax2.add_feature(cfeature.OCEAN)

    # Compute a circle in axes coordinates, which we can use as a boundary
    # for the map. We can pan/zoom as much as we like - the boundary will be
    # permanently circular.
    theta = np.linspace(0, 2*np.pi, 100)
    center, radius = [0.5, 0.5], 0.5
    verts = np.vstack([np.sin(theta), np.cos(theta)]).T
    circle = mpath.Path(verts * radius + center)

"""
    Plot basic TXLA shelf background: coastline, bathymetry, meridians, etc
    Can optionally input grid (so it doesn't have to be loaded again)

    Args:
        proj
        extent: west lon, east lon, south lat, north lat for axes bounds
        pars: parallels to plot
        mers: meridians to plot
        hlevs: which depth contours to plot
        outline: west, east, north, south lines (left, right, top, bottom)
        res: 10m, 50m, 110m
    """

    mpl.rcParams.update({'font.size': fontsize})#,'font.weight': 'bold'})
    land = cfeature.NaturalEarthFeature('physical', 'land', res,
                                        edgecolor='face',
                                        facecolor=cfeature.COLORS['land'])
    states = cfeature.NaturalEarthFeature('cultural',
                                          name='admin_1_states_provinces_lines',
                                          scale=res, facecolor='none')

    if fig is None and figsize is None:
        fig = plt.figure(figsize=(9.4, 7.7), dpi=100)
    elif fig is None and figsize is not None:
        fig = plt.figure(figsize=figsize, dpi=100)

    if ax is None:
        ax = fig.add_subplot(111, projection=proj)

    ax.set_extent(extent, pc)
    gl = ax.gridlines(linewidth=0.2, color='gray', alpha=0.5, linestyle='-', draw_labels=True)

min_lon = np.min(lon)
        max_lon = np.max(lon)
        min_lat = np.min(lat)
        max_lat = np.max(lat)
    else:
        min_lon, max_lon, min_lat, max_lat = extend

    ax.set_aspect("equal")
    cmaps = plt.get_cmap(cmap)
    norm = BoundaryNorm(bounds, ncolors=cmaps.N, clip=True)
    pm = ax.pcolormesh(lon, lat, data, transform=transform, cmap=cmap, norm=norm, zorder=4, **kwargs)
    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(), 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=transform)
    ax.set_yticks(parallels, crs=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 = bounds

ax1.set_extent(extent1, crs=tifproj)

    #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='blue',facecolor='none')
    BORDERS2_10m = cartopy.feature.NaturalEarthFeature('cultural', 'admin_1_states_provinces',
                                                       '10m', edgecolor='red', facecolor='none',
                                                       linestyle='-')
    #COASTS_10m = cartopy.feature.COASTLINE('10m', edgecolor='grey', facecolor='none')

    #ax1.add_feature(LAND_10m, edgecolor='grey', zorder=1.2)
    #ax1.coastlines(resolution='10m', color='grey', zorder=2)
    ax1.add_feature(RIVERS_10m, zorder=1.2)
    ax1.add_feature(cartopy.feature.COASTLINE, edgecolor='gray', color='none', zorder=1.2)
    ax1.add_feature(BORDERS2_10m, zorder=1.2)
    ax1.stock_img()
    # stock image is good enough for example, but OCEAN_10m could be used, but very slow
    #       ax.add_feature(OCEAN_10m)

    print('mapextent given to get gridlines:')
    print(mapextent)

    # work out gridline positions
    xticks, yticks = get_gridlines(mapextent[0], mapextent[1], mapextent[2], mapextent[3], nticks=6)

    # plot the gridlines
    gl = ax1.gridlines(crs=tifproj, xlocs=xticks, ylocs=yticks, linestyle='--', color='grey',
                       alpha=1, linewidth=1, zorder=1.3)
    # add ticks
    ax1.set_xticks(xticks[1:-1], crs=tifproj)

cyclic_data, cyclic_lons = add_cyclic_point(coriolis, coord=lons)

#####################################

# sphinx_gallery_thumbnail_number = 2

# Works with matplotlib's built-in transform support.
fig = plt.figure(figsize=(10, 8))
ax = fig.add_subplot(1, 1, 1, projection=ccrs.Robinson())

# Sets the extent to cover the whole globe
ax.set_global()

# Add variety of features
ax.add_feature(cfeature.LAND)
ax.add_feature(cfeature.OCEAN)
ax.add_feature(cfeature.COASTLINE)

# Set negative contours to be solid instead of dashed
matplotlib.rcParams['contour.negative_linestyle'] = 'solid'
CS = ax.contour(cyclic_lons, lats, cyclic_data, 20, colors='tab:brown',
                transform=ccrs.PlateCarree())
ax.clabel(CS, inline=1, fontsize=10, fmt='%1.1f')

plt.show()

ds = xr.open_dataset(get_test_data('narr_example.nc', as_file_obj=False))
data_var = ds.metpy.parse_cf('Temperature')

x = data_var.x
y = data_var.y
im_data = data_var.isel(time=0).sel(isobaric=1000.)

fig = plt.figure(figsize=(14, 14))
ax = fig.add_subplot(1, 1, 1, projection=data_var.metpy.cartopy_crs)

ax.imshow(im_data, extent=(x.min(), x.max(), y.min(), y.max()),
          cmap='RdBu', origin='lower' if y[0] < y[-1] else 'upper')
ax.coastlines(color='tab:green', resolution='10m')
ax.add_feature(cfeature.LAKES.with_scale('10m'), facecolor='none', edgecolor='tab:blue')
ax.add_feature(cfeature.RIVERS.with_scale('10m'), edgecolor='tab:blue')

plt.show()