Secure your code as it's written. Use Snyk Code to scan source code in minutes - no build needed - and fix issues immediately.
###############################################################################
# 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,