How to use earthpy - 10 common examples

def get_path(dataset):
    """ Construct a file path to a dataset.

    Parameters
    ----------
    dataset: string
        Name of a dataset to access (e.g., "epsg.json", or "RGB.byte.tif")

    Returns
    -------
    A file path (string) to the dataset
    """
    earthpy_path = os.path.split(earthpy.__file__)[0]
    data_dir = os.path.join(earthpy_path, "data")
    data_files = os.listdir(data_dir)
    if dataset not in data_files:
        raise KeyError(dataset + " not found in earthpy example data.")
    return os.path.join(data_dir, dataset)

----------
    dataset: string
        Name of a dataset to access (e.g., "epsg.json", or "RGB.byte.tif")

    Returns
    -------
    A file path (string) to the dataset

    Example
    -------

        >>> import earthpy.io as eio
        >>> eio.path_to_example('rmnp-dem.tif')
        '...rmnp-dem.tif'
    """
    earthpy_path = os.path.split(earthpy.__file__)[0]
    data_dir = os.path.join(earthpy_path, "example-data")
    data_files = os.listdir(data_dir)
    if dataset not in data_files:
        raise KeyError(dataset + " not found in earthpy example data.")
    return os.path.join(data_dir, dataset)

# Landsat 8 red band is band 4 at [3]
# Landsat 8 near-infrared band is band 5 at [4]
ndvi = es.normalized_diff(arr_st[4], arr_st[3])


###############################################################################
# Plot NDVI With Colorbar Legend of Continuous Values
# ----------------------------------------------------

# You can plot NDVI with a colorbar legend of continuous values using the
# ``plot_bands`` function from the ``earthpy.plot`` module.

titles = ["Landsat 8 - Normalized Difference Vegetation Index (NDVI)"]

# Turn off bytescale scaling due to float values for NDVI
ep.plot_bands(
    ndvi, cmap="RdYlGn", cols=1, title=titles, scale=False, vmin=-1, vmax=1
)


###############################################################################
# Classify NDVI
# -------------

# Next, you can classify the NDVI to categorize the results into useful classes.
# Values under 0 will be classified together as no vegetation. Additional classes
# will be created for bare area # and low, moderate, and high vegetation areas.

# Create classes and apply to NDVI results
ndvi_class_bins = [-np.inf, 0, 0.1, 0.25, 0.4, np.inf]
ndvi_landsat_class = np.digitize(ndvi, ndvi_class_bins)

# ---------------
# To begin, open your DEM layer as a numpy array using Rasterio. Below you set all
# terrain values < 0 to ``nan``. Next, plot the data using ``ep.plot_bands()``.

# Set the home directory and get the data for the exercise
os.chdir(os.path.join(et.io.HOME, "earth-analytics"))
dtm = "data/vignette-elevation/pre_DTM.tif"

# Open the DEM with Rasterio
with rio.open(dtm) as src:
    elevation = src.read(1)
    # Set masked values to np.nan
    elevation[elevation < 0] = np.nan

# Plot the data
ep.plot_bands(
    elevation,
    scale=False,
    cmap="gist_earth",
    title="DTM Without Hillshade",
    figsize=(10, 6),
)
plt.show()

####################################################################################
# Create the Hillshade
# --------------------
# Once the DEM is read in, call ``es.hillshade()`` to create the hillshade.

# Create and plot the hillshade with earthpy
hillshade = es.hillshade(elevation)

# -----------------------------------------------
# Next, have a look at the data, it looks like there is a large cloud that you
# may want to mask out.

ep.plot_bands(arr_st)
plt.show()


###############################################################################
# Mask the Data
# -----------------------------------------------
# You can use the EarthPy ``mask()`` function to handle this cloud.
# To begin you need to have a layer that defines the pixels that
# you wish to mask. In this case, the ``landsat_qa`` layer will be used.

ep.plot_bands(
    landsat_qa,
    title="The Landsat QA Layer Comes with Landsat Data\n It can be used to remove clouds and shadows",
)
plt.show()


###############################################################################
# Plot The Masked Data
# ~~~~~~~~~~~~~~~~~~~~~
# 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]

# The list must contain the same number of strings as there are bands in the stack.

titles = ["Ultra Blue", "Blue", "Green", "Red", "NIR", "SWIR 1", "SWIR 2"]
# sphinx_gallery_thumbnail_number = 1
ep.plot_bands(array_stack, title=titles)
plt.show()

##################################################################################
# Plot One Band in a Stack
# ------------------------
#
# If you give ``ep.plot_bands()`` a one dimensional numpy array,
# it will only plot that single band. You can turn off the
# colorbar using the ``cbar`` parameter (``cbar=False``).

ep.plot_bands(array_stack[4], cbar=False)
plt.show()

##################################################################################
# Turn Off Scaling
# -----------------
#
# ``ep.plot_bands()`` scales the imagery to a 0-255 scale by default. This range
# of values makes it easier for matplotlib to plot the data. To turn off
# scaling, set the scale parameter to ``False``. Below you
# plot NDVI with scaling turned off in order for the proper range of values
# (-1 to 1) to be displayed. You can use the ``cmap=`` parameter to adjust
# the colormap for the plot

NDVI = es.normalized_diff(array_stack[4], array_stack[3])
ep.plot_bands(NDVI, scale=False, cmap="RdYlGn")
plt.show()

)
landsat_path.sort()
array_stack, meta_data = es.stack(landsat_path, nodata=-9999)

###############################################################################
# Plot All Bands in a Stack
# --------------------------
#
# When you give ``ep.plot_bands()`` a three dimensional numpy array,
# it will plot all layers in the numpy array. You can create unique titles for
# each image by providing a list of titles using the ``title=`` parameter.
# The list must contain the same number of strings as there are bands in the stack.

titles = ["Ultra Blue", "Blue", "Green", "Red", "NIR", "SWIR 1", "SWIR 2"]
# sphinx_gallery_thumbnail_number = 1
ep.plot_bands(array_stack, title=titles)
plt.show()

##################################################################################
# Plot One Band in a Stack
# ------------------------
#
# If you give ``ep.plot_bands()`` a one dimensional numpy array,
# it will only plot that single band. You can turn off the
# colorbar using the ``cbar`` parameter (``cbar=False``).

ep.plot_bands(array_stack[4], cbar=False)
plt.show()

##################################################################################
# Turn Off Scaling
# -----------------

ep.plot_bands(array_stack[4], cbar=False)
plt.show()

##################################################################################
# Turn Off Scaling
# -----------------
#
# ``ep.plot_bands()`` scales the imagery to a 0-255 scale by default. This range
# of values makes it easier for matplotlib to plot the data. To turn off
# scaling, set the scale parameter to ``False``. Below you
# plot NDVI with scaling turned off in order for the proper range of values
# (-1 to 1) to be displayed. You can use the ``cmap=`` parameter to adjust
# the colormap for the plot

NDVI = es.normalized_diff(array_stack[4], array_stack[3])
ep.plot_bands(NDVI, scale=False, cmap="RdYlGn")
plt.show()

##################################################################################
# Adjust the Number of Columns for a Multi Band Plot
# ---------------------------------------------------
#
# The number of columns used while plotting multiple bands can be changed in order
# to change the arrangement of the images overall.

ep.plot_bands(array_stack, cols=2)
plt.show()

plt.show()

####################################################################################
# Change the Azimuth of the Sun
# -------------------------------
# The angle that sun light hits the landscape, impacts the shadows and highlights
# created on the landscape. You can adjust the azimuth values to adjust angle of the
# highlights and shadows that are created in your output hillshade. Azimuth numbers can
# range from 0 to 360 degrees, where 0 is due North. The default value for azimuth
# in ``es.hillshade()`` is 30 degrees.

# Change the azimuth of the hillshade layer
hillshade_azimuth_210 = es.hillshade(elevation, azimuth=210)

# Plot the hillshade layer with the modified azimuth
ep.plot_bands(
    hillshade_azimuth_210,
    scale=False,
    cbar=False,
    title="Hillshade with Azimuth set to 210 Degrees",
    figsize=(10, 6),
)
plt.show()

####################################################################################
# Change the Angle Altitude of the Sun
# -------------------------------------
# Another variable you can adjust for hillshade is what angle of the sun.
# The ``angle_altitude`` parameter values range from 0 to 90. 90 represents the sun
# shining from directly above the scene. The default value for ``angle_altitude`` in
# ``es.hillshade()`` is 30 degrees.

scale=False,
    cmap="gist_earth",
    title="DTM Without Hillshade",
    figsize=(10, 6),
)
plt.show()

####################################################################################
# Create the Hillshade
# --------------------
# Once the DEM is read in, call ``es.hillshade()`` to create the hillshade.

# Create and plot the hillshade with earthpy
hillshade = es.hillshade(elevation)

ep.plot_bands(
    hillshade,
    scale=False,
    cbar=False,
    title="Hillshade made from DTM",
    figsize=(10, 6),
)
plt.show()

####################################################################################
# Change the Azimuth of the Sun
# -------------------------------
# The angle that sun light hits the landscape, impacts the shadows and highlights
# created on the landscape. You can adjust the azimuth values to adjust angle of the
# highlights and shadows that are created in your output hillshade. Azimuth numbers can
# range from 0 to 360 degrees, where 0 is due North. The default value for azimuth
# in ``es.hillshade()`` is 30 degrees.