PK vN7d " plot_calculate_classify_ndvi.ipynb{
"cells": [
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"%matplotlib inline"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\nCalculate and Classify Normalized Difference Results with EarthPy\n==================================================================\n\nLearn how to calculate and classify normalized difference indices in Python using EarthPy.\nThis example shows how to calculate and classify the normalized difference vegetation\nindex (NDVI) using Landsat 8 data.\n\n\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Calculating Normalized Difference in Python Using EarthPy\n---------------------------------------------------------\n\n
Note
The examples below will show you how to use the ``normalized_diff()`` function\n to calculate the normalized difference vegetation index (NDVI), a commonly\n used remotely sensed index for quantifying vegetation health.
\n\nThe example below walks you through a typical workflow for calculating the normalized\ndifference vegetation index (NDVI) using Landsat 8 data with EarthPy. NDVI provides\na measure of healthy vegetation and ranges in value from -1 to 1. Values closer to\n1 represent healthy, green vegetation. NDVI can be calculated from Landsat 8 data\nusing band 4 (red) and band 5 (near-infrared).\n\nFirst, you will create a stack of bands using Landsat 8 data and then calculate\nNDVI using the ``normalized_diff()`` function. Then, you will plot the NDVI results\nusing a colorbar legend with continuous values. Last, you will classify the NDVI\nresults using threshold values and plot the classified data with a categorical legend.\n\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Import Packages\n---------------\n\nTo begin, import the needed packages. You will use a combination of several EarthPy\nmodules including spatial and plot.\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"import os\nfrom glob import glob\nimport numpy as np\nimport matplotlib.pyplot as plt\nfrom matplotlib.colors import ListedColormap\nimport earthpy as et\nimport earthpy.spatial as es\nimport earthpy.plot as ep\n\n# Get data and set your home working directory\ndata = et.data.get_data(\"vignette-landsat\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Import Example Data\n-------------------\n\nTo get started, make sure your directory is set. Then, create a stack from all of\nthe Landsat .tif files (one per band). The nodata value for Landsat 8 is\n``-9999`` so you can use the ``nodata=`` parameter when you call the\n``stack()`` function.\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 8 bands\n# This creates a numpy array with each \"layer\" representing a single band\nlandsat_path = glob(\n \"data/vignette-landsat/LC08_L1TP_034032_20160621_20170221_01_T1_sr_band*_crop.tif\"\n)\nlandsat_path.sort()\narr_st, meta = es.stack(landsat_path, nodata=-9999)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Calculate Normalized Difference Vegetation Index (NDVI)\n-------------------------------------------------------\n\nYou can calculate NDVI for your dataset using the\n``normalized_diff`` function from the ``earthpy.spatial`` module.\nMath will be calculated (b1-b2) / (b1 + b2).\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"# Landsat 8 red band is band 4 at [3]\n# Landsat 8 near-infrared band is band 5 at [4]\nndvi = es.normalized_diff(arr_st[4], arr_st[3])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Plot NDVI With Colorbar Legend of Continuous Values\n----------------------------------------------------\n\nYou can plot NDVI with a colorbar legend of continuous values using the\n``plot_bands`` function from the ``earthpy.plot`` module.\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"titles = [\"Landsat 8 - Normalized Difference Vegetation Index (NDVI)\"]\n\n# Turn off bytescale scaling due to float values for NDVI\nep.plot_bands(\n ndvi, cmap=\"RdYlGn\", cols=1, title=titles, scale=False, vmin=-1, vmax=1\n)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Classify NDVI\n-------------\n\nNext, you can categorize (or classify) the NDVI results into useful classes.\nValues under 0 will be classified together as no vegetation. Additional classes\nwill be created for bare area and low, moderate, and high vegetation areas.\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"# Create classes and apply to NDVI results\nndvi_class_bins = [-np.inf, 0, 0.1, 0.25, 0.4, np.inf]\nndvi_landsat_class = np.digitize(ndvi, ndvi_class_bins)\n\n# Apply the nodata mask to the newly classified NDVI data\nndvi_landsat_class = np.ma.masked_where(np.ma.getmask(ndvi), ndvi_landsat_class)\nnp.unique(ndvi_landsat_class)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Plot Classified NDVI With Categorical Legend - EarthPy Draw_Legend()\n--------------------------------------------------------------------\n\nYou can plot the classified NDVI with a categorical legend using the\n``draw_legend()`` function from the ``earthpy.plot`` module.\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"# Define color map\nnbr_colors = [\"gray\", \"y\", \"yellowgreen\", \"g\", \"darkgreen\"]\nnbr_cmap = ListedColormap(nbr_colors)\n\n# Define class names\nndvi_cat_names = [\n \"No Vegetation\",\n \"Bare Area\",\n \"Low Vegetation\",\n \"Moderate Vegetation\",\n \"High Vegetation\",\n]\n\n# Get list of classes\nclasses = np.unique(ndvi_landsat_class)\nclasses = classes.tolist()\n# The mask returns a value of none in the classes. remove that\nclasses = classes[0:5]\n\n# Plot your data\nfig, ax = plt.subplots(figsize=(12, 12))\nim = ax.imshow(ndvi_landsat_class, cmap=nbr_cmap)\n\nep.draw_legend(im_ax=im, classes=classes, titles=ndvi_cat_names)\nax.set_title(\n \"Landsat 8 - Normalized Difference Vegetation Index (NDVI) Classes\",\n fontsize=14,\n)\nax.set_axis_off()\n\n# Auto adjust subplot to fit figure size\nplt.tight_layout()"
]
}
],
"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.3"
}
},
"nbformat": 4,
"nbformat_minor": 0
}PK kN plot_draw_legend_docs.ipynb{
"cells": [
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"%matplotlib inline"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\nCustom Legends with EarthPy\n===========================\n\nLearn how to create discrete legends for raster plots with classes that you define in Python using EarthPy.\n\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Plotting Data in Python Using EarthPy\n-------------------------------------\n\nNote
The examples below will show you how to use the ``draw_legend()`` function for creating plots
\n\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Plot Continuous Data\n---------------------\n\nLet's explore a simple plot using EarthPy. To begin, import the needed packages\nand create an array to be plotted. Below we plot the data as continuous with a colorbar\nusing the ``plot_bands()`` function.\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"import matplotlib.pyplot as plt\nimport earthpy.plot as ep\nimport numpy as np\n\n# Create a numpy array. Let's pretend this is what you want to plot.\narr = np.random.randint(4, size=(5, 5))\n\n# When plot_bands is updated a cbar will be here as well\nep.plot_bands(arr)\nplt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Create Custom Discrete Legends with Earthpy\n-------------------------------------------\nIf you want to create a custom categorical legend, you can use the ``ep.draw_legend()`` function.\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"f, ax = plt.subplots(figsize=(8, 5))\nim = ax.imshow(arr)\nep.draw_legend(im)\nplt.tight_layout()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Customize Discrete Legends\n---------------------------\nBy default the draw_legend function creates a legend with default categories.\nYou can customize the legend by adding titles.\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"f, ax = plt.subplots(figsize=(8, 5))\nim = ax.imshow(arr)\nep.draw_legend(im, titles=[\"Small\", \"Bigger\", \"Even Bigger\", \"Ginormous\"])\nplt.tight_layout()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Discrete Legends With Missing Values in the Array\n~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\n\nNow let's pretend that you are creating a series of classified plots. You may have a range of\nvalues that are expected between 0-4. However not all of your data has all values\nIn this case, your legend won't be able to by default create 4 categories because one\ndoesn't exist in your data. In this instance, you can specify the values explicitly:\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"new_arr = arr.copy()\nnew_arr[new_arr == 0] = 1\n\nf, ax = plt.subplots(figsize=(8, 5))\nim = ax.imshow(new_arr)\nep.draw_legend(\n im,\n titles=[\"Small\", \"Bigger\", \"Even Bigger\", \"Ginormous\"],\n classes=[0, 1, 2, 3],\n)\nplt.tight_layout()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Custom Colormaps and Ensuring Cmaps Apply to All Valid Classes\n~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\n\nYou can customize the color map used in your plot too. Notice that in this example,\n4 category colors are rendered. Yet, the image only contains three values and thus will be\nrendered using three colors. The three colors used to render the image are incorrect by default.\nThe colors begin at white and end at black even though the value of 0 which should be black\nis missing from the data.\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"f, ax = plt.subplots(figsize=(8, 5))\nim = ax.imshow(new_arr, cmap=\"Greys_r\")\nep.draw_legend(\n im,\n titles=[\"Small\", \"Bigger\", \"Even Bigger\", \"Ginormous\"],\n classes=[0, 1, 2, 3],\n)\nplt.tight_layout()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Specify vmin and vmax to set the colormap range\n~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~\n\nIn this case, you can use the ``vmin`` and ``vmax`` arguments to set the range of values to use for the colormap.\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"f, ax = plt.subplots(figsize=(8, 5))\nim = ax.imshow(new_arr, cmap=\"Greys_r\", vmin=0, vmax=3)\nep.draw_legend(\n im,\n titles=[\"Small\", \"Bigger\", \"Even Bigger\", \"Ginormous\"],\n classes=[0, 1, 2, 3],\n)\nplt.tight_layout()"
]
}
],
"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.3"
}
},
"nbformat": 4,
"nbformat_minor": 0
}PK |NZ1 1 plot_rgb.ipynb{
"cells": [
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"%matplotlib inline"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\nPlot Combinations of Raster Bands Using EarthPy\n===================================================\nLearn how to plot different band combinations from satellite imagery using EarthPy\n\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Create Composite Images Using EarthPy\n---------------------------------------------\n\nNote
The examples below will show you how to use the ``plot_rgb()``\n function from EarthPy to plot different band combinations and\n create composite images from satellite imagery such as Landsat.
\n\nA color image consists of three band channels displayed together: red, green,\nand blue. Some remote sensing datasets such as Landsat collect data across\nmany spectral bands and are stored with each band in a separate file. These\nbands can be combined to create true color composites (or RGB images) as well\nas color infrared (CIR) composites, which highlight vegetation in the image.\nTo create an RGB image, you need to plot the red, green, and blue bands,\nrespectively, such that red band is plotted on the red channel, green on the\ngreen, etc. In a color infrared image (sometimes called a false color image),\nthe near-infrared (NIR) band is plotted on the red channel. Thus,\nvegetation, which reflects strongly in the NIR part of the electromagnetic\nspectrum, will appear \"red\" in the resulting composite image.\n\nThe example below walks you through a typical workflow for plotting different\nband combinations to create RGB (true color) and CIR (false color) composite\nimages. You will create a stack of bands using Landsat 8 data and then plot\nthe different band combinations using the ``plot_rgb()`` function. You will\nalso learn how to stretch the image values to brighten dark images, how to\noverlay polygon boundaries on the images, and how to create figures that\ncontain multiple plots.\n\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Import Packages\n------------------------------\n\nTo begin, import the needed packages. You will use a combination of several\nEarthPy modules including spatial and plot.\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"import os\nfrom glob import glob\nimport matplotlib.pyplot as plt\nimport rasterio as rio\nfrom rasterio.plot import plotting_extent\nimport geopandas as gpd\nimport earthpy as et\nimport earthpy.spatial as es\nimport earthpy.plot as ep"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Import and Stack Example Data\n------------------------------\nTo plot band combinations and create composite images using ``plot_rgb()``,\nthe bands of a satellite image such as those from Landsat, need to be stacked\ninto one file. With EarthPy, you can create an image stack from all of the\nLandsat .tif files (one per band) in a directory using the ``stack()`` function\nfrom the ``earthpy.spatial`` module.\n\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Error found on Windows systems\n-------------------------------\nNote
If you are running this script on a Windows system, there is a\n known bug with ``.to_crs()``, which is used in this script. If an error\n occurs, you have to reset your os environment with the command\n ``os.environ[\"PROJ_LIB\"] = r\"path-to-share-folder-in-environment\"``.
\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"# Get sample data from EarthPy and set working directory\ndata_path = et.data.get_data(\"vignette-landsat\")\nos.chdir(os.path.join(et.io.HOME, \"earth-analytics\"))\n\n# Get list of bands and sort by ascending band number\nlandsat_bands_data_path = \"data/vignette-landsat/LC08_L1TP_034032_20160621_20170221_01_T1_sr_band*[1-7]*.tif\"\nstack_band_paths = glob(landsat_bands_data_path)\nstack_band_paths.sort()\n\n# Create image stack and apply nodata value for Landsat\narr_st, meta = es.stack(stack_band_paths, nodata=-9999)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Plot RGB Composite Image\n--------------------------\nYou can use the ``plot_rgb()`` function from the ``earthpy.plot`` module to quickly\nplot three band composite images. For RGB composite images, you will plot the red,\ngreen, and blue bands, which are bands 4, 3, and 2, respectively, in the image\nstack you created. Python uses a zero-based index system, so you need to subtract\na value of 1 from each index. Thus, the index for the red band is 3, green is 2,\nand blue is 1. These index values are provided to the ``rgb`` argument to identify\nthe bands for the composite image.\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"# Create figure with one plot\nfig, ax = plt.subplots(figsize=(12, 12))\n\n# Plot red, green, and blue bands, respectively\nep.plot_rgb(arr_st, rgb=(3, 2, 1), ax=ax, title=\"Landsat 8 RGB Image\")\nplt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Stretch Composite Images\n-------------------------\nComposite images can sometimes be dark if the pixel brightness values are\nskewed toward the value of zero. You can stretch the pixel brightness values\nin an image using the argument ``stretch=True`` to extend the values to the\nfull 0-255 range of potential values to increase the visual contrast of the\nimage. In addition, the ``str_clip`` argument allows you to specify how much of\nthe tails of the data that you want to clip off. The larger the number, the\nmore the data will be stretched or brightened.\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"# Create figure with one plot\nfig, ax = plt.subplots(figsize=(12, 12))\n\n# Plot bands with stretched applied\nep.plot_rgb(\n arr_st,\n rgb=(3, 2, 1),\n ax=ax,\n stretch=True,\n str_clip=0.5,\n title=\"Landsat 8 RGB Image with Stretch Applied\",\n)\nplt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Plot Color Infrared (CIR) Composite Image\n------------------------------------------\nFor color infrared (CIR) composite images, you will plot the near-infrared (NIR),\nred, and green bands, which are bands 5, 4, 2, respectively. Once again, the\nPython index values will be the original band number minus 1, thus, 4, 3, and 2\nfor NIR, red, and green, respectively.\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"# Create figure with one plot\nfig, ax = plt.subplots(figsize=(12, 12))\n\n# Plot NIR, red, and green bands, respectively, with stretch\nep.plot_rgb(\n arr_st,\n rgb=(4, 3, 2),\n ax=ax,\n stretch=True,\n str_clip=0.5,\n title=\"Landsat 8 CIR Image with Stretch Applied\",\n)\nplt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Plot Boundary Over Composite Image\n-----------------------------------\nNote
If you are on windows, you may need to add the crs issue discussed above\n here!
\n\nYou can overlay a polygon boundary on top of an RGB plot created with EarthPy.\nTo begin, the raster data and the boundary need to be in the same\nCoordinate Reference System (CRS). You can reproject the boundary layer to\nmatch the CRS of the image by getting the CRS of the image from the Rasterio\nprofile object and passing that CRS to the ``to_crs`` method from GeoPandas.\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# Open polygon boundary using GeoPandas\nbound = gpd.read_file(\n \"data/vignette-landsat/vector_layers/fire-boundary-geomac/co_cold_springs_20160711_2200_dd83.shp\"\n)\n\n# Reproject boundary to match CRS of the Landsat images\nwith rio.open(stack_band_paths[0]) as raster_crs:\n raster_profile = raster_crs.profile\n bound_utm13N = bound.to_crs(raster_profile[\"crs\"])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Create a Plot With the Boundary overlayed on the RGB Image \n----------------------------------------------------------\nYou can plot a polygon boundary over an image by creating a raster extent\nfor the plot using the ``plotting_extent`` function from ``rasterio.plot``.\nThe function needs the Rasterio profile of the image and a single layer of a\nnumpy array, which can be specified with ``arr_str[0]``. The function also\nneeds the spatial transformation for the Rasterio object, which can be acquired\nby accessing the ``\"transform\"`` key within the Rasterio profile.\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"# Create raster extent for the plot\nextent = plotting_extent(arr_st[0], raster_profile[\"transform\"])\n\n# Create figure with one plot\nfig, ax = plt.subplots(figsize=(12, 12))\n\n# Plot boundary with high zorder for contrast\nbound_utm13N.boundary.plot(ax=ax, color=\"black\", zorder=10)\n\n# Plot CIR image using the raster extent\nep.plot_rgb(\n arr_st,\n rgb=(4, 3, 2),\n ax=ax,\n stretch=True,\n extent=extent,\n str_clip=0.5,\n title=\"Landsat 8 CIR Image with Polygon Boundary\",\n)\nplt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Create Figure with Multiple Axes or Subplots\n--------------------------------------------\n```plot_rgb()`` has an ``ax=`` parameter which supports subplots. You can \ncreate figures that contain multiple plots by creating multiple ax\nobjects, one for each plot. You can also specify the number of rows and\ncolumns in which to display the plots. In the example below, the two plots\nwill be displayed side-by-side along one row with two columns.\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"# Create figure with two plots\nfig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 6))\n\n# Plot of RGB composite image with polygon boundary\nbound_utm13N.boundary.plot(ax=ax1, color=\"black\", zorder=10)\nax1 = ep.plot_rgb(\n arr_st,\n rgb=(3, 2, 1),\n ax=ax1,\n stretch=True,\n extent=extent,\n str_clip=0.5,\n title=\"Landsat 8 RGB Image with Polygon Boundary\",\n)\n\n# Plot of CIR composite image with polygon boundary\nbound_utm13N.boundary.plot(ax=ax2, color=\"black\", zorder=10)\nax2 = ep.plot_rgb(\n arr_st,\n rgb=(4, 3, 2),\n ax=ax2,\n stretch=True,\n extent=extent,\n str_clip=0.5,\n title=\"Landsat 8 CIR Image with Polygon Boundary\",\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.3"
}
},
"nbformat": 4,
"nbformat_minor": 0
}PK jN; plot_bands_functionality.ipynb{
"cells": [
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"%matplotlib inline"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\nPlot Bands of Satellite Imagery with EarthPy\n==================================================================\n\nLearn how to use the EarthPy ``plot_bands()`` function to quickly plot\nraster bands for an image. ``Plot_bands()`` can be used to both \nplot many bands with one command with custom titles and legends OR\nto plot a single raster layer with (or without) a legend.\n\n\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Plot Raster Data Layers Using EarthPy\n--------------------------------------\n\nNote
The examples below will show you how to use the ``plot_bands()`` function\n to plot individual raster layers in images using python. To plot rgb data,\n read help documentation related to ``ep.plot_rgb()``.
\n\nIn this vignette, you will use Landsat 8 data. To begin, you will create a\nstack of bands using Landsat 8 data. You will then plot the raster layers.\n\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Import Packages\n---------------\n\nIn order to use the ``plot_bands()`` function with Landsat 8 data, the\nfollowing packages need to be imported.\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"import os\nfrom glob import glob\nimport matplotlib.pyplot as plt\nimport earthpy as et\nimport earthpy.spatial as es\nimport earthpy.plot as ep"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Import Example Data\n-------------------\n\nTo get started, make sure your directory is set. Then, create a stack from all of\nthe Landsat .tif files (one per band).\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"# Get data for example\ndata = et.data.get_data(\"vignette-landsat\")\n\n# Set working directory\nos.chdir(os.path.join(et.io.HOME, \"earth-analytics\"))\n\n# Stack the Landsat 8 bands\n# This creates a numpy array with each \"layer\" representing a single band\n# You can use the nodata= parameter to mask nodata values\nlandsat_path = glob(\n \"data/vignette-landsat/LC08_L1TP_034032_20160621_20170221_01_T1_sr_band*_crop.tif\"\n)\nlandsat_path.sort()\narray_stack, meta_data = es.stack(landsat_path, nodata=-9999)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Plot All Bands in a Stack\n--------------------------\n\nWhen you give ``ep.plot_bands()`` a three dimensional numpy array,\nit will plot all layers in the numpy array. You can create unique titles for\neach image by providing a list of titles using the ``title=`` parameter.\nThe list must contain the same number of strings as there are bands in the stack.\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"titles = [\"Ultra Blue\", \"Blue\", \"Green\", \"Red\", \"NIR\", \"SWIR 1\", \"SWIR 2\"]\n# sphinx_gallery_thumbnail_number = 5\nep.plot_bands(array_stack, title=titles)\nplt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Plot One Band in a Stack\n------------------------\n\nIf you give ``ep.plot_bands()`` a one dimensional numpy array,\nit will only plot that single band. You can turn off the\ncolorbar using the ``cbar`` parameter (``cbar=False``).\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"ep.plot_bands(array_stack[4], cbar=False)\nplt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Turn Off Scaling\n-----------------\n\n``ep.plot_bands()`` scales the imagery to a 0-255 scale by default. This range\nof values makes it easier for matplotlib to plot the data. To turn off\nscaling, set the scale parameter to ``False``. Below you\nplot NDVI with scaling turned off in order for the proper range of values\n(-1 to 1) to be displayed. You can use the ``cmap=`` parameter to adjust\nthe colormap for the plot\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"NDVI = es.normalized_diff(array_stack[4], array_stack[3])\nep.plot_bands(NDVI, scale=False, cmap=\"RdYlGn\")\nplt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Adjust the Number of Columns for a Multi Band Plot\n---------------------------------------------------\n\nThe number of columns used while plotting multiple bands can be changed in order\nto change the arrangement of the images overall.\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"ep.plot_bands(array_stack, cols=2)\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.3"
}
},
"nbformat": 4,
"nbformat_minor": 0
}PK qN^w徹 plot_stack_masks.ipynb{
"cells": [
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"%matplotlib inline"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\nMask and Plot Remote Sensing Data with EarthPy\n==============================================\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\n\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Plotting with EarthPy\n---------------------\n\nNote
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------------------------------\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\n------------------------------\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\n----------------------------------------\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\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\n-----------------------------------------------\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\n-----------------------------------------------\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\n~~~~~~~~~~~~~~~~~~~~~\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.3"
}
},
"nbformat": 4,
"nbformat_minor": 0
}PK lNݍ{} plot_hist_functionality.ipynb{
"cells": [
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"%matplotlib inline"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\nPlot Histograms of Pixel Values from Multi-band Imagery with EarthPy\n====================================================================\nLearn how to quickly plot distributions of pixel values in Python\nusing Earthpy. This examples shows you how to create histogram\nplots for each raster band in a multi-band image such as\nLandsat 8 data.\n\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Plot Raster Data Histograms Using EarthPy\n------------------------------------------\n\nNote
This example shows you how to create histogram plots of pixel values\n for each raster band of a multi-band image using the ``ep.hist()``\n function from the ``earthpy.spatial`` module.
\n\nIn this example, you will learn how to plot histograms from multi-band\nimagery such as Landsat 8. Multi-band images store data as individual raster\nbands that contain reflectance values for various sections of the\nelectromagnetic spectrum. For example, the second band of Landsat 8 provides\nsurface reflectance within the blue wavelength of the electromagnetic\nspectrum, while the third and fourth bands of Landsat 8 provides surface\nreflectance within the green and red wavelengths, respectively.\nHistogram plots provide a quick way to review the distribution of pixel\nvalues for a specific band, which can be helpful to identify potential\ndata quality issues or simply learn more about the surface area that was\ncaptured.\n\nTo begin, you will create a stack of these individuals bands from\nLandsat 8 data and then use the ``ep.hist()`` function to plot the histograms\nfor each band in the stack.\n\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Import Packages\n----------------\n\nIn order to use the ``ep.hist()`` function from EarthPy, the following\npackages need to be imported.\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"import os\nfrom glob import glob\nimport matplotlib.pyplot as plt\nimport earthpy as et\nimport earthpy.spatial as es\nimport earthpy.plot as ep"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Import Example Data\n-------------------\n\nTo get started, make sure your directory is set. Then, create a stack from all\nof the Landsat .tif files (one per band).\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"# Get data for example\ndata = et.data.get_data(\"vignette-landsat\")\n\n# Set working directory\nos.chdir(os.path.join(et.io.HOME, \"earth-analytics\"))\n\n# Stack the Landsat 8 bands\n# This creates a numpy array with each \"layer\" representing a single band\n# You can use the nodata parameter to mask nodata values\nlandsat_path = glob(os.path.join(\"data\", \"vignette-landsat\", \"*band*.tif\"))\nlandsat_path.sort()\narray_stack, meta_data = es.stack(landsat_path, nodata=-9999)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Plot All Histograms in a Stack With Custom Titles and Colors\n-------------------------------------------------------------\n\nYou can create histograms for each band with unique colors and titles by\nfirst creating a list of colors and titles that will be provided to the\nep.hist() function. The list for titles must contain the same number of\nstrings as there are bands in the stack. You can also modify the colors for\neach image with a list of Matplotlib colors. If one color is provided in\nthe list, every band will inherit that color. Otherwise, the list must\ncontain the same number of strings as there are bands in the stack.\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"# Create the list of color names for each band\ncolors_list = [\n \"midnightblue\",\n \"Blue\",\n \"Green\",\n \"Red\",\n \"Maroon\",\n \"Purple\",\n \"Violet\",\n]\n\n# Create the list of titles for each band. The titles and colors listed\n# in this example reflect the order and wavelengths of the Landsat 8 bands\ntitles = [\"Ultra Blue\", \"Blue\", \"Green\", \"Red\", \"NIR\", \"SWIR 1\", \"SWIR 2\"]\n\n# Plot the histograms with the color and title lists you just created\n# sphinx_gallery_thumbnail_number = 5\nep.hist(array_stack, colors=colors_list, title=titles)\nplt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Customize Bin Size and Arrangement of Histograms\n-------------------------------------------------\n\nYou can customize the number of bins each histogram plot uses to group the\ndata it is plotting. The default number is 20. This can be adjusted to match\nthe data you are trying to display. Additionally, you can change the\narrangement of the image overall by modifying the number of columns used\nto plot the data.\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"# Plot each histogram with 50 bins, arranged across three columns\nep.hist(array_stack, bins=50, cols=3)\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.3"
}
},
"nbformat": 4,
"nbformat_minor": 0
}PK Nv)8 8 plot_raster_stack_crop.ipynb{
"cells": [
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"%matplotlib inline"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\nStack and Crop Raster Data Using EarthPy\n========================================\n\nLearn how to stack and crop satellite imagery using EarthPy\n\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Stack and Crop Raster Data Using EarthPy\n---------------------------------------------\n\nNote
The examples below will show you how to use the ``es.stack()`` and\n ``es.crop_image()`` functions from EarthPy.
\n\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Stack Multi Band Imagery\n-----------------------------\nSome remote sensing datasets are stored with each band in a separate file. However,\noften you want to use all of the bands together in your analysis. For example\nyou need all of the bands together in the same file or \"stack\" in order to plot a color\nRGB image. EarthPy has a ``stack()`` function that allows you\nto take a set of ``.tif`` files that are all in the same spatial extent, CRS and resolution\nand either export them together a single stacked ``.tif`` file or work with them in Python\ndirectly as a stacked numpy array.\n\nTo begin using the EarthPy ``stack()`` function, import the needed packages\nand create an array to be plotted. Below you plot the data as continuous with a colorbar\nusing the ``plot_bands()`` function.\n\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Import Packages\n------------------------------\n\nYou will need several packages to stack your raster. You will use GeoPandas to\nopen up a shapefile that will be used to crop your data. You will primarily be\nusing the EarthPy spatial module in this vignette.\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"import os\nfrom glob import glob\nimport matplotlib.pyplot as plt\nimport rasterio as rio\nfrom rasterio.plot import plotting_extent\nimport geopandas as gpd\nimport earthpy as et\nimport earthpy.spatial as es\nimport earthpy.plot as ep"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Get Example Data Ready for Stack\n----------------------------------\nWith EarthPy you can create a stack from all of the Landsat .tif files (one per band)\nin a folder with the ``es.stack()`` function.\n\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Error found on Windows systems\n-------------------------------\nNote
If you are running this script on a Windows system, there is a\n known bug with ``.to_crs()``, which is used in this script. If an error\n occurs, you have to reset your os environment with the command\n ``os.environ[\"PROJ_LIB\"] = r\"path-to-share-folder-in-environment\"``.
\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"# Get sample data from EarthPy and setting your home working directory\n\ndata_path = et.data.get_data(\"vignette-landsat\")\nos.chdir(os.path.join(et.io.HOME, \"earth-analytics\"))\n\n# Prepare the landsat bands to be stacked using glob and sort\n\nlandsat_bands_data_path = \"data/vignette-landsat/LC08_L1TP_034032_20160621_20170221_01_T1_sr_band*[2-4]*_crop.tif\"\nstack_band_paths = glob(landsat_bands_data_path)\nstack_band_paths.sort()\n\n# Create output directory and the output path\n\noutput_dir = os.path.join(\"data\", \"outputs\")\nif os.path.isdir(output_dir) == False:\n os.mkdir(output_dir)\n\nraster_out_path = os.path.join(output_dir, \"raster.tiff\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Stack the Bands\n---------------------------\nThe stack function has an optional output argument, where you can write the raster\nto a tiff file in a folder. If you want to use this functionality, make sure there\nis a folder to write your tiff file to.\nThe Stack function also returns two object, an array and a RasterIO profile. Make\nsure to be catch both in variables.\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"# Stack Landsat bands\n\nos.chdir(os.path.join(et.io.HOME, \"earth-analytics\"))\narray, raster_prof = es.stack(stack_band_paths, out_path=raster_out_path)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Create Extent Object\n--------------------------------\nTo get the raster extent, use the ``plotting_extent`` function on the\narray from ``es.stack()`` and the Rasterio profile or metadata object. The function\nneeds a single\nlayer of a numpy array, which is why we use ``arr[0]``. The function also\nneeds the spatial transformation for the Rasterio object, which can be acquired by accessing\nthe ``\"transform\"`` key within the Rasterio Profile.\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"extent = plotting_extent(array[0], raster_prof[\"transform\"])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Plot Un-cropped Data\n------------------------------\nYou can see the boundary and the raster before the crop using ``ep.plot_rgb()``\nNotice that the data appear washed out.\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"fig, ax = plt.subplots(figsize=(12, 12))\nep.plot_rgb(\n array,\n ax=ax,\n stretch=True,\n extent=extent,\n str_clip=0.5,\n title=\"RGB Image of Un-cropped Raster\",\n)\nplt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Explore the Range of Values in the Data\n---------------------------------------\nYou can explore the range of values found in the data using the EarthPy ``hist()``\nfunction. Do you notice any extreme values that may be impacting the stretch\nof the image?\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"ep.hist(array, title=[\"Band 1\", \"Band 2\", \"Band 3\"])\nplt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"No Data Option\n---------------\n``es.stack()`` can handle ``nodata`` values in a raster. To use this\nparameter, specify ``nodata=``. This will mask every pixel that contains\nthe specified ``nodata`` value. The output will be a numpy masked array.\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"os.chdir(os.path.join(et.io.HOME, \"earth-analytics\"))\narray_nodata, raster_prof_nodata = es.stack(stack_band_paths, nodata=-9999)\n\n# View hist of data with nodata values removed\nep.hist(\n array_nodata,\n title=[\n \"Band 1 - No Data Values Removed\",\n \"Band 2 - No Data Values Removed\",\n \"Band 3 - No Data Values Removed\",\n ],\n)\nplt.show()\n\n# Recreate extent object for the No Data array\n\nextent_nodata = plotting_extent(\n array_nodata[0], raster_prof_nodata[\"transform\"]\n)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Plot Un-cropped Data\n------------------------------\nPlot the data again after the nodata values are removed.\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"fig, ax = plt.subplots(figsize=(12, 12))\nep.plot_rgb(\n array_nodata,\n ax=ax,\n stretch=True,\n extent=extent,\n str_clip=0.5,\n title=\"RGB image of Un-cropped Raster, No Data Value Selected\",\n)\nplt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Crop the Data\n------------------\nSometimes you have data for an area that is larger than your study area.\nIt is more efficient to first crop the data to your study area before processing\nit in Python. The fastest and most efficient option is to crop each file\nindividually, write out the cropped raster to a new file, and then stack\nthe new files together. To do this, make sure you have a ShapeFile boundary\nin the form of a GeoPandas object you can use as the cropping object.\nThen, loop through every file you wish to crop and crop the image, then\nwrite it out to a file. Take the rasters created and stack them like\nyou stacked bands in the previous examples.\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# Open the crop boundary using GeoPandas.\n\ncrop_bound = gpd.read_file(\n \"data/vignette-landsat/vector_layers/fire-boundary-geomac/co_cold_springs_20160711_2200_dd83.shp\"\n)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Reproject the data\n------------------\nNote
If you are on windows, make sure to set your environment here!
\n\nThe crop function won't work properly if the data are in different Coordinate\nReference Systems (CRS). To fix this, be sure to reproject the crop layer to match\nthe CRS of your raster data.\nTo reproject your data, first get the CRS of the raster from the rasterio profile\nobject. Then use that to reproject using geopandas ``.to_crs`` method.\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"os.chdir(os.path.join(et.io.HOME, \"earth-analytics\"))\n\nwith rio.open(stack_band_paths[0]) as raster_crs:\n crop_raster_profile = raster_crs.profile\n crop_bound_utm13N = crop_bound.to_crs(crop_raster_profile[\"crs\"])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Crop Each Band\n--------------\nWhen you need to crop and stack a set of images, it is most efficient to first\ncrop each image, and then stack it.\n``es.crop_all()`` is an efficient way to crop all bands in an image quickly.\nThe function will write out cropped rasters to a\ndirectory and return a list of file paths that can then be used with\n``es.stack()``.\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"os.chdir(os.path.join(et.io.HOME, \"earth-analytics\"))\n\nband_paths_list = es.crop_all(\n stack_band_paths, output_dir, crop_bound_utm13N, overwrite=True\n)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Stack All Bands\n---------------\nOnce the data are cropped, you are ready to create a new stack.\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"os.chdir(os.path.join(et.io.HOME, \"earth-analytics\"))\n\ncropped_array, array_raster_profile = es.stack(band_paths_list, nodata=-9999)\ncrop_extent = plotting_extent(\n cropped_array[0], array_raster_profile[\"transform\"]\n)\n\n# Plotting the cropped image\n# sphinx_gallery_thumbnail_number = 5\nfig, ax = plt.subplots(figsize=(12, 6))\ncrop_bound_utm13N.boundary.plot(ax=ax, color=\"red\", zorder=10)\nep.plot_rgb(\n cropped_array,\n ax=ax,\n stretch=True,\n extent=crop_extent,\n title=\"Cropped Raster and Fire Boundary\",\n)\nplt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Crop Individual Bands\n---------------------\nIf you only need to crop one raster image, you can use EarthPy's \n``es.crop_image()`` function.\nThis function takes a Rasterio object and crops it to the provided \nspatial extent.\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"# Open Landsat image as a Rasterio object in order to crop it\nos.chdir(os.path.join(et.io.HOME, \"earth-analytics\"))\n\nwith rio.open(stack_band_paths[0]) as src:\n single_cropped_image, single_cropped_meta = es.crop_image(\n src, crop_bound_utm13N\n )\n\n# Create the extent object\nsingle_crop_extent = plotting_extent(\n single_cropped_image[0], single_cropped_meta[\"transform\"]\n)\n\n# Plot the newly cropped image\nfig, ax = plt.subplots(figsize=(12, 6))\ncrop_bound_utm13N.boundary.plot(ax=ax, color=\"red\", zorder=10)\nep.plot_bands(\n single_cropped_image,\n ax=ax,\n extent=single_crop_extent,\n title=\"Single Cropped Raster and Fire Boundary\",\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.3"
}
},
"nbformat": 4,
"nbformat_minor": 0
}PK uNz:|ψ plot_dem_hillshade.ipynb{
"cells": [
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"%matplotlib inline"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\nCreate A Hillshade From a Digital Elevation Model Using EarthPy\n===============================================================\n\nLearn how to create a hillshade from a DEM using the EarthPy\nes.hillshade() function.\n\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Create A Hillshade Layer Using EarthPy\n---------------------------------------\n\nNote
The examples below show you how to use the ``es.hillshade()``\n function from EarthPy.
\n\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Create a Hillshade from a Digital Elevation Model (DEM)\n-------------------------------------------------------\nA hillshade is a 3D representation of a surface. Hillshades are generally\nrendered in greyscale. The darker and lighter colors represent the shadows\nand highlights that you would visually expect to see in a terrain model.\nHillshades are often used as an underlay in a map, to make the data appear\nmore 3-Dimensional and thus visually interesting. This vignette will show\nyou how to create a hillshade from a DEM using EarthPy. It will highlight\nhow to adjust the sun's azimuth, altitude and other settings that will impact\nhow the hillshade shadows are modeled in the data.\n\nThe hillshade function is a part of the spatial module in EarthPy.\n\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Import Packages\n----------------\n\nYou will need several packages to get started. You will use Rasterio to\nopen up the DEM file needed to create the hillshade layer.\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"import os\nimport numpy as np\nimport matplotlib.pyplot as plt\nimport earthpy as et\nimport earthpy.spatial as es\nimport earthpy.plot as ep\nimport rasterio as rio\n\n# Download the data needed for this vignette\ndata = et.data.get_data(\"vignette-elevation\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Open up the DEM\n---------------\nTo begin, open your DEM layer as a numpy array using Rasterio. Below you set all\nterrain values < 0 to ``nan``. Next, plot the data using ``ep.plot_bands()``.\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"# Set the home directory and get the data for the exercise\nos.chdir(os.path.join(et.io.HOME, \"earth-analytics\"))\ndtm = \"data/vignette-elevation/pre_DTM.tif\"\n\n# Open the DEM with Rasterio\nwith rio.open(dtm) as src:\n elevation = src.read(1)\n # Set masked values to np.nan\n elevation[elevation < 0] = np.nan\n\n# Plot the data\nep.plot_bands(\n elevation,\n scale=False,\n cmap=\"gist_earth\",\n title=\"DTM Without Hillshade\",\n figsize=(10, 6),\n)\nplt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Create the Hillshade\n--------------------\nOnce the DEM is read in, call ``es.hillshade()`` to create the hillshade.\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"# Create and plot the hillshade with earthpy\nhillshade = es.hillshade(elevation)\n\nep.plot_bands(\n hillshade,\n scale=False,\n cbar=False,\n title=\"Hillshade made from DTM\",\n figsize=(10, 6),\n)\nplt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Change the Azimuth of the Sun\n-------------------------------\nThe angle that sun light hits the landscape, impacts the shadows and highlights\ncreated on the landscape. You can adjust the azimuth values to adjust angle of the\nhighlights and shadows that are created in your output hillshade. Azimuth numbers can\nrange from 0 to 360 degrees, where 0 is due North. The default value for azimuth\nin ``es.hillshade()`` is 30 degrees.\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"# Change the azimuth of the hillshade layer\nhillshade_azimuth_210 = es.hillshade(elevation, azimuth=210)\n\n# Plot the hillshade layer with the modified azimuth\nep.plot_bands(\n hillshade_azimuth_210,\n scale=False,\n cbar=False,\n title=\"Hillshade with Azimuth set to 210 Degrees\",\n figsize=(10, 6),\n)\nplt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Change the Angle Altitude of the Sun\n-------------------------------------\nAnother variable you can adjust for hillshade is what angle of the sun.\nThe ``angle_altitude`` parameter values range from 0 to 90. 90 represents the sun\nshining from directly above the scene. The default value for ``angle_altitude`` in\n``es.hillshade()`` is 30 degrees.\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"# Adjust the azimuth value\nhillshade_angle_10 = es.hillshade(elevation, angle_altitude=10)\n\n# Plot the hillshade layer with the modified angle altitude\nep.plot_bands(\n hillshade_angle_10,\n scale=False,\n cbar=False,\n title=\"Hillshade with Angle Altitude set to 10 Degrees\",\n figsize=(10, 6),\n)\nplt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Overlay a DEM on top of the Hillshade\n-------------------------------------\nA hillshade can be used to visually enhance a DEM.\nTo overlay the data, use the ``ep.plot_bands()`` function in EarthPy combined with\n``ax.imshow()``. The alpha setting sets the tranparency value for the hillshade layer.\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"# Plot the DEM and hillshade at the same time\n# sphinx_gallery_thumbnail_number = 5\nfig, ax = plt.subplots(figsize=(10, 6))\nep.plot_bands(\n elevation,\n ax=ax,\n scale=False,\n cmap=\"terrain\",\n title=\"Lidar Digital Elevation Model (DEM)\\n overlayed on top of a hillshade\",\n)\nax.imshow(hillshade, cmap=\"Greys\", alpha=0.5)\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.3"
}
},
"nbformat": 4,
"nbformat_minor": 0
}PK vN7d " plot_calculate_classify_ndvi.ipynbPK kN plot_draw_legend_docs.ipynbPK |NZ1 1 27 plot_rgb.ipynbPK jN; ih plot_bands_functionality.ipynbPK qN^w徹 À plot_stack_masks.ipynbPK lNݍ{} plot_hist_functionality.ipynbPK Nv)8 8 ݵ plot_raster_stack_crop.ipynbPK uNz:|ψ plot_dem_hillshade.ipynbPK @