NDVI after the Las Conchas fire in New Mexico¶
This notebook will be used to further investigate NDVI rasters before and after the Las Conchas fire in New Mexico. We'll also create plot with a slider so that we can see how NDVI recovers in the years after the fire.
Step 0: Restore variables and import libraries¶
In [1]:
%store -r ndvi_da fire_bound_gdf
In [2]:
# Import libraries
import earthpy
import hvplot.pandas
import hvplot.xarray
import geopandas as gpd
import matplotlib.pyplot as plt
import pandas as pd
import panel as pn
import rioxarray as rxr
import xarray as xr
from shapely.geometry import mapping
Step 1: Assess if NDVI is different inside and outside the fire boundary¶
In [3]:
# Check CRSs, make them match if needed.
if ndvi_da.NDVI.rio.crs == fire_bound_gdf.crs:
print('The CRSs match!')
else:
print(ndvi_da.NDVI.rio.crs) # EPSG:4326
print(fire_bound_gdf.crs) # EPSG:4269
# The CRSs match, so we don't need to do anything more here.
The CRSs match!
We need to clip the DataArray into two arrays, so that we can look at NDVI from inside the fire boundary, and outside.
In [4]:
# Clip data to both inside and outside the boundary
# Inner
ndvi_in = ndvi_da.rio.clip(fire_bound_gdf.geometry.apply(mapping))
# Outer
ndvi_out = ndvi_da.rio.clip(
fire_bound_gdf.geometry.apply(mapping), invert=True
)
In [5]:
# Quick plots
# We should see the two plots match the clipped area like puzzle pieces.
# Make plots
plot_in = ndvi_in.isel(date=0).hvplot(x='x', y='y', cmap=plt.cm.PiYG, geo=True,
title='NDVI inside Fire Boundary')
plot_out = ndvi_out.isel(date=0).hvplot(x='x', y='y', cmap=plt.cm.PiYG, geo=True,
title='NDVI outside Fire Boundary')
# Show plots
(plot_in + plot_out)
Out[5]:
In [6]:
# Compute mean annual NDVI
# Inside
ndvi_in_mean_df = (ndvi_in
#.sel(date=ndvi_in['date'].dt.month.isin([7]))
# Select NDVI data in July; I've commented this out so the
# analysis works
# Now group into years
.groupby('date.year')
# Now take the mean across all dimensions
.mean(...)
.to_dataframe()
)
# Outside
ndvi_out_mean_df = (
ndvi_out#.sel(date=ndvi_out['date'].dt.month.isin([7]))
.groupby('date.year')
.mean(...)
.to_dataframe()
)
In [7]:
# Join the inside and outside DFs
ndvi_df = (ndvi_in_mean_df[['NDVI']]
.join(ndvi_out_mean_df[['NDVI']],
lsuffix=' Inside Fire',
rsuffix=' Outside Fire'))
ndvi_df
Out[7]:
| NDVI Inside Fire | NDVI Outside Fire | |
|---|---|---|
| year | ||
| 2006 | 0.536539 | 0.470742 |
| 2007 | 0.575677 | 0.502179 |
| 2008 | 0.560684 | 0.486118 |
| 2009 | 0.569864 | 0.495966 |
| 2010 | 0.576190 | 0.491205 |
| 2011 | 0.426936 | 0.458000 |
| 2012 | 0.377027 | 0.455526 |
| 2013 | 0.390443 | 0.447340 |
| 2014 | 0.439131 | 0.479393 |
| 2015 | 0.472262 | 0.500086 |
| 2016 | 0.467678 | 0.486759 |
| 2017 | 0.473199 | 0.486518 |
| 2018 | 0.435072 | 0.445602 |
| 2019 | 0.484973 | 0.487404 |
| 2020 | 0.440681 | 0.445364 |
| 2021 | 0.454713 | 0.460829 |
In [8]:
# Plot annual mean NDVI inside and outside on one plot
ins_out_plot_ma_plot = ndvi_df.hvplot(
by='column', color=['red', 'green'],
# set inside to red, outside to blue
title='Mean Annual NDVI inside and outside\n' \
'Las Conchas Fire Boundary')
# Save plot
hvplot.save(ins_out_plot_ma_plot, 'Mean_NDVI_In_Out_Fire.html')
# Show the plot
ins_out_plot_ma_plot
Out[8]:
Now, we'll plot just the difference in NDVI between the two areas on an annual basis.
In [9]:
# Calculate the difference inside and outside
ndvi_df['Difference'] = ndvi_df['NDVI Outside Fire'] - ndvi_df['NDVI Inside Fire']
ndvi_df
Out[9]:
| NDVI Inside Fire | NDVI Outside Fire | Difference | |
|---|---|---|---|
| year | |||
| 2006 | 0.536539 | 0.470742 | -0.065797 |
| 2007 | 0.575677 | 0.502179 | -0.073498 |
| 2008 | 0.560684 | 0.486118 | -0.074566 |
| 2009 | 0.569864 | 0.495966 | -0.073898 |
| 2010 | 0.576190 | 0.491205 | -0.084985 |
| 2011 | 0.426936 | 0.458000 | 0.031064 |
| 2012 | 0.377027 | 0.455526 | 0.078499 |
| 2013 | 0.390443 | 0.447340 | 0.056897 |
| 2014 | 0.439131 | 0.479393 | 0.040262 |
| 2015 | 0.472262 | 0.500086 | 0.027824 |
| 2016 | 0.467678 | 0.486759 | 0.019082 |
| 2017 | 0.473199 | 0.486518 | 0.013319 |
| 2018 | 0.435072 | 0.445602 | 0.010529 |
| 2019 | 0.484973 | 0.487404 | 0.002431 |
| 2020 | 0.440681 | 0.445364 | 0.004683 |
| 2021 | 0.454713 | 0.460829 | 0.006116 |
In [10]:
# Plot annual mean NDVI inside and outside on one plot
difference_ma_plot = ndvi_df['Difference'].hvplot(
title=('Difference in Mean Annual NDVI\n'
'outside and inside Las Conchas Fire Boundary'),
ylabel='NDVI')
# Save plot
hvplot.save(difference_ma_plot, 'NDVI_Mean_Diff_line.html')
# Show the plot
difference_ma_plot
Out[10]:
PLOT INTERPRETATION HERE
Step 2: Create a dynamic plot to see how NDVI changes each year after the fire.¶
In [11]:
# Compute annual spatial means
ndvi_mean_da = (
ndvi_da
.groupby('date.year')
.mean(dim='date')
)
ndvi_mean_da
Out[11]:
<xarray.Dataset> Size: 2MB
Dimensions: (year: 16, y: 197, x: 137)
Coordinates:
band int64 8B 1
* x (x) float64 1kB -106.6 -106.6 -106.6 ... -106.3 -106.3 -106.3
* y (y) float64 2kB 36.07 36.07 36.07 36.07 ... 35.67 35.67 35.67
spatial_ref int64 8B 0
* year (year) int64 128B 2006 2007 2008 2009 ... 2018 2019 2020 2021
Data variables:
NDVI (year, y, x) float32 2MB 0.5589 0.5676 0.5676 ... 0.259 0.2621In [12]:
# Set up a DataArray for visualization with a slider
# Set the first slice of the array to a baseline NDVI
# We'll average the years 2006-2010
ndvi_base_da = ndvi_mean_da.sel(year = slice(2006, 2010)).mean("year")
# Reassign a year coordinate to the DA so we can concatenate later
ndvi_base_da = ndvi_base_da.expand_dims({"year": ["2006-2010"]})
# Grab each year after 2010
ndvi_post_da = ndvi_mean_da.sel(year = slice(2011, 2021))
# Now combine the two DAs
ndvi_vis_da = xr.concat([ndvi_base_da, ndvi_post_da], dim = 'year')
# Check it out
ndvi_vis_da
Out[12]:
<xarray.Dataset> Size: 1MB
Dimensions: (year: 12, y: 197, x: 137)
Coordinates:
* year (year) object 96B '2006-2010' 2011 2012 2013 ... 2019 2020 2021
band int64 8B 1
* x (x) float64 1kB -106.6 -106.6 -106.6 ... -106.3 -106.3 -106.3
* y (y) float64 2kB 36.07 36.07 36.07 36.07 ... 35.67 35.67 35.67
spatial_ref int64 8B 0
Data variables:
NDVI (year, y, x) float32 1MB 0.5961 0.5952 0.5952 ... 0.259 0.2621In [13]:
ndvi_vis_da.year
Out[13]:
<xarray.DataArray 'year' (year: 12)> Size: 96B
array(['2006-2010', 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018, 2019, 2020,
2021], dtype=object)
Coordinates:
* year (year) object 96B '2006-2010' 2011 2012 2013 ... 2019 2020 2021
band int64 8B 1
spatial_ref int64 8B 0In [14]:
ndvi_vis_da.year.values.astype(str)
Out[14]:
array(['2006-2010', '2011', '2012', '2013', '2014', '2015', '2016',
'2017', '2018', '2019', '2020', '2021'], dtype='<U9')
This DataArray looks like what we need. Now, we'll make a plot with a dynamic slider.
In [15]:
# Set the year coordinate to string so that our slider works
ndvi_vis_da = ndvi_vis_da.assign_coords(
year=ndvi_vis_da.year.astype(str)
)
# Plot interactive baseline NDVI and post fire years
ndvi_plot = ndvi_vis_da.hvplot(
x='x', y='y',
#groupby='year',
geo=True,
cmap=plt.cm.PiYG,
#widget_location = 'bottom',
title='2006-2010 Baseline and Post-Fire Recovery NDVI',
xlabel='Longitude',
ylabel='Latitude'
)
# Overlay the fire boundary
ndvi_int_plot = ndvi_plot * fire_bound_gdf.hvplot(
geo=True,
fill_color=None,
line_color='black')
# Save plot to be put on the website
panel_plot = pn.panel(ndvi_int_plot)
panel_plot.save('NDVI_Recovery.html', embed=True)
# This is going to look a bit clunky on the website,
# but it's the best we can do.
ndvi_int_plot
# I've used ChatGPT to help me troubleshoot this plotting.
WARNING:bokeh.core.validation.check:W-1005 (FIXED_SIZING_MODE): 'fixed' sizing mode requires width and height to be set: figure(id='p1579', ...)
Out[15]:
BokehModel(combine_events=True, render_bundle={'docs_json': {'60a1d1fc-a1f3-439e-b3ca-bead1fddc5cb': {'version…