The Cameron Peak Fire, Colorado, USA

The Cameron Peak Fire was the largest fire in Colorado history, with 326 square miles burned.

Observing vegetation health from space

We will look at the destruction and recovery of vegetation using NDVI (Normalized Difference Vegetation Index). How does it work? First, we need to learn about spectral reflectance signatures.

Every object reflects some wavelengths of light more or less than others. We can see this with our eyes, since, for example, plants reflect a lot of green in the summer, and then as that green diminishes in the fall they look more yellow or orange. The image below shows spectral signatures for water, soil, and vegetation:

> Image source: SEOS Project

Healthy vegetation reflects a lot of Near-InfraRed (NIR) radiation. Dead ve Less healthy vegetation reflects a similar amounts of the visible light spectra, but less NIR radiation. We don’t see a huge drop in Green radiation until the plant is very stressed or dead. That means that NIR allows us to get ahead of what we can see with our eyes.

> Image source: Spectral signature literature review by px39n

Different species of plants reflect different spectral signatures, but the pattern of the signatures are similar. NDVI compares the amount of NIR reflectance to the amount of Red reflectance, thus accounting for many of the species differences and isolating the health of the plant. The formula for calculating NDVI is:

\[NDVI = \frac{(NIR - Red)}{(NIR + Red)}\]

Read more about NDVI and other vegetation indices: * earthdatascience.org * USGS

Import necessary libraries

In the cell below, making sure to keep the packages in order, add packages for:

• Working with DataFrames
• Working with GeoDataFrames

What are we using the rest of these packages for? See if you can figure it out as you complete the notebook.

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 getpass
import json
import os
import pathlib
from glob import glob

import earthpy.appeears as eaapp
import geopandas as gpd
import hvplot.pandas
import hvplot.xarray
import pandas as pd
import rioxarray as rxr
import xarray as xr

We have one more setup task. We’re not going to be able to load all our data directly from the web to Python this time. That means we need to set up a place for it.

GOTCHA ALERT

A lot of times in Python we say “directory” to mean a “folder” on your computer. The two words mean the same thing in this context.

Your task

In the cell below, replace ‘my-data-folder’ with a descriptive directory name.

data_dir = os.path.join(pathlib.Path.home(), 'my-data-folder')
# Make the data directory
os.makedirs(data_dir, exist_ok=True)

data_dir = os.path.join(
    pathlib.Path.home(), 'earth-analytics', 'data', 'cameron-peak')
# Make the data directory
os.makedirs(data_dir, exist_ok=True)

Study Area: Cameron Peak Fire Boundary

Earth Data Science data formats

In Earth Data Science, we get data in three main formats:

Data type	Descriptions	Common file formats	Python type
Time Series	The same data points (e.g. streamflow) collected multiple times over time	Tabular formats (e.g. .csv, or .xlsx)	pandas DataFrame
Vector	Points, lines, and areas (with coordinates)	Shapefile (often an archive like a .zip file because a Shapefile is actually a collection of at least 3 files)	geopandas GeoDataFrame
Raster	Evenly spaced spatial grid (with coordinates)	GeoTIFF (.tif), NetCDF (.nc), HDF (.hdf)	rioxarray DataArray

Read more

Check out the sections about about vector data and raster data in the textbook.

What do you think?

For this coding challenge, we are interested in the boundary of the Cameron Peak Fire. In the cell below, answer the following question: What data type do you think the reservation boundaries will be?

Your task:

• Search the National Interagency Fire Center Wildfire Boundary catalog for and incident name “Cameron Peak”
• Copy the API results to your clipboard.
• Load the data into Python using the geopandas library, e.g.: python gpd.read_file(url)
• Call your data at the end of the cell for testing.

# Download the Cameron Peak fire boundary

url = (
    "https://services3.arcgis.com/T4QMspbfLg3qTGWY/arcgis/rest/services"
    "/WFIGS_Interagency_Perimeters/FeatureServer/0/query"
    "?where=poly_IncidentName%20%3D%20'CAMERON%20PEAK'"
    "&outFields=*&outSR=4326&f=json")

gdf = gpd.read_file(url)
gdf

	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	14659	10393	Cameron Peak	Wildfire Final Fire Perimeter	Infrared Image	208913.3	NaT	2023-03-14 15:13:42.810000+00:00	2020-10-06 21:30:00+00:00	{53741A13-D269-4CD5-AF91-02E094B944DA}	...	2021-12-06 19:48:12.647000+00:00	FODR	None	None	None	{5F5EC9BA-5F85-495B-8B97-1E0969A1434E}	e8e4ffd4-db84-4d47-bd32-d4b0a8381fff	0.089957	5.541765	MULTIPOLYGON (((-105.88333 40.54598, -105.8837...

1 rows × 114 columns

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

Site Map

We always want to create a site map when working with geospatial data. This helps us see that we’re looking at the right location, and learn something about the context of the analysis.

Your task

• Plot your Cameron Peak Fire shapefile on an interactive map
• Make sure to add a title
• Add ESRI World Imagery as the basemap/background using the tiles=... parameter

# Plot the Cameron Peak Fire boundary

gdf.hvplot(
    title='Cameron Peak Fire, 2020',
    tiles='EsriImagery')

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 gdf with 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?

# Initialize AppeearsDownloader for MODIS NDVI data
ndvi_downloader = eaapp.AppeearsDownloader(
    download_key='cp-ndvi',
    ea_dir=data_dir,
    product='MOD13Q1.061',
    layer='_250m_16_days_NDVI',
    start_date="01-01",
    end_date="01-31",
    recurring=True,
    year_range=[2021, 2021],
    polygon=gdf
)

# Initialize AppeearsDownloader for MODIS NDVI data
ndvi_downloader = eaapp.AppeearsDownloader(
    download_key='cp-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=gdf
)
ndvi_downloader.download_files(cache=True)

** Message: 15:58:57.832: Remote error from secret service: org.freedesktop.DBus.Error.ServiceUnknown: The name org.freedesktop.secrets was not provided by any .service files

Putting it together: Working with multi-file raster datasets in Python

Now you need to load all the downloaded files into Python. Let’s start by getting all the file names. You will also need to extract the date from the filename. Check out the lesson on getting information from filenames in the textbook.

GOTCHA ALERT

glob doesn’t necessarily find files in the order you would expect. Make sure to sort your file names like it says in the textbook.

# Get a list of NDVI tif file paths

# Get list of NDVI tif file paths
ndvi_paths = sorted(glob(os.path.join(data_dir, 'cp-ndvi', '*', '*NDVI*.tif')))
len(ndvi_paths)

18

Repeating tasks in Python

Now you should have a few dozen files! For each file, you need to:

• Load the file in using the rioxarray library
• 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_factor variable to be the correct scale factor for this NDVI dataset (HINT: NDVI should range between 0 and 1)
• Using indices or regular

scale_factor = 1
doy_start = -1
doy_end = -1

scale_factor = 10000
doy_start = -19
doy_end = -12

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)

18

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.

ndvi_da = xr.combine_by_coords(ndvi_das, coords=['date'])
ndvi_da

<xarray.Dataset> Size: 4MB
Dimensions:      (date: 18, y: 153, x: 339)
Coordinates:
    band         int64 8B 1
  * x            (x) float64 3kB -105.9 -105.9 -105.9 ... -105.2 -105.2 -105.2
  * y            (y) float64 1kB 40.78 40.78 40.77 40.77 ... 40.47 40.46 40.46
    spatial_ref  int64 8B 0
  * date         (date) datetime64[ns] 144B 2018-06-26 2018-07-12 ... 2023-07-28
Data variables:
    NDVI         (date, y, x) float32 4MB 0.5672 0.5711 0.6027 ... 0.6327 0.6327

xarray.Dataset

• Dimensions:
  • date: 18
  • y: 153
  • x: 339
• Coordinates: (5)
  • band
    ()
    int64
    1

    array(1)

  • x
    (x)
    float64
    -105.9 -105.9 ... -105.2 -105.2

    array([-105.903125, -105.901042, -105.898958, ..., -105.203125, -105.201042,
           -105.198958])

  • y
    (y)
    float64
    40.78 40.78 40.77 ... 40.46 40.46

    array([40.778125, 40.776042, 40.773958, 40.771875, 40.769792, 40.767708,
           40.765625, 40.763542, 40.761458, 40.759375, 40.757292, 40.755208,
           40.753125, 40.751042, 40.748958, 40.746875, 40.744792, 40.742708,
           40.740625, 40.738542, 40.736458, 40.734375, 40.732292, 40.730208,
           40.728125, 40.726042, 40.723958, 40.721875, 40.719792, 40.717708,
           40.715625, 40.713542, 40.711458, 40.709375, 40.707292, 40.705208,
           40.703125, 40.701042, 40.698958, 40.696875, 40.694792, 40.692708,
           40.690625, 40.688542, 40.686458, 40.684375, 40.682292, 40.680208,
           40.678125, 40.676042, 40.673958, 40.671875, 40.669792, 40.667708,
           40.665625, 40.663542, 40.661458, 40.659375, 40.657292, 40.655208,
           40.653125, 40.651042, 40.648958, 40.646875, 40.644792, 40.642708,
           40.640625, 40.638542, 40.636458, 40.634375, 40.632292, 40.630208,
           40.628125, 40.626042, 40.623958, 40.621875, 40.619792, 40.617708,
           40.615625, 40.613542, 40.611458, 40.609375, 40.607292, 40.605208,
           40.603125, 40.601042, 40.598958, 40.596875, 40.594792, 40.592708,
           40.590625, 40.588542, 40.586458, 40.584375, 40.582292, 40.580208,
           40.578125, 40.576042, 40.573958, 40.571875, 40.569792, 40.567708,
           40.565625, 40.563542, 40.561458, 40.559375, 40.557292, 40.555208,
           40.553125, 40.551042, 40.548958, 40.546875, 40.544792, 40.542708,
           40.540625, 40.538542, 40.536458, 40.534375, 40.532292, 40.530208,
           40.528125, 40.526042, 40.523958, 40.521875, 40.519792, 40.517708,
           40.515625, 40.513542, 40.511458, 40.509375, 40.507292, 40.505208,
           40.503125, 40.501042, 40.498958, 40.496875, 40.494792, 40.492708,
           40.490625, 40.488542, 40.486458, 40.484375, 40.482292, 40.480208,
           40.478125, 40.476042, 40.473958, 40.471875, 40.469792, 40.467708,
           40.465625, 40.463542, 40.461458])

  • spatial_ref
    ()
    int64
    0

    crs_wkt :

    GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]]

    semi_major_axis :

    6378137.0

    semi_minor_axis :

    6356752.314245179

    inverse_flattening :

    298.257223563

    reference_ellipsoid_name :

    WGS 84

    longitude_of_prime_meridian :

    0.0

    prime_meridian_name :

    Greenwich

    geographic_crs_name :

    WGS 84

    horizontal_datum_name :

    World Geodetic System 1984

    grid_mapping_name :

    latitude_longitude

    spatial_ref :

    GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AXIS["Latitude",NORTH],AXIS["Longitude",EAST],AUTHORITY["EPSG","4326"]]

    GeoTransform :

    -105.90416665717922 0.0020833333331466974 0.0 40.779166663013456 0.0 -0.0020833333331466974

    array(0)

  • date
    (date)
    datetime64[ns]
    2018-06-26 ... 2023-07-28

    array(['2018-06-26T00:00:00.000000000', '2018-07-12T00:00:00.000000000',
           '2018-07-28T00:00:00.000000000', '2019-06-26T00:00:00.000000000',
           '2019-07-12T00:00:00.000000000', '2019-07-28T00:00:00.000000000',
           '2020-06-25T00:00:00.000000000', '2020-07-11T00:00:00.000000000',
           '2020-07-27T00:00:00.000000000', '2021-06-26T00:00:00.000000000',
           '2021-07-12T00:00:00.000000000', '2021-07-28T00:00:00.000000000',
           '2022-06-26T00:00:00.000000000', '2022-07-12T00:00:00.000000000',
           '2022-07-28T00:00:00.000000000', '2023-06-26T00:00:00.000000000',
           '2023-07-12T00:00:00.000000000', '2023-07-28T00:00:00.000000000'],
          dtype='datetime64[ns]')

• Data variables: (1)
  • NDVI
    (date, y, x)
    float32
    0.5672 0.5711 ... 0.6327 0.6327

    array([[[0.5672, 0.5711, 0.6027, ..., 0.3888, 0.3888, 0.381 ],
            [0.5793, 0.6917, 0.6917, ..., 0.3902, 0.3683, 0.381 ],
            [0.6174, 0.6174, 0.6917, ..., 0.3867, 0.3902, 0.4094],
            ...,
            [0.2816, 0.2816, 0.1774, ..., 0.4708, 0.3777, 0.4174],
            [0.2966, 0.2966, 0.2701, ..., 0.3807, 0.4774, 0.4307],
            [0.363 , 0.2543, 0.2701, ..., 0.3878, 0.4797, 0.4797]],
    
           [[0.5806, 0.5789, 0.5865, ..., 0.3451, 0.3451, 0.3483],
            [0.5633, 0.6341, 0.6341, ..., 0.3777, 0.3484, 0.3483],
            [0.6144, 0.6144, 0.6309, ..., 0.3511, 0.3484, 0.3376],
            ...,
            [0.2836, 0.2836, 0.2292, ..., 0.4965, 0.4373, 0.4082],
            [0.249 , 0.2925, 0.3079, ..., 0.4385, 0.4865, 0.4492],
            [0.4257, 0.2689, 0.2137, ..., 0.4532, 0.4682, 0.4682]],
    
           [[0.6154, 0.5971, 0.5883, ..., 0.3806, 0.3806, 0.3809],
            [0.6154, 0.6525, 0.6525, ..., 0.4257, 0.363 , 0.3619],
            [0.5944, 0.5944, 0.6709, ..., 0.3582, 0.3541, 0.3617],
            ...,
    ...
            [0.1624, 0.1624, 0.0984, ..., 0.6586, 0.6443, 0.5771],
            [0.2924, 0.2428, 0.0984, ..., 0.659 , 0.7146, 0.678 ],
            [0.1628, 0.2924, 0.2267, ..., 0.659 , 0.7149, 0.7149]],
    
           [[0.6071, 0.5808, 0.5888, ..., 0.4637, 0.4637, 0.4901],
            [0.6071, 0.6501, 0.6501, ..., 0.4922, 0.4738, 0.4317],
            [0.6792, 0.6792, 0.6792, ..., 0.4516, 0.4334, 0.4317],
            ...,
            [0.213 , 0.213 , 0.1423, ..., 0.6524, 0.6861, 0.5645],
            [0.3889, 0.3428, 0.0974, ..., 0.6346, 0.686 , 0.6473],
            [0.4555, 0.3806, 0.3308, ..., 0.6316, 0.6616, 0.6616]],
    
           [[0.5787, 0.6146, 0.6053, ..., 0.427 , 0.427 , 0.427 ],
            [0.6303, 0.6303, 0.6303, ..., 0.4454, 0.4249, 0.4264],
            [0.6415, 0.6415, 0.6334, ..., 0.4316, 0.4117, 0.4594],
            ...,
            [0.24  , 0.24  , 0.2518, ..., 0.6657, 0.6761, 0.5383],
            [0.2111, 0.1931, 0.1973, ..., 0.6778, 0.7129, 0.603 ],
            [0.3925, 0.364 , 0.2705, ..., 0.6235, 0.6327, 0.6327]]],
          dtype=float32)

• Indexes: (3)
  • x
    PandasIndex

    PandasIndex(Index([-105.90312499051265,  -105.9010416571795, -105.89895832384636,
           -105.89687499051321, -105.89479165718006, -105.89270832384692,
           -105.89062499051377, -105.88854165718062, -105.88645832384748,
           -105.88437499051433,
           ...
           -105.21770832390739, -105.21562499057424,  -105.2135416572411,
           -105.21145832390795,  -105.2093749905748, -105.20729165724165,
           -105.20520832390851, -105.20312499057536, -105.20104165724221,
           -105.19895832390907],
          dtype='float64', name='x', length=339))

  • y
    PandasIndex

    PandasIndex(Index([ 40.77812499634688, 40.776041663013736,  40.77395832968059,
            40.77187499634744, 40.769791663014296,  40.76770832968115,
              40.765624996348, 40.763541663014855,  40.76145832968171,
            40.75937499634856,
           ...
           40.480208329706905,  40.47812499637376,  40.47604166304061,
           40.473958329707465,  40.47187499637432,  40.46979166304117,
           40.467708329708024,  40.46562499637488,  40.46354166304173,
           40.461458329708584],
          dtype='float64', name='y', length=153))

  • date
    PandasIndex

    PandasIndex(DatetimeIndex(['2018-06-26', '2018-07-12', '2018-07-28', '2019-06-26',
                   '2019-07-12', '2019-07-28', '2020-06-25', '2020-07-11',
                   '2020-07-27', '2021-06-26', '2021-07-12', '2021-07-28',
                   '2022-06-26', '2022-07-12', '2022-07-28', '2023-06-26',
                   '2023-07-12', '2023-07-28'],
                  dtype='datetime64[ns]', name='date', freq=None))

• Attributes: (0)

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.PiYG

There 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

# Compute the difference in NDVI before and after the fire

# Plot the difference
(
    ndvi_diff.hvplot(x='', y='', cmap='', geo=True)
    *
    gdf.hvplot(geo=True, fill_color=None, line_color='black')
)

NameError: name 'ndvi_diff' is not defined

ndvi_diff = (
    ndvi_da
        .sel(date=slice('2021', '2023'))
        .mean('date')
        .NDVI 
   - ndvi_da
        .sel(date=slice('2018', '2020'))
        .mean('date')
        .NDVI
)
(
    ndvi_diff.hvplot(x='x', y='y', cmap='PiYG', geo=True)
    *
    gdf.hvplot(geo=True, fill_color=None, line_color='black')
)

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.

out_gdf = (
    gpd.GeoDataFrame(geometry=gdf.envelope)
    .overlay(gdf, how='difference'))

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=True when 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

# Clip data to both inside and outside the boundary

ndvi_cp_da = ndvi_da.rio.clip(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 xarray what 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 the lsuffix= and rsuffix= parameters

GOTCHA ALERT

The DateIndex in pandas is a little different from the Datetime Dimension in xarray. You will need to use the .dt.year syntax 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.

:::

# Compute mean annual July NDVI

# Compute mean annual July NDVI
jul_ndvi_cp_df = (
    ndvi_cp_da
    .groupby(ndvi_cp_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_cp_df[['NDVI']]
    .join(
        jul_ndvi_out_df[['NDVI']], 
        lsuffix=' Burned Area', rsuffix=' Unburned Area')
)

jul_ndvi_df.hvplot(
    title='NDVI before and after the Cameron Peak Fire'
)

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!

# Plot difference inside and outside the reservation

# 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 Cameron Peak Fire'
)

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