{ "cells": [ { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "%matplotlib inline" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n# Mask and Plot Remote Sensing Data with EarthPy\n\nLearn how to mask out pixels in a raster dataset. This example shows how to apply a cloud mask to\nLandsat 8 data.\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Plotting with EarthPy\n\n

Note

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

\n\nThe example below uses Landsat 8 data. In the example below, the landsat_qa layer is the\nquality assurance data layer that comes with Landsat 8 to identify pixels that may represent\ncloud, shadow and water. The mask values used below are suggested values associated with the\nlandsat_qa layer that represent pixels with clouds and cloud shadows.\n\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Import Packages\n\nTo begin, import the needed packages. You will use a combination of several EarthPy\nmodules including spatial, plot and mask.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "from glob import glob\nimport os\nimport matplotlib.pyplot as plt\nimport rasterio as rio\nfrom rasterio.plot import plotting_extent\nimport earthpy as et\nimport earthpy.spatial as es\nimport earthpy.plot as ep\nimport earthpy.mask as em\n\n# Get data and set your home working directory\ndata = et.data.get_data(\"cold-springs-fire\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Import Example Data\nTo get started, make sure your directory is set. Create a stack from all of the\nLandsat .tif files (one per band) and import the ``landsat_qa`` layer which provides\nthe locations of cloudy and shadowed pixels in the scene.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "os.chdir(os.path.join(et.io.HOME, \"earth-analytics\"))\n\n# Stack the landsat bands\n# This creates a numpy array with each \"layer\" representing a single band\nlandsat_paths_pre = glob(\n \"data/cold-springs-fire/landsat_collect/LC080340322016070701T1-SC20180214145604/crop/*band*.tif\"\n)\nlandsat_paths_pre.sort()\narr_st, meta = es.stack(landsat_paths_pre)\n\n# Import the landsat qa layer\nwith rio.open(\n \"data/cold-springs-fire/landsat_collect/LC080340322016070701T1-SC20180214145604/crop/LC08_L1TP_034032_20160707_20170221_01_T1_pixel_qa_crop.tif\"\n) as landsat_pre_cl:\n landsat_qa = landsat_pre_cl.read(1)\n landsat_ext = plotting_extent(landsat_pre_cl)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Plot Histogram of Each Band in Your Data\nYou can view a histogram for each band in your dataset by using the\n``hist()`` function from the ``earthpy.plot`` module.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "ep.hist(arr_st)\nplt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Customize Histogram Plot with Titles and Colors\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "ep.hist(\n arr_st,\n colors=[\"blue\"],\n title=[\n \"Band 1\",\n \"Band 2\",\n \"Band 3\",\n \"Band 4\",\n \"Band 5\",\n \"Band 6\",\n \"Band 7\",\n ],\n)\nplt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## View Single Band Plots\nNext, have a look at the data, it looks like there is a large cloud that you\nmay want to mask out.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "ep.plot_bands(arr_st)\nplt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Mask the Data\nYou can use the EarthPy ``mask()`` function to handle this cloud.\nTo begin you need to have a layer that defines the pixels that\nyou wish to mask. In this case, the ``landsat_qa`` layer will be used.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "ep.plot_bands(\n landsat_qa,\n title=\"The Landsat QA Layer Comes with Landsat Data\\n It can be used to remove clouds and shadows\",\n)\nplt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Plot The Masked Data\nNow apply the mask and plot the masked data. The mask applies to every band in your data.\nThe mask values below are values documented in the Landsat 8 documentation that represent\nclouds and cloud shadows.\n\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": false }, "outputs": [], "source": [ "# Generate array of all possible cloud / shadow values\ncloud_shadow = [328, 392, 840, 904, 1350]\ncloud = [352, 368, 416, 432, 480, 864, 880, 928, 944, 992]\nhigh_confidence_cloud = [480, 992]\n\n# Mask the data\nall_masked_values = cloud_shadow + cloud + high_confidence_cloud\narr_ma = em.mask_pixels(arr_st, landsat_qa, vals=all_masked_values)\n\n# sphinx_gallery_thumbnail_number = 5\nep.plot_rgb(\n arr_ma, rgb=[4, 3, 2], title=\"Array with Clouds and Shadows Masked\"\n)\nplt.show()" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.7.9" } }, "nbformat": 4, "nbformat_minor": 0 }