The Grizzly Creek Fire, Colorado, USA¶
In [4]:
# Install the development version of the earthpy package
!pip install git+https://github.com/earthlab/earthpy@apppears
Collecting git+https://github.com/earthlab/earthpy@apppears Cloning https://github.com/earthlab/earthpy (to revision apppears) to /tmp/pip-req-build-co0_z09e Running command git clone --filter=blob:none --quiet https://github.com/earthlab/earthpy /tmp/pip-req-build-co0_z09e Running command git checkout -b apppears --track origin/apppears Switched to a new branch 'apppears' Branch 'apppears' set up to track remote branch 'apppears' from 'origin'. Resolved https://github.com/earthlab/earthpy to commit 7241165d59af510d62bba312e48c7f513bc9dc05 Installing build dependencies ... done Getting requirements to build wheel ... done Preparing metadata (pyproject.toml) ... done Requirement already satisfied: geopandas in /opt/conda/lib/python3.11/site-packages (from earthpy==0.10.0) (0.14.2) Requirement already satisfied: matplotlib>=2.0.0 in /opt/conda/lib/python3.11/site-packages (from earthpy==0.10.0) (3.8.4) Requirement already satisfied: numpy>=1.14.0 in /opt/conda/lib/python3.11/site-packages (from earthpy==0.10.0) (1.24.3) Requirement already satisfied: rasterio in /opt/conda/lib/python3.11/site-packages (from earthpy==0.10.0) (1.3.9) Requirement already satisfied: scikit-image in /opt/conda/lib/python3.11/site-packages (from earthpy==0.10.0) (0.22.0) Requirement already satisfied: requests in /opt/conda/lib/python3.11/site-packages (from earthpy==0.10.0) (2.31.0) Requirement already satisfied: keyring in /opt/conda/lib/python3.11/site-packages (from earthpy==0.10.0) (25.1.0) Requirement already satisfied: contourpy>=1.0.1 in /opt/conda/lib/python3.11/site-packages (from matplotlib>=2.0.0->earthpy==0.10.0) (1.2.0) Requirement already satisfied: cycler>=0.10 in /opt/conda/lib/python3.11/site-packages (from matplotlib>=2.0.0->earthpy==0.10.0) (0.11.0) Requirement already satisfied: fonttools>=4.22.0 in /opt/conda/lib/python3.11/site-packages (from matplotlib>=2.0.0->earthpy==0.10.0) (4.51.0) Requirement already satisfied: kiwisolver>=1.3.1 in /opt/conda/lib/python3.11/site-packages (from matplotlib>=2.0.0->earthpy==0.10.0) (1.4.4) Requirement already satisfied: packaging>=20.0 in /opt/conda/lib/python3.11/site-packages (from matplotlib>=2.0.0->earthpy==0.10.0) (24.0) Requirement already satisfied: pillow>=8 in /opt/conda/lib/python3.11/site-packages (from matplotlib>=2.0.0->earthpy==0.10.0) (10.3.0) Requirement already satisfied: pyparsing>=2.3.1 in /opt/conda/lib/python3.11/site-packages (from matplotlib>=2.0.0->earthpy==0.10.0) (3.0.9) Requirement already satisfied: python-dateutil>=2.7 in /opt/conda/lib/python3.11/site-packages (from matplotlib>=2.0.0->earthpy==0.10.0) (2.9.0) Requirement already satisfied: fiona>=1.8.21 in /opt/conda/lib/python3.11/site-packages (from geopandas->earthpy==0.10.0) (1.9.6) Requirement already satisfied: pandas>=1.4.0 in /opt/conda/lib/python3.11/site-packages (from geopandas->earthpy==0.10.0) (2.2.1) Requirement already satisfied: pyproj>=3.3.0 in /opt/conda/lib/python3.11/site-packages (from geopandas->earthpy==0.10.0) (3.6.1) Requirement already satisfied: shapely>=1.8.0 in /opt/conda/lib/python3.11/site-packages (from geopandas->earthpy==0.10.0) (2.0.4) Requirement already satisfied: jaraco.classes in /opt/conda/lib/python3.11/site-packages (from keyring->earthpy==0.10.0) (3.4.0) Requirement already satisfied: jaraco.functools in /opt/conda/lib/python3.11/site-packages (from keyring->earthpy==0.10.0) (4.0.1) Requirement already satisfied: jaraco.context in /opt/conda/lib/python3.11/site-packages (from keyring->earthpy==0.10.0) (5.3.0) Requirement already satisfied: importlib-metadata>=4.11.4 in /opt/conda/lib/python3.11/site-packages (from keyring->earthpy==0.10.0) (7.1.0) Requirement already satisfied: SecretStorage>=3.2 in /opt/conda/lib/python3.11/site-packages (from keyring->earthpy==0.10.0) (3.3.3) Requirement already satisfied: jeepney>=0.4.2 in /opt/conda/lib/python3.11/site-packages (from keyring->earthpy==0.10.0) (0.8.0) Requirement already satisfied: affine in /opt/conda/lib/python3.11/site-packages (from rasterio->earthpy==0.10.0) (2.3.0) Requirement already satisfied: attrs in /opt/conda/lib/python3.11/site-packages (from rasterio->earthpy==0.10.0) (23.2.0) Requirement already satisfied: certifi in /opt/conda/lib/python3.11/site-packages (from rasterio->earthpy==0.10.0) (2024.2.2) Requirement already satisfied: click>=4.0 in /opt/conda/lib/python3.11/site-packages (from rasterio->earthpy==0.10.0) (8.1.7) Requirement already satisfied: cligj>=0.5 in /opt/conda/lib/python3.11/site-packages (from rasterio->earthpy==0.10.0) (0.7.2) Requirement already satisfied: snuggs>=1.4.1 in /opt/conda/lib/python3.11/site-packages (from rasterio->earthpy==0.10.0) (1.4.7) Requirement already satisfied: click-plugins in /opt/conda/lib/python3.11/site-packages (from rasterio->earthpy==0.10.0) (1.1.1) Requirement already satisfied: setuptools in /opt/conda/lib/python3.11/site-packages (from rasterio->earthpy==0.10.0) (69.5.1) Requirement already satisfied: charset-normalizer<4,>=2 in /opt/conda/lib/python3.11/site-packages (from requests->earthpy==0.10.0) (3.3.2) Requirement already satisfied: idna<4,>=2.5 in /opt/conda/lib/python3.11/site-packages (from requests->earthpy==0.10.0) (3.7) Requirement already satisfied: urllib3<3,>=1.21.1 in /opt/conda/lib/python3.11/site-packages (from requests->earthpy==0.10.0) (2.2.1) Requirement already satisfied: scipy>=1.8 in /opt/conda/lib/python3.11/site-packages (from scikit-image->earthpy==0.10.0) (1.12.0) Requirement already satisfied: networkx>=2.8 in /opt/conda/lib/python3.11/site-packages (from scikit-image->earthpy==0.10.0) (3.1) Requirement already satisfied: imageio>=2.27 in /opt/conda/lib/python3.11/site-packages (from scikit-image->earthpy==0.10.0) (2.33.1) Requirement already satisfied: tifffile>=2022.8.12 in /opt/conda/lib/python3.11/site-packages (from scikit-image->earthpy==0.10.0) (2023.4.12) Requirement already satisfied: lazy_loader>=0.3 in /opt/conda/lib/python3.11/site-packages (from scikit-image->earthpy==0.10.0) (0.3) Requirement already satisfied: six in /opt/conda/lib/python3.11/site-packages (from fiona>=1.8.21->geopandas->earthpy==0.10.0) (1.16.0) Requirement already satisfied: zipp>=0.5 in /opt/conda/lib/python3.11/site-packages (from importlib-metadata>=4.11.4->keyring->earthpy==0.10.0) (3.17.0) Requirement already satisfied: pytz>=2020.1 in /opt/conda/lib/python3.11/site-packages (from pandas>=1.4.0->geopandas->earthpy==0.10.0) (2024.1) Requirement already satisfied: tzdata>=2022.7 in /opt/conda/lib/python3.11/site-packages (from pandas>=1.4.0->geopandas->earthpy==0.10.0) (2023.3) Requirement already satisfied: cryptography>=2.0 in /opt/conda/lib/python3.11/site-packages (from SecretStorage>=3.2->keyring->earthpy==0.10.0) (42.0.5) Requirement already satisfied: more-itertools in /opt/conda/lib/python3.11/site-packages (from jaraco.classes->keyring->earthpy==0.10.0) (10.2.0) Requirement already satisfied: backports.tarfile in /opt/conda/lib/python3.11/site-packages (from jaraco.context->keyring->earthpy==0.10.0) (1.1.1) Requirement already satisfied: cffi>=1.12 in /opt/conda/lib/python3.11/site-packages (from cryptography>=2.0->SecretStorage>=3.2->keyring->earthpy==0.10.0) (1.16.0) Requirement already satisfied: pycparser in /opt/conda/lib/python3.11/site-packages (from cffi>=1.12->cryptography>=2.0->SecretStorage>=3.2->keyring->earthpy==0.10.0) (2.22)
In [5]:
import getpass
import json
import os
import pathlib
from glob import glob
import earthpy.appeears as eaapp
import hvplot.pandas
import hvplot.xarray
import rioxarray as rxr
import xarray as xr
import geopandas as gpd
import pandas as pd
In [6]:
data_dir = os.path.join(pathlib.Path.home(), 'grizzly-creek-fire')
# Make the data directory
os.makedirs(data_dir, exist_ok=True)
In [7]:
# Download the Cameron Peak fire boundary
fire_url_gdf = gpd.read_file(
"https://services3.arcgis.com/T4QMspbfLg3qTGWY/arcgis/rest/services/"
"WFIGS_Interagency_Perimeters/FeatureServer/0/query?where"
"=poly_IncidentName%20%3D%20'GRIZZLY%20CREEK'&outFields"
"=*&outSR=4326&f=json"
)
fire_url_gdf
/opt/conda/lib/python3.11/site-packages/geopandas/io/file.py:399: FutureWarning: errors='ignore' is deprecated and will raise in a future version. Use to_datetime without passing `errors` and catch exceptions explicitly instead as_dt = pd.to_datetime(df[k], errors="ignore") /opt/conda/lib/python3.11/site-packages/geopandas/io/file.py:399: FutureWarning: errors='ignore' is deprecated and will raise in a future version. Use to_datetime without passing `errors` and catch exceptions explicitly instead as_dt = pd.to_datetime(df[k], errors="ignore") /opt/conda/lib/python3.11/site-packages/geopandas/io/file.py:399: FutureWarning: errors='ignore' is deprecated and will raise in a future version. Use to_datetime without passing `errors` and catch exceptions explicitly instead as_dt = pd.to_datetime(df[k], errors="ignore") /opt/conda/lib/python3.11/site-packages/geopandas/io/file.py:399: FutureWarning: errors='ignore' is deprecated and will raise in a future version. Use to_datetime without passing `errors` and catch exceptions explicitly instead as_dt = pd.to_datetime(df[k], errors="ignore") /opt/conda/lib/python3.11/site-packages/geopandas/io/file.py:399: FutureWarning: errors='ignore' is deprecated and will raise in a future version. Use to_datetime without passing `errors` and catch exceptions explicitly instead as_dt = pd.to_datetime(df[k], errors="ignore") /opt/conda/lib/python3.11/site-packages/geopandas/io/file.py:399: FutureWarning: errors='ignore' is deprecated and will raise in a future version. Use to_datetime without passing `errors` and catch exceptions explicitly instead as_dt = pd.to_datetime(df[k], errors="ignore") /opt/conda/lib/python3.11/site-packages/geopandas/io/file.py:403: FutureWarning: errors='ignore' is deprecated and will raise in a future version. Use to_datetime without passing `errors` and catch exceptions explicitly instead as_dt = pd.to_datetime(df[k], errors="ignore", utc=True) /opt/conda/lib/python3.11/site-packages/geopandas/io/file.py:399: FutureWarning: errors='ignore' is deprecated and will raise in a future version. Use to_datetime without passing `errors` and catch exceptions explicitly instead as_dt = pd.to_datetime(df[k], errors="ignore") /opt/conda/lib/python3.11/site-packages/geopandas/io/file.py:399: FutureWarning: errors='ignore' is deprecated and will raise in a future version. Use to_datetime without passing `errors` and catch exceptions explicitly instead as_dt = pd.to_datetime(df[k], errors="ignore") /opt/conda/lib/python3.11/site-packages/geopandas/io/file.py:399: FutureWarning: errors='ignore' is deprecated and will raise in a future version. Use to_datetime without passing `errors` and catch exceptions explicitly instead as_dt = pd.to_datetime(df[k], errors="ignore") /opt/conda/lib/python3.11/site-packages/geopandas/io/file.py:399: FutureWarning: errors='ignore' is deprecated and will raise in a future version. Use to_datetime without passing `errors` and catch exceptions explicitly instead as_dt = pd.to_datetime(df[k], errors="ignore") /opt/conda/lib/python3.11/site-packages/geopandas/io/file.py:399: FutureWarning: errors='ignore' is deprecated and will raise in a future version. Use to_datetime without passing `errors` and catch exceptions explicitly instead as_dt = pd.to_datetime(df[k], errors="ignore") /opt/conda/lib/python3.11/site-packages/geopandas/io/file.py:399: FutureWarning: errors='ignore' is deprecated and will raise in a future version. Use to_datetime without passing `errors` and catch exceptions explicitly instead as_dt = pd.to_datetime(df[k], errors="ignore") /opt/conda/lib/python3.11/site-packages/geopandas/io/file.py:399: FutureWarning: errors='ignore' is deprecated and will raise in a future version. Use to_datetime without passing `errors` and catch exceptions explicitly instead as_dt = pd.to_datetime(df[k], errors="ignore") /opt/conda/lib/python3.11/site-packages/geopandas/io/file.py:399: FutureWarning: errors='ignore' is deprecated and will raise in a future version. Use to_datetime without passing `errors` and catch exceptions explicitly instead as_dt = pd.to_datetime(df[k], errors="ignore") /opt/conda/lib/python3.11/site-packages/geopandas/io/file.py:399: FutureWarning: errors='ignore' is deprecated and will raise in a future version. Use to_datetime without passing `errors` and catch exceptions explicitly instead as_dt = pd.to_datetime(df[k], errors="ignore")
Out[7]:
| OBJECTID | poly_SourceOID | poly_IncidentName | poly_FeatureCategory | poly_MapMethod | poly_GISAcres | poly_CreateDate | poly_DateCurrent | poly_PolygonDateTime | poly_IRWINID | ... | attr_Source | attr_IsCpxChild | attr_CpxName | attr_CpxID | attr_SourceGlobalID | GlobalID | Shape__Area | Shape__Length | attr_IncidentComplexityLevel | geometry | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 11275 | 13459 | Grizzly Creek | Wildfire Final Fire Perimeter | GPS-Walked | 322.00 | NaT | 2023-03-14 15:13:42.810000+00:00 | 2020-09-05 03:13:00+00:00 | {3478A04D-5A7E-4E89-9342-80DE610BA9B3} | ... | FODR | None | None | None | {6668D67B-68EE-4CF8-AC3E-EAAD67743DA1} | d5a38ca3-2c28-45eb-b06c-ca4394ecd056 | 0.000143 | 0.079356 | None | POLYGON ((-122.40703 42.27382, -122.40700 42.2... |
| 1 | 14453 | 13272 | Grizzly Creek | Wildfire Final Fire Perimeter | Infrared Image | 32431.62 | NaT | 2023-03-14 15:13:42.810000+00:00 | 2020-09-04 00:37:00+00:00 | {BC150C8C-D9C8-4C14-8725-2B84D7695302} | ... | FODR | None | None | None | {E4E21DCE-C8EE-4133-B9F7-5CA2E5FBF79C} | c9133643-c537-4030-a059-04d80c6c3a96 | 0.013763 | 1.346520 | None | MULTIPOLYGON (((-107.19151 39.66880, -107.1914... |
2 rows × 115 columns
In [8]:
ans_gdf = _
gdf_pts = 0
if isinstance(ans_gdf, gpd.GeoDataFrame):
print('\u2705 Great work! You downloaded and opened a GeoDataFrame')
gdf_pts +=2
else:
print('\u274C Hmm, your answer is not a GeoDataFrame')
print('\u27A1 You earned {} of 2 points for downloading data'.format(gdf_pts))
✅ Great work! You downloaded and opened a GeoDataFrame ➡ You earned 2 of 2 points for downloading data
In [9]:
# Selecting the infared image instance as our geodataframe
fire_url_gdf = fire_url_gdf.loc[[1]]
fire_url_gdf
Out[9]:
| OBJECTID | poly_SourceOID | poly_IncidentName | poly_FeatureCategory | poly_MapMethod | poly_GISAcres | poly_CreateDate | poly_DateCurrent | poly_PolygonDateTime | poly_IRWINID | ... | attr_Source | attr_IsCpxChild | attr_CpxName | attr_CpxID | attr_SourceGlobalID | GlobalID | Shape__Area | Shape__Length | attr_IncidentComplexityLevel | geometry | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 1 | 14453 | 13272 | Grizzly Creek | Wildfire Final Fire Perimeter | Infrared Image | 32431.62 | NaT | 2023-03-14 15:13:42.810000+00:00 | 2020-09-04 00:37:00+00:00 | {BC150C8C-D9C8-4C14-8725-2B84D7695302} | ... | FODR | None | None | None | {E4E21DCE-C8EE-4133-B9F7-5CA2E5FBF79C} | c9133643-c537-4030-a059-04d80c6c3a96 | 0.013763 | 1.34652 | None | MULTIPOLYGON (((-107.19151 39.66880, -107.1914... |
1 rows × 115 columns
In [10]:
# Plot the Grizzly Creek Fire boundary
fire_url_gdf.hvplot(
tiles = 'EsriImagery',
geo = True,
title = 'Grizzly Creek Fire'
)
Out[10]:
Exploring the AppEEARS API for NASA Earthdata access¶
Over the next four cells, you will download MODIS NDVI data for the study period. MODIS is a multispectral instrument that measures Red and NIR data (and so can be used for NDVI). There are two MODIS sensors on two different platforms: satellites Terra and Aqua.
Read more
Learn more about MODIS datasets and the science they support
Since we’re asking for a special download that only covers our study
area, we can’t just find a link to the data - we have to negotiate with
the data server. We’re doing this using the
APPEEARS API
(Application Programming Interface). The API makes it possible for you
to request data using code. You can use code from the earthpy library
to handle the API request.
Your task
Often when we want to do something more complex in coding we find an example and modify it. This download code is already almost a working example. Your task will be:
- Enter your NASA Earthdata username and password when prompted
- Replace the start and end dates in the task parameters. Download data from July, when greenery is at its peak in the Northern Hemisphere.
- Replace the year range. You should get 3 years before and after the fire so you can see the change!
- Replace
gdfwith the name of your site geodataframe.What would the product and layer name be if you were trying to download Landsat Surface Temperature Analysis Ready Data (ARD) instead of MODIS NDVI?
In [11]:
# Initialize AppeearsDownloader for MODIS NDVI data
ndvi_downloader = eaapp.AppeearsDownloader(
download_key='gc-ndvi',
ea_dir=data_dir,
product='MOD13Q1.061', #satellite it is accessing.version
layer='_250m_16_days_NDVI',
start_date="08-01",
end_date="09-15",
recurring=True,
year_range=[2017, 2023],
polygon=fire_url_gdf
)
In [12]:
# Initialize AppeearsDownloader for MODIS NDVI data
ndvi_downloader.download_files(cache=True)
In [13]:
# Get a list of NDVI tif file paths
ndvi_paths = sorted(glob(os.path.join(data_dir, 'gc-ndvi', '*', '*NDVI*.tif')))
ndvi_paths
#len(ndvi_paths)
Out[13]:
['/home/jovyan/grizzly-creek-fire/gc-ndvi/MOD13Q1.061_2017198_to_2023258/MOD13Q1.061__250m_16_days_NDVI_doy2017209_aid0001.tif', '/home/jovyan/grizzly-creek-fire/gc-ndvi/MOD13Q1.061_2017198_to_2023258/MOD13Q1.061__250m_16_days_NDVI_doy2017225_aid0001.tif', '/home/jovyan/grizzly-creek-fire/gc-ndvi/MOD13Q1.061_2017198_to_2023258/MOD13Q1.061__250m_16_days_NDVI_doy2017241_aid0001.tif', '/home/jovyan/grizzly-creek-fire/gc-ndvi/MOD13Q1.061_2017198_to_2023258/MOD13Q1.061__250m_16_days_NDVI_doy2017257_aid0001.tif', '/home/jovyan/grizzly-creek-fire/gc-ndvi/MOD13Q1.061_2017198_to_2023258/MOD13Q1.061__250m_16_days_NDVI_doy2018209_aid0001.tif', '/home/jovyan/grizzly-creek-fire/gc-ndvi/MOD13Q1.061_2017198_to_2023258/MOD13Q1.061__250m_16_days_NDVI_doy2018225_aid0001.tif', '/home/jovyan/grizzly-creek-fire/gc-ndvi/MOD13Q1.061_2017198_to_2023258/MOD13Q1.061__250m_16_days_NDVI_doy2018241_aid0001.tif', '/home/jovyan/grizzly-creek-fire/gc-ndvi/MOD13Q1.061_2017198_to_2023258/MOD13Q1.061__250m_16_days_NDVI_doy2018257_aid0001.tif', '/home/jovyan/grizzly-creek-fire/gc-ndvi/MOD13Q1.061_2017198_to_2023258/MOD13Q1.061__250m_16_days_NDVI_doy2019209_aid0001.tif', '/home/jovyan/grizzly-creek-fire/gc-ndvi/MOD13Q1.061_2017198_to_2023258/MOD13Q1.061__250m_16_days_NDVI_doy2019225_aid0001.tif', '/home/jovyan/grizzly-creek-fire/gc-ndvi/MOD13Q1.061_2017198_to_2023258/MOD13Q1.061__250m_16_days_NDVI_doy2019241_aid0001.tif', '/home/jovyan/grizzly-creek-fire/gc-ndvi/MOD13Q1.061_2017198_to_2023258/MOD13Q1.061__250m_16_days_NDVI_doy2019257_aid0001.tif', '/home/jovyan/grizzly-creek-fire/gc-ndvi/MOD13Q1.061_2017198_to_2023258/MOD13Q1.061__250m_16_days_NDVI_doy2020209_aid0001.tif', '/home/jovyan/grizzly-creek-fire/gc-ndvi/MOD13Q1.061_2017198_to_2023258/MOD13Q1.061__250m_16_days_NDVI_doy2020225_aid0001.tif', '/home/jovyan/grizzly-creek-fire/gc-ndvi/MOD13Q1.061_2017198_to_2023258/MOD13Q1.061__250m_16_days_NDVI_doy2020241_aid0001.tif', '/home/jovyan/grizzly-creek-fire/gc-ndvi/MOD13Q1.061_2017198_to_2023258/MOD13Q1.061__250m_16_days_NDVI_doy2020257_aid0001.tif', '/home/jovyan/grizzly-creek-fire/gc-ndvi/MOD13Q1.061_2017198_to_2023258/MOD13Q1.061__250m_16_days_NDVI_doy2021209_aid0001.tif', '/home/jovyan/grizzly-creek-fire/gc-ndvi/MOD13Q1.061_2017198_to_2023258/MOD13Q1.061__250m_16_days_NDVI_doy2021225_aid0001.tif', '/home/jovyan/grizzly-creek-fire/gc-ndvi/MOD13Q1.061_2017198_to_2023258/MOD13Q1.061__250m_16_days_NDVI_doy2021241_aid0001.tif', '/home/jovyan/grizzly-creek-fire/gc-ndvi/MOD13Q1.061_2017198_to_2023258/MOD13Q1.061__250m_16_days_NDVI_doy2021257_aid0001.tif', '/home/jovyan/grizzly-creek-fire/gc-ndvi/MOD13Q1.061_2017198_to_2023258/MOD13Q1.061__250m_16_days_NDVI_doy2022209_aid0001.tif', '/home/jovyan/grizzly-creek-fire/gc-ndvi/MOD13Q1.061_2017198_to_2023258/MOD13Q1.061__250m_16_days_NDVI_doy2022225_aid0001.tif', '/home/jovyan/grizzly-creek-fire/gc-ndvi/MOD13Q1.061_2017198_to_2023258/MOD13Q1.061__250m_16_days_NDVI_doy2022241_aid0001.tif', '/home/jovyan/grizzly-creek-fire/gc-ndvi/MOD13Q1.061_2017198_to_2023258/MOD13Q1.061__250m_16_days_NDVI_doy2022257_aid0001.tif', '/home/jovyan/grizzly-creek-fire/gc-ndvi/MOD13Q1.061_2017198_to_2023258/MOD13Q1.061__250m_16_days_NDVI_doy2023209_aid0001.tif', '/home/jovyan/grizzly-creek-fire/gc-ndvi/MOD13Q1.061_2017198_to_2023258/MOD13Q1.061__250m_16_days_NDVI_doy2023225_aid0001.tif', '/home/jovyan/grizzly-creek-fire/gc-ndvi/MOD13Q1.061_2017198_to_2023258/MOD13Q1.061__250m_16_days_NDVI_doy2023241_aid0001.tif', '/home/jovyan/grizzly-creek-fire/gc-ndvi/MOD13Q1.061_2017198_to_2023258/MOD13Q1.061__250m_16_days_NDVI_doy2023257_aid0001.tif', '/home/jovyan/grizzly-creek-fire/gc-ndvi/MOD13Q1.061_2017198_to_2023304/MOD13Q1.061__250m_16_days_NDVI_doy2017209_aid0001.tif', '/home/jovyan/grizzly-creek-fire/gc-ndvi/MOD13Q1.061_2017198_to_2023304/MOD13Q1.061__250m_16_days_NDVI_doy2017225_aid0001.tif', '/home/jovyan/grizzly-creek-fire/gc-ndvi/MOD13Q1.061_2017198_to_2023304/MOD13Q1.061__250m_16_days_NDVI_doy2017241_aid0001.tif', '/home/jovyan/grizzly-creek-fire/gc-ndvi/MOD13Q1.061_2017198_to_2023304/MOD13Q1.061__250m_16_days_NDVI_doy2017257_aid0001.tif', '/home/jovyan/grizzly-creek-fire/gc-ndvi/MOD13Q1.061_2017198_to_2023304/MOD13Q1.061__250m_16_days_NDVI_doy2017273_aid0001.tif', '/home/jovyan/grizzly-creek-fire/gc-ndvi/MOD13Q1.061_2017198_to_2023304/MOD13Q1.061__250m_16_days_NDVI_doy2017289_aid0001.tif', '/home/jovyan/grizzly-creek-fire/gc-ndvi/MOD13Q1.061_2017198_to_2023304/MOD13Q1.061__250m_16_days_NDVI_doy2018209_aid0001.tif', '/home/jovyan/grizzly-creek-fire/gc-ndvi/MOD13Q1.061_2017198_to_2023304/MOD13Q1.061__250m_16_days_NDVI_doy2018225_aid0001.tif', '/home/jovyan/grizzly-creek-fire/gc-ndvi/MOD13Q1.061_2017198_to_2023304/MOD13Q1.061__250m_16_days_NDVI_doy2018241_aid0001.tif', '/home/jovyan/grizzly-creek-fire/gc-ndvi/MOD13Q1.061_2017198_to_2023304/MOD13Q1.061__250m_16_days_NDVI_doy2018257_aid0001.tif', '/home/jovyan/grizzly-creek-fire/gc-ndvi/MOD13Q1.061_2017198_to_2023304/MOD13Q1.061__250m_16_days_NDVI_doy2018273_aid0001.tif', '/home/jovyan/grizzly-creek-fire/gc-ndvi/MOD13Q1.061_2017198_to_2023304/MOD13Q1.061__250m_16_days_NDVI_doy2018289_aid0001.tif', '/home/jovyan/grizzly-creek-fire/gc-ndvi/MOD13Q1.061_2017198_to_2023304/MOD13Q1.061__250m_16_days_NDVI_doy2019209_aid0001.tif', '/home/jovyan/grizzly-creek-fire/gc-ndvi/MOD13Q1.061_2017198_to_2023304/MOD13Q1.061__250m_16_days_NDVI_doy2019225_aid0001.tif', '/home/jovyan/grizzly-creek-fire/gc-ndvi/MOD13Q1.061_2017198_to_2023304/MOD13Q1.061__250m_16_days_NDVI_doy2019241_aid0001.tif', '/home/jovyan/grizzly-creek-fire/gc-ndvi/MOD13Q1.061_2017198_to_2023304/MOD13Q1.061__250m_16_days_NDVI_doy2019257_aid0001.tif', '/home/jovyan/grizzly-creek-fire/gc-ndvi/MOD13Q1.061_2017198_to_2023304/MOD13Q1.061__250m_16_days_NDVI_doy2019273_aid0001.tif', '/home/jovyan/grizzly-creek-fire/gc-ndvi/MOD13Q1.061_2017198_to_2023304/MOD13Q1.061__250m_16_days_NDVI_doy2019289_aid0001.tif', '/home/jovyan/grizzly-creek-fire/gc-ndvi/MOD13Q1.061_2017198_to_2023304/MOD13Q1.061__250m_16_days_NDVI_doy2020209_aid0001.tif', '/home/jovyan/grizzly-creek-fire/gc-ndvi/MOD13Q1.061_2017198_to_2023304/MOD13Q1.061__250m_16_days_NDVI_doy2020225_aid0001.tif', '/home/jovyan/grizzly-creek-fire/gc-ndvi/MOD13Q1.061_2017198_to_2023304/MOD13Q1.061__250m_16_days_NDVI_doy2020241_aid0001.tif', '/home/jovyan/grizzly-creek-fire/gc-ndvi/MOD13Q1.061_2017198_to_2023304/MOD13Q1.061__250m_16_days_NDVI_doy2020257_aid0001.tif', '/home/jovyan/grizzly-creek-fire/gc-ndvi/MOD13Q1.061_2017198_to_2023304/MOD13Q1.061__250m_16_days_NDVI_doy2020273_aid0001.tif', '/home/jovyan/grizzly-creek-fire/gc-ndvi/MOD13Q1.061_2017198_to_2023304/MOD13Q1.061__250m_16_days_NDVI_doy2020289_aid0001.tif', '/home/jovyan/grizzly-creek-fire/gc-ndvi/MOD13Q1.061_2017198_to_2023304/MOD13Q1.061__250m_16_days_NDVI_doy2020305_aid0001.tif', '/home/jovyan/grizzly-creek-fire/gc-ndvi/MOD13Q1.061_2017198_to_2023304/MOD13Q1.061__250m_16_days_NDVI_doy2021209_aid0001.tif', '/home/jovyan/grizzly-creek-fire/gc-ndvi/MOD13Q1.061_2017198_to_2023304/MOD13Q1.061__250m_16_days_NDVI_doy2021225_aid0001.tif', '/home/jovyan/grizzly-creek-fire/gc-ndvi/MOD13Q1.061_2017198_to_2023304/MOD13Q1.061__250m_16_days_NDVI_doy2021241_aid0001.tif', '/home/jovyan/grizzly-creek-fire/gc-ndvi/MOD13Q1.061_2017198_to_2023304/MOD13Q1.061__250m_16_days_NDVI_doy2021257_aid0001.tif', '/home/jovyan/grizzly-creek-fire/gc-ndvi/MOD13Q1.061_2017198_to_2023304/MOD13Q1.061__250m_16_days_NDVI_doy2021273_aid0001.tif', '/home/jovyan/grizzly-creek-fire/gc-ndvi/MOD13Q1.061_2017198_to_2023304/MOD13Q1.061__250m_16_days_NDVI_doy2021289_aid0001.tif', '/home/jovyan/grizzly-creek-fire/gc-ndvi/MOD13Q1.061_2017198_to_2023304/MOD13Q1.061__250m_16_days_NDVI_doy2022209_aid0001.tif', '/home/jovyan/grizzly-creek-fire/gc-ndvi/MOD13Q1.061_2017198_to_2023304/MOD13Q1.061__250m_16_days_NDVI_doy2022225_aid0001.tif', '/home/jovyan/grizzly-creek-fire/gc-ndvi/MOD13Q1.061_2017198_to_2023304/MOD13Q1.061__250m_16_days_NDVI_doy2022241_aid0001.tif', '/home/jovyan/grizzly-creek-fire/gc-ndvi/MOD13Q1.061_2017198_to_2023304/MOD13Q1.061__250m_16_days_NDVI_doy2022257_aid0001.tif', '/home/jovyan/grizzly-creek-fire/gc-ndvi/MOD13Q1.061_2017198_to_2023304/MOD13Q1.061__250m_16_days_NDVI_doy2022273_aid0001.tif', '/home/jovyan/grizzly-creek-fire/gc-ndvi/MOD13Q1.061_2017198_to_2023304/MOD13Q1.061__250m_16_days_NDVI_doy2022289_aid0001.tif', '/home/jovyan/grizzly-creek-fire/gc-ndvi/MOD13Q1.061_2017198_to_2023304/MOD13Q1.061__250m_16_days_NDVI_doy2023209_aid0001.tif', '/home/jovyan/grizzly-creek-fire/gc-ndvi/MOD13Q1.061_2017198_to_2023304/MOD13Q1.061__250m_16_days_NDVI_doy2023225_aid0001.tif', '/home/jovyan/grizzly-creek-fire/gc-ndvi/MOD13Q1.061_2017198_to_2023304/MOD13Q1.061__250m_16_days_NDVI_doy2023241_aid0001.tif', '/home/jovyan/grizzly-creek-fire/gc-ndvi/MOD13Q1.061_2017198_to_2023304/MOD13Q1.061__250m_16_days_NDVI_doy2023257_aid0001.tif', '/home/jovyan/grizzly-creek-fire/gc-ndvi/MOD13Q1.061_2017198_to_2023304/MOD13Q1.061__250m_16_days_NDVI_doy2023273_aid0001.tif', '/home/jovyan/grizzly-creek-fire/gc-ndvi/MOD13Q1.061_2017198_to_2023304/MOD13Q1.061__250m_16_days_NDVI_doy2023289_aid0001.tif']
Repeating tasks in Python¶
Now you should have a few dozen files! For each file, you need to:
- Load the file in using the
rioxarraylibrary - Get the date from the file name
- Add the date as a dimension coordinate
- Give your data variable a name
- Divide by the scale factor of 10000
You don’t want to write out the code for each file! That’s a recipe for
copy pasta. Luckily, Python has tools for doing similar tasks
repeatedly. In this case, you’ll use one called a for loop.
There’s some code below that uses a for loop in what is called an
accumulation pattern to process each file. That means that you will
save the results of your processing to a list each time you process the
files, and then merge all the arrays in the list.
Your task
- Look at the file names. How many characters from the end is the date?
- Replace any required variable names with your chosen variable names
- Change the
scale_factorvariable to be the correct scale factor for this NDVI dataset (HINT: NDVI should range between 0 and 1)- Using indices or regular
In [15]:
scale_factor = 10000
doy_start = -19
doy_end = -12
In [16]:
ndvi_das = []
for ndvi_path in ndvi_paths:
# Get date from file name
doy = ndvi_path[doy_start:doy_end]
date = pd.to_datetime(doy, format='%Y%j')
# Open dataset
da = rxr.open_rasterio(ndvi_path, masked=True).squeeze()
# Add date dimension and clean up metadata
da = da.assign_coords({'date': date})
da = da.expand_dims({'date': 1})
da.name = 'NDVI'
# Multiple by scale factor
da = da / scale_factor
# Prepare for concatenation
ndvi_das.append(da)
len(ndvi_das)
Out[16]:
71
Next, stack your arrays by date into a time series using the
xr.combine_by_coords() function. You will have to tell it which
dimension you want to stack your data in.
Plot the change in NDVI spatially
Complete the following: * Select data from 2021 to 2023 (3 years after the fire) * Take the temporal mean (over the date, not spatially) * Get the NDVI variable (should be a DataArray, not a Dataset) * Repeat for the data from 2018 to 2020 (3 years before the fire) * Subtract the 2018-2020 time period from the 2021-2023 time period * Plot the result using a diverging color map like
cmap=plt.cm.PiYGThere are different types of color maps for different types of data. In this case, we want decreases to be a different color from increases, so we should use a diverging color map. Check out available colormaps in the matplotlib documentation.
For an extra challenge, add the fire boundary to the plot
In [18]:
ndvi_da = xr.combine_by_coords(ndvi_das, coords=['date'])
ndvi_da
Out[18]:
<xarray.Dataset>
Dimensions: (x: 94, y: 63, date: 43)
Coordinates:
band int64 1
* x (x) float64 -107.3 -107.3 -107.3 ... -107.1 -107.1 -107.1
* y (y) float64 39.67 39.67 39.67 39.66 ... 39.55 39.54 39.54 39.54
spatial_ref int64 0
* date (date) datetime64[ns] 2017-07-28 2017-08-13 ... 2023-10-16
Data variables:
NDVI (date, y, x) float32 0.6496 0.7459 0.7599 ... 0.4019 0.3721In [25]:
# Compute the difference in NDVI before and after the fire
ndvi_diff = (
ndvi_da
.sel(date=slice('2021', '2023'))
.mean('date')
.NDVI
- ndvi_da
.sel(date=slice('2018', '2020'))
.mean('date')
.NDVI
)
# Plot the difference
(
ndvi_diff.hvplot(x='x', y='y', cmap='PiYG', geo=True)
*
fire_url_gdf.hvplot(geo=True, fill_color=None, line_color='black')
)
Out[25]:
Is the NDVI lower within the fire boundary after the fire?¶
You will compute the mean NDVI inside and outside the fire boundary.
First, use the code below to get a GeoDataFrame of the area outside
the Reservation. Your task: * Check the variable names - Make sure that
the code uses your boundary GeoDataFrame * How could you test if the
geometry was modified correctly? Add some code to take a look at the
results.
In [26]:
out_gdf = (
gpd.GeoDataFrame(geometry=fire_url_gdf.envelope)
.overlay(fire_url_gdf, how='difference'))
# testing if geometry modified correctly
out_gdf.hvplot(geo=True, fill_color = 'blue', line_color='black')
Out[26]:
Next, clip your DataArray to the boundaries for both inside and outside
the reservation. You will need to replace the GeoDataFrame name with
your own. Check out the lesson on clipping data with the rioxarray
library in the
textbook.
GOTCHA ALERT
It’s important to use
from_disk=Truewhen clipping large arrays like this. It allows the computer to use less valuable memory resources when clipping - you will probably find that otherwise the cell below crashes the kernel
In [27]:
# Clip data to both inside and outside the boundary
ndvi_gc_da = ndvi_da.rio.clip(fire_url_gdf.geometry, from_disk=True)
ndvi_out_da = ndvi_da.rio.clip(out_gdf.geometry, from_disk=True)
Your Task
For both inside and outside the fire boundary:
- Group the data by year
- Take the mean. You always need to tell reducing methods in
xarraywhat dimensions you want to reduce. When you want to summarize data across all dimensions, you can use the...syntax, e.g..mean(...)as a shorthand.- Select the NDVI variable
- Convert to a DataFrame using the
to_dataframe()method- Join the two DataFrames for plotting using the
.join()method. You will need to rename the columns using thelsuffix=andrsuffix=parameters
GOTCHA ALERT
The DateIndex in pandas is a little different from the Datetime Dimension in xarray. You will need to use the
.dt.yearsyntax to access information about the year, not just.year.
Finally, plot annual July means for both inside and outside the Reservation on the same plot.
:::
In [28]:
# Compute mean annual July NDVI
jul_ndvi_gc_df = (
ndvi_gc_da
.groupby(ndvi_gc_da.date.dt.year)
.mean(...)
.NDVI.to_dataframe())
jul_ndvi_out_df = (
ndvi_out_da
.groupby(ndvi_out_da.date.dt.year)
.mean(...)
.NDVI.to_dataframe())
# Plot inside and outside the reservation
jul_ndvi_df = (
jul_ndvi_gc_df[['NDVI']]
.join(
jul_ndvi_out_df[['NDVI']],
lsuffix=' Burned Area', rsuffix=' Unburned Area')
)
jul_ndvi_df.hvplot(
title='NDVI before and after the Grizzly Creek Fire'
)
Out[28]:
Now, take the difference between outside and inside the Reservation and plot that. What do you observe? Don’t forget to write a headline and description of your plot!
In [29]:
# Plot difference inside and outside the reservation
jul_ndvi_df['difference'] = (
jul_ndvi_df['NDVI Burned Area']
- jul_ndvi_df['NDVI Unburned Area'])
jul_ndvi_df.difference.hvplot(
title='Difference between NDVI within and outside the Grizzly Creek Fire'
)
Out[29]:
Your turn! Repeat this workflow in a different time and place.¶
It’s not just fires that affect NDVI! You could look at:
- Recovery after a national disaster, like a wildfire or hurricane
- Drought
- Deforestation
- Irrigation
- Beaver reintroduction
In [24]:
%%capture
%%bash
jupyter nbconvert *.ipynb --to html