How to use the earthpy.spatial.stack function in earthpy

###############################################################################
# Import Example Data
# -------------------
#
# To get started, make sure your directory is set. Create a stack from all of the
# Landsat .tif files (one per band).

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

# Stack the Landsat 8 bands
# This creates a numpy array with each "layer" representing a single band
landsat_path = glob(
    "data/cold-springs-fire/landsat_collect/LC080340322016072301T1-SC20180214145802/crop/*band*.tif"
)
landsat_path.sort()
arr_st, meta = es.stack(landsat_path)


###############################################################################
# Calculate Normalized Difference Vegetation Index (NDVI)
# -------------------------------------------------------

# You can calculate NDVI for your dataset using the
# ``normalized_diff`` function from the ``earthpy.spatial`` module.
# Math will be calculated (b1-b2) / (b1 + b2).

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


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

raster_out_path = os.path.join(output_dir, "raster.tiff")

####################################################################################
# Stack the Bands
# ---------------------------
# The stack function has an optional output argument, where you can write the raster
# to a tiff file in a folder. If you want to use this functionality, make sure there
# is a folder to write your tiff file to.
# The Stack function also returns two object, an array and a RasterIO profile. Make
# sure to be catch both in variables.

# Stack Landsat bands

os.chdir(os.path.join(et.io.HOME, "earth-analytics"))
array, raster_prof = es.stack(stack_band_paths, out_path=raster_out_path)

####################################################################################
# Create Extent Object
# --------------------------------
# To get the raster extent, use the ``plotting_extent`` function on the
# 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 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?

ep.hist(array, title=["Band 1", "Band 2", "Band 3"])
plt.show()

###########################################################################
# No Data Option
# ---------------
# ``es.stack()`` can handle ``nodata`` values in a raster. To use this
# parameter, specify ``nodata=``. This will mask every pixel that contains
# the specified ``nodata`` value. The output will be a numpy masked array.

os.chdir(os.path.join(et.io.HOME, "earth-analytics"))
array_nodata, raster_prof_nodata = es.stack(stack_band_paths, nodata=-9999)

# View hist of data with nodata values removed
ep.hist(
    array_nodata,
    title=[
        "Band 1 - No Data Values Removed",
        "Band 2 - No Data Values Removed",
        "Band 3 - No Data Values Removed",
    ],
)
plt.show()

# Recreate extent object for the No Data array

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

#       If you are running this script on a Windows system, there is a
#       known bug with ``.to_crs()``, which is used in this script. If an error
#       occurs, you have to reset your os environment with the command
#       ``os.environ["PROJ_LIB"] = r"path-to-share-folder-in-environment"``.

# Get sample data from EarthPy and set working directory
data_path = et.data.get_data("vignette-landsat")
os.chdir(os.path.join(et.io.HOME, "earth-analytics"))

# Get list of bands and sort by ascending band number
landsat_bands_data_path = "data/vignette-landsat/LC08_L1TP_034032_20160621_20170221_01_T1_sr_band*[1-7]*.tif"
stack_band_paths = glob(landsat_bands_data_path)
stack_band_paths.sort()

# Create image stack and apply nodata value for Landsat
arr_st, meta = es.stack(stack_band_paths, nodata=-9999)

###############################################################################
# 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

###############################################################################
# Import Example Data
# ------------------------------
# To get started, make sure your directory is set. Create a stack from all of the
# Landsat .tif files (one per band) and import the ``landsat_qa`` layer which provides
# the locations of cloudy and shadowed pixels in the scene.

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

# Stack the landsat bands
# This creates a numpy array with each "layer" representing a single band
landsat_paths_pre = glob(
    "data/cold-springs-fire/landsat_collect/LC080340322016070701T1-SC20180214145604/crop/*band*.tif"
)
landsat_paths_pre.sort()
arr_st, meta = es.stack(landsat_paths_pre)

# Import the landsat qa layer
with rio.open(
    "data/cold-springs-fire/landsat_collect/LC080340322016070701T1-SC20180214145604/crop/LC08_L1TP_034032_20160707_20170221_01_T1_pixel_qa_crop.tif"
) as landsat_pre_cl:
    landsat_qa = landsat_pre_cl.read(1)
    landsat_ext = plotting_extent(landsat_pre_cl)

###############################################################################
# Plot Histogram of Each Band in Your Data
# ----------------------------------------
# You can view a histogram for each band in your dataset by using the
# ``hist()`` function from the ``earthpy.plot`` module.

ep.hist(arr_st)
plt.show()