Mask and Plot Remote Sensing Data with EarthPy

Learn how to mask out pixels in a raster dataset. This example shows how to apply a cloud mask to Landsat 8 data.

Plotting with EarthPy


Below we walk through a typical workflow using Landsat data with EarthPy.

The example below uses Landsat 8 data. In the example below, the landsat_qa layer is the quality assurance data layer that comes with Landsat 8 to identify pixels that may represent cloud, shadow and water. The mask values used below are suggested values associated with the landsat_qa layer that represent pixels with clouds and cloud shadows.

Import Packages

To begin, import the needed packages. You will use a combination of several EarthPy modules including spatial, plot and mask.

from glob import glob
import os
import matplotlib.pyplot as plt
import rasterio as rio
from rasterio.plot import plotting_extent
import earthpy as et
import earthpy.spatial as es
import earthpy.plot as ep
import earthpy.mask as em

# Get data and set your home working directory
data ="cold-springs-fire")

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(, "earth-analytics"))

# Stack the landsat bands
# This creates a numpy array with each "layer" representing a single band
landsat_paths_pre = glob(
arr_st, meta = es.stack(landsat_paths_pre)

# Import the landsat qa layer
) as landsat_pre_cl:
    landsat_qa =
    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.


Customize Histogram Plot with Titles and Colors

        "Band 1",
        "Band 2",
        "Band 3",
        "Band 4",
        "Band 5",
        "Band 6",
        "Band 7",

View Single Band Plots

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


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.

    title="The Landsat QA Layer Comes with Landsat Data\n It can be used to remove clouds and shadows",

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]
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
    arr_ma, rgb=[4, 3, 2], title="Array with Clouds and Shadows Masked"

Total running time of the script: ( 0 minutes 1.597 seconds)

Gallery generated by Sphinx-Gallery