How to use the earthpy.plot.plot_rgb function in earthpy

###############################################################################
# Plot RGB Composite Image
# --------------------------
# You can use the ``plot_rgb()`` function from the ``earthpy.plot`` module to quickly
# plot three band composite images. For RGB composite images, you will plot the red,
# green, and blue bands, which are bands 4, 3, and 2, respectively, in the image
# stack you created. Python uses a zero-based index system, so you need to subtract
# a value of 1 from each index. Thus, the index for the red band is 3, green is 2,
# and blue is 1. These index values are provided to the ``rgb`` argument to identify
# the bands for the composite image.

# Create figure with one plot
fig, ax = plt.subplots(figsize=(12, 12))

# Plot red, green, and blue bands, respectively
ep.plot_rgb(arr_st, rgb=(3, 2, 1), ax=ax, title="Landsat 8 RGB Image")
plt.show()

###############################################################################
# Stretch Composite Images
# -------------------------
# Composite images can sometimes be dark if the pixel brightness values are
# skewed toward the value of zero. You can stretch the pixel brightness values
# in an image using the argument ``stretch=True`` to extend the values to the
# full 0-255 range of potential values to increase the visual contrast of the
# image. In addition, the ``str_clip`` argument allows you to specify how much of
# the tails of the data that you want to clip off. The larger the number, the
# more the data will be stretched or brightened.

# Create figure with one plot
fig, ax = plt.subplots(figsize=(12, 12))

# array from ``es.stack()`` and the Rasterio profile or metadata object. The function
# needs a single
# layer of a numpy array, which is why we use ``arr[0]``. The function also
# needs the spatial transformation for the Rasterio object, which can be acquired by accessing
# the ``"transform"`` key within the Rasterio Profile.

extent = plotting_extent(array[0], raster_prof["transform"])

################################################################################
# Plot Un-cropped Data
# ------------------------------
# You can see the boundary and the raster before the crop using ``ep.plot_rgb()``
# Notice that the data appear washed out.

fig, ax = plt.subplots(figsize=(12, 12))
ep.plot_rgb(
    array,
    ax=ax,
    stretch=True,
    extent=extent,
    str_clip=0.5,
    title="RGB Image of Un-cropped Raster",
)
plt.show()


################################################################################
# Explore the Range of Values in the Data
# ---------------------------------------
# You can explore the range of values found in the data using the EarthPy ``hist()``
# function. Do you notice any extreme values that may be impacting the stretch
# of the image?

# Plot of RGB composite image with polygon boundary
bound_utm13N.boundary.plot(ax=ax1, color="black", zorder=10)
ax1 = ep.plot_rgb(
    arr_st,
    rgb=(3, 2, 1),
    ax=ax1,
    stretch=True,
    extent=extent,
    str_clip=0.5,
    title="Landsat 8 RGB Image with Polygon Boundary",
)

# Plot of CIR composite image with polygon boundary
bound_utm13N.boundary.plot(ax=ax2, color="black", zorder=10)
ax2 = ep.plot_rgb(
    arr_st,
    rgb=(4, 3, 2),
    ax=ax2,
    stretch=True,
    extent=extent,
    str_clip=0.5,
    title="Landsat 8 CIR Image with Polygon Boundary",
)
plt.show()

# Stack All Bands
# ---------------
# Once the data are cropped, you are ready to create a new stack.

os.chdir(os.path.join(et.io.HOME, "earth-analytics"))

cropped_array, array_raster_profile = es.stack(band_paths_list, nodata=-9999)
crop_extent = plotting_extent(
    cropped_array[0], array_raster_profile["transform"]
)

# Plotting the cropped image
# sphinx_gallery_thumbnail_number = 5
fig, ax = plt.subplots(figsize=(12, 6))
crop_bound_utm13N.boundary.plot(ax=ax, color="red", zorder=10)
ep.plot_rgb(
    cropped_array,
    ax=ax,
    stretch=True,
    extent=crop_extent,
    title="Cropped Raster and Fire Boundary",
)
plt.show()

)
plt.show()

# Recreate extent object for the No Data array

extent_nodata = plotting_extent(
    array_nodata[0], raster_prof_nodata["transform"]
)

################################################################################
# Plot Un-cropped Data
# ------------------------------
# Plot the data again after the nodata values are removed.

fig, ax = plt.subplots(figsize=(12, 12))
ep.plot_rgb(
    array_nodata,
    ax=ax,
    stretch=True,
    extent=extent,
    str_clip=0.5,
    title="RGB image of Un-cropped Raster, No Data Value Selected",
)
plt.show()

#############################################################################
# Crop the Data
# ------------------
# Sometimes you have data for an area that is larger than your study area.
# It is more efficient to first crop the data to your study area before processing
# it in Python. The fastest and most efficient option is to crop each file
# individually, write out the cropped raster to a new file, and then stack

# ~~~~~~~~~~~~~~~~~~~~~
# Now apply the mask and plot the masked data. The mask applies to every band in your data.
# The mask values below are values documented in the Landsat 8 documentation that represent
# clouds and cloud shadows.

# Generate array of all possible cloud / shadow values
cloud_shadow = [328, 392, 840, 904, 1350]
cloud = [352, 368, 416, 432, 480, 864, 880, 928, 944, 992]
high_confidence_cloud = [480, 992]

# Mask the data
all_masked_values = cloud_shadow + cloud + high_confidence_cloud
arr_ma = em.mask_pixels(arr_st, landsat_qa, vals=all_masked_values)

# sphinx_gallery_thumbnail_number = 5
ep.plot_rgb(
    arr_ma, rgb=[4, 3, 2], title="Array with Clouds and Shadows Masked"
)
plt.show()

###############################################################################
# Stretch Composite Images
# -------------------------
# Composite images can sometimes be dark if the pixel brightness values are
# skewed toward the value of zero. You can stretch the pixel brightness values
# in an image using the argument ``stretch=True`` to extend the values to the
# full 0-255 range of potential values to increase the visual contrast of the
# image. In addition, the ``str_clip`` argument allows you to specify how much of
# the tails of the data that you want to clip off. The larger the number, the
# more the data will be stretched or brightened.

# Create figure with one plot
fig, ax = plt.subplots(figsize=(12, 12))

# Plot bands with stretched applied
ep.plot_rgb(
    arr_st,
    rgb=(3, 2, 1),
    ax=ax,
    stretch=True,
    str_clip=0.5,
    title="Landsat 8 RGB Image with Stretch Applied",
)
plt.show()

###############################################################################
# Plot Color Infrared (CIR) Composite Image
# ------------------------------------------
# For color infrared (CIR) composite images, you will plot the near-infrared (NIR),
# red, and green bands, which are bands 5, 4, 2, respectively. Once again, the
# Python index values will be the original band number minus 1, thus, 4, 3, and 2
# for NIR, red, and green, respectively.

)
plt.show()

###############################################################################
# Plot Color Infrared (CIR) Composite Image
# ------------------------------------------
# For color infrared (CIR) composite images, you will plot the near-infrared (NIR),
# red, and green bands, which are bands 5, 4, 2, respectively. Once again, the
# Python index values will be the original band number minus 1, thus, 4, 3, and 2
# for NIR, red, and green, respectively.

# Create figure with one plot
fig, ax = plt.subplots(figsize=(12, 12))

# Plot NIR, red, and green bands, respectively, with stretch
ep.plot_rgb(
    arr_st,
    rgb=(4, 3, 2),
    ax=ax,
    stretch=True,
    str_clip=0.5,
    title="Landsat 8 CIR Image with Stretch Applied",
)
plt.show()

#############################################################################
# Plot Boundary Over Composite Image
# -----------------------------------
# .. note::
#       If you are on windows, you may need to add the crs issue discussed above
#       here!
#

###############################################################################
# Create Figure with Multiple Axes or Subplots
# --------------------------------------------
# ```plot_rgb()`` has an ``ax=`` parameter which supports subplots. You can 
# create figures that contain multiple plots by creating multiple ax
# objects, one for each plot. You can also specify the number of rows and
# columns in which to display the plots. In the example below, the two plots
# will be displayed side-by-side along one row with two columns.

# Create figure with two plots
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 6))

# Plot of RGB composite image with polygon boundary
bound_utm13N.boundary.plot(ax=ax1, color="black", zorder=10)
ax1 = ep.plot_rgb(
    arr_st,
    rgb=(3, 2, 1),
    ax=ax1,
    stretch=True,
    extent=extent,
    str_clip=0.5,
    title="Landsat 8 RGB Image with Polygon Boundary",
)

# Plot of CIR composite image with polygon boundary
bound_utm13N.boundary.plot(ax=ax2, color="black", zorder=10)
ax2 = ep.plot_rgb(
    arr_st,
    rgb=(4, 3, 2),
    ax=ax2,
    stretch=True,