Assignment : Vegetation Coding Challenge 🌏 📈¶
Study Area: The Cherry Canyon Fire, Colorado, USA
Cherry Canyon Fire broke out on 24/5/2020, 13 miles north of Kim, Colorado. As stated, it was caused due to Strong winds and extremely low humidity.
In this assignment, three years pre 2018-2020 and post fire 2020-2023 is analyzed using MODIS NDVI data.
Data: MODIS NDVI Data
Duration: 2018-2023 (3 years pre 2018-2020 and post fires 2020-2023)
Fire_Event: 2020
API Used: AppEEARS API for NASA Earthdata access
url: https://appeears.earthdatacloud.nasa.gov/api/
In [1]:
# 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-hrgp4vnk Running command git clone --filter=blob:none --quiet https://github.com/earthlab/earthpy /tmp/pip-req-build-hrgp4vnk 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 20e0b17e563da17ba46dc32832cf4ed9173a7531 Installing build dependencies ... done Getting requirements to build wheel ... done Preparing metadata (pyproject.toml) ... done Requirement already satisfied: geopandas in /opt/conda/lib/python3.10/site-packages (from earthpy==0.10.0) (0.13.2) Requirement already satisfied: matplotlib>=2.0.0 in /opt/conda/lib/python3.10/site-packages (from earthpy==0.10.0) (3.7.2) Requirement already satisfied: numpy>=1.14.0 in /opt/conda/lib/python3.10/site-packages (from earthpy==0.10.0) (1.25.2) Requirement already satisfied: rasterio in /opt/conda/lib/python3.10/site-packages (from earthpy==0.10.0) (1.3.8) Requirement already satisfied: scikit-image in /opt/conda/lib/python3.10/site-packages (from earthpy==0.10.0) (0.21.0) Requirement already satisfied: requests in /opt/conda/lib/python3.10/site-packages (from earthpy==0.10.0) (2.31.0) Requirement already satisfied: keyring in /opt/conda/lib/python3.10/site-packages (from earthpy==0.10.0) (25.2.0) Requirement already satisfied: contourpy>=1.0.1 in /opt/conda/lib/python3.10/site-packages (from matplotlib>=2.0.0->earthpy==0.10.0) (1.1.0) Requirement already satisfied: cycler>=0.10 in /opt/conda/lib/python3.10/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.10/site-packages (from matplotlib>=2.0.0->earthpy==0.10.0) (4.42.1) Requirement already satisfied: kiwisolver>=1.0.1 in /opt/conda/lib/python3.10/site-packages (from matplotlib>=2.0.0->earthpy==0.10.0) (1.4.5) Requirement already satisfied: packaging>=20.0 in /opt/conda/lib/python3.10/site-packages (from matplotlib>=2.0.0->earthpy==0.10.0) (23.1) Requirement already satisfied: pillow>=6.2.0 in /opt/conda/lib/python3.10/site-packages (from matplotlib>=2.0.0->earthpy==0.10.0) (10.0.0) Requirement already satisfied: pyparsing<3.1,>=2.3.1 in /opt/conda/lib/python3.10/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.10/site-packages (from matplotlib>=2.0.0->earthpy==0.10.0) (2.8.2) Requirement already satisfied: fiona>=1.8.19 in /opt/conda/lib/python3.10/site-packages (from geopandas->earthpy==0.10.0) (1.9.4) Requirement already satisfied: pandas>=1.1.0 in /opt/conda/lib/python3.10/site-packages (from geopandas->earthpy==0.10.0) (2.1.0) Requirement already satisfied: pyproj>=3.0.1 in /opt/conda/lib/python3.10/site-packages (from geopandas->earthpy==0.10.0) (3.6.0) Requirement already satisfied: shapely>=1.7.1 in /opt/conda/lib/python3.10/site-packages (from geopandas->earthpy==0.10.0) (2.0.1) Requirement already satisfied: jaraco.classes in /opt/conda/lib/python3.10/site-packages (from keyring->earthpy==0.10.0) (3.4.0) Requirement already satisfied: jaraco.functools in /opt/conda/lib/python3.10/site-packages (from keyring->earthpy==0.10.0) (4.0.1) Requirement already satisfied: jaraco.context in /opt/conda/lib/python3.10/site-packages (from keyring->earthpy==0.10.0) (5.3.0) Requirement already satisfied: importlib-metadata>=4.11.4 in /opt/conda/lib/python3.10/site-packages (from keyring->earthpy==0.10.0) (6.8.0) Requirement already satisfied: SecretStorage>=3.2 in /opt/conda/lib/python3.10/site-packages (from keyring->earthpy==0.10.0) (3.3.3) Requirement already satisfied: jeepney>=0.4.2 in /opt/conda/lib/python3.10/site-packages (from keyring->earthpy==0.10.0) (0.8.0) Requirement already satisfied: affine in /opt/conda/lib/python3.10/site-packages (from rasterio->earthpy==0.10.0) (2.4.0) Requirement already satisfied: attrs in /opt/conda/lib/python3.10/site-packages (from rasterio->earthpy==0.10.0) (23.1.0) Requirement already satisfied: certifi in /opt/conda/lib/python3.10/site-packages (from rasterio->earthpy==0.10.0) (2023.7.22) Requirement already satisfied: click>=4.0 in /opt/conda/lib/python3.10/site-packages (from rasterio->earthpy==0.10.0) (8.1.7) Requirement already satisfied: cligj>=0.5 in /opt/conda/lib/python3.10/site-packages (from rasterio->earthpy==0.10.0) (0.7.2) Requirement already satisfied: snuggs>=1.4.1 in /opt/conda/lib/python3.10/site-packages (from rasterio->earthpy==0.10.0) (1.4.7) Requirement already satisfied: click-plugins in /opt/conda/lib/python3.10/site-packages (from rasterio->earthpy==0.10.0) (1.1.1) Requirement already satisfied: setuptools in /opt/conda/lib/python3.10/site-packages (from rasterio->earthpy==0.10.0) (68.1.2) Requirement already satisfied: charset-normalizer<4,>=2 in /opt/conda/lib/python3.10/site-packages (from requests->earthpy==0.10.0) (3.2.0) Requirement already satisfied: idna<4,>=2.5 in /opt/conda/lib/python3.10/site-packages (from requests->earthpy==0.10.0) (3.4) Requirement already satisfied: urllib3<3,>=1.21.1 in /opt/conda/lib/python3.10/site-packages (from requests->earthpy==0.10.0) (2.0.4) Requirement already satisfied: scipy>=1.8 in /opt/conda/lib/python3.10/site-packages (from scikit-image->earthpy==0.10.0) (1.11.2) Requirement already satisfied: networkx>=2.8 in /opt/conda/lib/python3.10/site-packages (from scikit-image->earthpy==0.10.0) (3.1) Requirement already satisfied: imageio>=2.27 in /opt/conda/lib/python3.10/site-packages (from scikit-image->earthpy==0.10.0) (2.31.1) Requirement already satisfied: tifffile>=2022.8.12 in /opt/conda/lib/python3.10/site-packages (from scikit-image->earthpy==0.10.0) (2023.8.30) Requirement already satisfied: PyWavelets>=1.1.1 in /opt/conda/lib/python3.10/site-packages (from scikit-image->earthpy==0.10.0) (1.4.1) Requirement already satisfied: lazy_loader>=0.2 in /opt/conda/lib/python3.10/site-packages (from scikit-image->earthpy==0.10.0) (0.3) Requirement already satisfied: six in /opt/conda/lib/python3.10/site-packages (from fiona>=1.8.19->geopandas->earthpy==0.10.0) (1.16.0) Requirement already satisfied: zipp>=0.5 in /opt/conda/lib/python3.10/site-packages (from importlib-metadata>=4.11.4->keyring->earthpy==0.10.0) (3.16.2) Requirement already satisfied: pytz>=2020.1 in /opt/conda/lib/python3.10/site-packages (from pandas>=1.1.0->geopandas->earthpy==0.10.0) (2023.3.post1) Requirement already satisfied: tzdata>=2022.1 in /opt/conda/lib/python3.10/site-packages (from pandas>=1.1.0->geopandas->earthpy==0.10.0) (2023.3) Requirement already satisfied: cryptography>=2.0 in /opt/conda/lib/python3.10/site-packages (from SecretStorage>=3.2->keyring->earthpy==0.10.0) (41.0.3) Requirement already satisfied: more-itertools in /opt/conda/lib/python3.10/site-packages (from jaraco.classes->keyring->earthpy==0.10.0) (10.2.0) Requirement already satisfied: backports.tarfile in /opt/conda/lib/python3.10/site-packages (from jaraco.context->keyring->earthpy==0.10.0) (1.1.1) Requirement already satisfied: cffi>=1.12 in /opt/conda/lib/python3.10/site-packages (from cryptography>=2.0->SecretStorage>=3.2->keyring->earthpy==0.10.0) (1.15.1) Requirement already satisfied: pycparser in /opt/conda/lib/python3.10/site-packages (from cffi>=1.12->cryptography>=2.0->SecretStorage>=3.2->keyring->earthpy==0.10.0) (2.21)
In [2]:
import getpass
import json
import os
import pathlib
from glob import glob
# Library for tabular data
import pandas as pd
#Library for vector data
import geopandas as gpd
import earthpy.appeears as eaapp
import hvplot.pandas
import hvplot.xarray
import rioxarray as rxr
import xarray as xr
In [3]:
data_dir = os.path.join(pathlib.Path.home(), 'cherry-data')
# Make the data directory
os.makedirs(data_dir, exist_ok=True)
#To visualize the directory path
data_dir
Out[3]:
'/home/jovyan/cherry-data'
In [4]:
# Download the Cherry Canyon fire boundary
cherry_url = ("https://services3.arcgis.com/T4QMspbfLg3qTGWY/arcgis"
"/rest/services/WFIGS_Interagency_Perimeters/FeatureServer/0/"
"query?where=poly_IncidentName%20%3D%20'CHERRY%20CANYON'&outFields=*&outSR=4326&f=json")
cherry_url
#Opening data with geopandas
cherry_gdf = gpd.read_file(cherry_url)
cherry_gdf
Out[4]:
| OBJECTID | poly_SourceOID | poly_IncidentName | poly_FeatureCategory | poly_MapMethod | poly_GISAcres | poly_CreateDate | poly_DateCurrent | poly_PolygonDateTime | poly_IRWINID | ... | attr_ModifiedOnDateTime_dt | attr_Source | attr_IsCpxChild | attr_CpxName | attr_CpxID | attr_SourceGlobalID | GlobalID | Shape__Area | Shape__Length | geometry | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 15518 | NaN | Cherry Canyon | Wildfire Daily Fire Perimeter | Hand Sketch | 11756.34 | 1679512187728 | NaN | NaN | {F32D0475-5315-4448-9C2A-9272C5B0F6FE} | ... | 1593543458777 | IRWIN | 0 | NaN | NaN | {F32D0475-5315-4448-9C2A-9272C5B0F6FE} | 5c500faf-4dad-4819-b5bb-ebba2d85a578 | 0.004843 | 0.812748 | MULTIPOLYGON (((-103.43758 37.47623, -103.4375... |
1 rows × 114 columns
In [10]:
cherry_gdf.explore()
Out[10]:
Make this Notebook Trusted to load map: File -> Trust Notebook
In [11]:
# Initialize AppeearsDownloader for MODIS NDVI data
ndvi_downloader = eaapp.AppeearsDownloader(
download_key='cherry-ndvi',
ea_dir=data_dir,
product='MOD13Q1.061',
layer='_250m_16_days_NDVI',
start_date="07-01",
end_date="07-31",
recurring=True,
year_range=[2018, 2023],
polygon=cherry_gdf
)
ndvi_downloader.download_files(cache=True)
In [5]:
# Get a list of NDVI tif file paths
ndvi_files = sorted(glob(os.path.join(data_dir, 'cherry-ndvi', '*', '*NDVI*.tif')))
len(ndvi_files)
Out[5]:
18
In [6]:
scale_factor = 10000
d_start = -19
d_end = -12
In [7]:
ndvi_das = []
for ndvi_path in ndvi_files:
# Get date from file name
doy = ndvi_path[d_start:d_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[7]:
18
In [8]:
ndvi_da = xr.combine_by_coords(ndvi_das, coords=['date'])
ndvi_da
Out[8]:
<xarray.Dataset>
Dimensions: (x: 49, y: 55, date: 18)
Coordinates:
band int64 1
* x (x) float64 -103.5 -103.5 -103.5 ... -103.4 -103.4 -103.4
* y (y) float64 37.48 37.48 37.48 37.47 ... 37.37 37.37 37.37 37.37
spatial_ref int64 0
* date (date) datetime64[ns] 2018-06-26 2018-07-12 ... 2023-07-28
Data variables:
NDVI (date, y, x) float32 0.2481 0.2481 0.2036 ... 0.2931 0.2931In [9]:
# Compute the difference in NDVI before and after the fire
# Plot the difference
import matplotlib.pyplot as plt
fig, ax = plt.subplots()
ndvi_diff = (
ndvi_da
.sel(date=slice('2021', '2023'))
.mean('date')
.NDVI
- ndvi_da
.sel(date=slice('2018', '2020'))
.mean('date')
.NDVI
)
ndvi_diff.plot(x='x',y='y',cmap='PiYG', ax= ax)
cherry_gdf.plot(color = 'None' , edgecolor='black', ax= ax)
Out[9]:
<Axes: title={'center': 'band = 1, spatial_ref = 0'}, xlabel='x', ylabel='y'>
Higher values of NDVI that is towards postive values shows a healthy vegetation whereas negative, near zero values shows stressed vegetation and may represent features like water body, exposed soil, urban area. In the above figure, healthy vegetation are shown in green colors and stressed or other features with pink and white colors. Negative values are shown using dark pink, near zero values with white, and positive values with green.
In [10]:
area_out_gdf = (
gpd.GeoDataFrame(geometry=cherry_gdf.envelope)
.overlay(cherry_gdf, how='difference'))
In [11]:
# Clip data to both inside and outside the boundary
ndvi_cherry_da = ndvi_da.rio.clip(cherry_gdf.geometry, from_disk=True)
ndvi_out_da = ndvi_da.rio.clip(area_out_gdf.geometry, from_disk=True)
In [12]:
# Compute mean annual July NDVI
jul_ndvi_cherry_df = (
ndvi_cherry_da
.groupby(ndvi_cherry_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_cherry_df[['NDVI']]
.join(
jul_ndvi_out_df[['NDVI']],
lsuffix=' Burned Area', rsuffix=' Unburned Area')
)
jul_ndvi_df.hvplot(
title='NDVI before and after the Cherry Canyon Fire'
)
/opt/conda/lib/python3.10/site-packages/holoviews/core/data/pandas.py:39: FutureWarning: Series.__getitem__ treating keys as positions is deprecated. In a future version, integer keys will always be treated as labels (consistent with DataFrame behavior). To access a value by position, use `ser.iloc[pos]` return dataset.data.dtypes[idx].type /opt/conda/lib/python3.10/site-packages/holoviews/core/data/pandas.py:39: FutureWarning: Series.__getitem__ treating keys as positions is deprecated. In a future version, integer keys will always be treated as labels (consistent with DataFrame behavior). To access a value by position, use `ser.iloc[pos]` return dataset.data.dtypes[idx].type
Out[12]:
Analysis of the Graph¶
In 2019, the NDVI values of burned and Unburned areas indicates a healthy vegetation with values greater than 0.35. However, in 2020, lower NDVI values of burned area indicates stressed vegetation due to fire. Post fire event, the vegetation of burned area shows an upward trend indicating improvement in vegetation health with time in 2021 and afterwards.
In [13]:
# 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 Cherry Canyon Fire'
)
/opt/conda/lib/python3.10/site-packages/holoviews/core/data/pandas.py:39: FutureWarning: Series.__getitem__ treating keys as positions is deprecated. In a future version, integer keys will always be treated as labels (consistent with DataFrame behavior). To access a value by position, use `ser.iloc[pos]` return dataset.data.dtypes[idx].type /opt/conda/lib/python3.10/site-packages/holoviews/core/data/pandas.py:39: FutureWarning: Series.__getitem__ treating keys as positions is deprecated. In a future version, integer keys will always be treated as labels (consistent with DataFrame behavior). To access a value by position, use `ser.iloc[pos]` return dataset.data.dtypes[idx].type
Out[13]:
In the year 2020, as it is seen from the graph the value of NDVI is negative (less than -0.06). This suggests a very stressed vegetation. If compared from 2018 to 2020, there is depreciating value of NDVI from greater then 0.02 to -0.06, whereas the value of NDVI rises from -0.06 to -0.03 from 2020 to 2023.
In [14]:
%%capture
%%bash
#Covert to HTML
jupyter nbconvert vegetation.ipynb --to html