Analysis of Vegetation Health in Rocky Mountain National Park Using NDVI Data from 2000 to 2024¶
Introduction¶
Rocky Mountain National Park (RMNP) in Colorado is known for its impressive mountains, forests, alpine landscapes, and wildlife. While there are many important ecological features in the park, pine forests are especially important for a variety of reasons. Pine trees in the park help to reduce erosion, regulate temperatures, and provide important habitat for many birds and animals (Colorado State Forest Service). However, in recent decades, these trees have been under threat from pine bark beetles and wildfires. To analyse the impact of these threats, I will be using the Normalized Difference Vegetation Index (NDVI) from the timeframe of 2000 to 2024. This index measures "greenness" of vegetation using remotely sensed data, which helps to quantify vegetation health.
Image Source: USGS. RMNP. https://www.usgs.gov/geology-and-ecology-of-national-parks/geology-rocky-mountain-national-park
Pine Bark Beetles in Rocky Mountain National Park¶
According to the National Park Service, there are 17 native species of pine bark beetles common in the park, and for decades, periodic outbreaks were relatively common but of small scale. Pine bark beetle populations are affected by winter temperatures, with consistent periods of low temperatures needed to kill off beetles and prevent their populations for growing too rapidly. However, warmer average temperatures in the park in recent decades have allowed beetle populations to dramatically increase, with a peak in severity around 2009 when a large percentage of trees in the park were killed (Mountain Pine Beetle. NPS). The pine beetles kill lodgepole pine, ponderosa pine, limber pine, Engelmann spruce, subalpine fir, and Colorado blue spruce (Mountain Pine Beetle. NPS).
Wildfires in Rocky Mountain National Park¶
Fires are a natural part of forest ecosystems. However, with impacts from climate change and pine bark beetles, fires in the park and surrounding area have been increasing in severity and size in recent years (Fire History. NPS). In 2012, the Fern Lake Fire burned over 3,500 acres in the park, with beetle kill trees exacerbating conditions. In 2020, the Cameron Peak and East Troublesome Fires burned roughly 30,000 acres within the park (Fire History. NPS).
Timeframe¶
To analyze the impacts of pine beetles and wildfires, I will be using NDVI data from 2000 to 2024. This time frame will encompass the beginnings of the pine beetle outbreaks starting in the early 2000s through 2009, and it will also capture the more recent wildfires from 2012 and 2020, as well as any recent regrowth from these impacts. Additionally, I'm only interested in the vegetation health of pine trees, so I will be downloading data from January to account for this. In January, the only vegetation cover should be pine trees, with other vegetation either covered in snow or dormant.
Data¶
I will be using Moderate Resolution Imaging Spectroradiometer (MODIS) NDVI data downloaded from NASA's EarthData catalog. Specifically, I will be downloading via NASA's Application for Extracting and Exploring Analysis Ready Samples (AppEEARS) using the earthpy package. This will allow me to access the data quickly in Github Codespaces.
Python Workflow¶
Below is the Python workflow that I'm using for this analysis.
# Import packages
# File structure
import os
import pathlib
# Vector data
import geopandas as gpd
import pandas as pd
# Raster data
import rioxarray as rxr
import xarray as xr
# Arrays
import numpy as np
# Maps and plots
import holoviews as hv
import hvplot.pandas
import hvplot.xarray
# Open street map
from osmnx import features as osm
import osmnx as ox
# Earthpy
import earthpy
import earthpy.api.appeears as eaapp
# Plotting
import matplotlib.pyplot as plt
from sklearn.linear_model import LinearRegression
# Set up a project and folder for data using earthpy
project = earthpy.Project("Gila River Vegetation", dirname = 'ndvi_project')
# Search for RMNP using OpenStreetMaps
rmnp_gdf = ox.geocode_to_gdf(
'Rocky Mountain National Park')
# Check it out
rmnp_gdf
| geometry | bbox_west | bbox_south | bbox_east | bbox_north | place_id | osm_type | osm_id | lat | lon | class | type | place_rank | importance | addresstype | name | display_name | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | MULTIPOLYGON (((-105.53362 40.29996, -105.5338... | -105.913714 | 40.15777 | -105.493583 | 40.553787 | 416089268 | relation | 390960 | 40.355392 | -105.717696 | boundary | protected_area | 25 | 0.520509 | protected_area | Rocky Mountain National Park | Rocky Mountain National Park, Larimer County, ... |
# Make a plot of the boundary
rmnp_gdf.plot()
<Axes: >
Figure 1: This plot is showing the general footprint of the park. Data from OpenStreetMaps.
# Initialize AppeearsDownloader for MODIS NDVI data
# Set parameters
ndvi_downloader = eaapp.AppeearsDownloader(
# Give the download a name
download_key = "rmnp_ndvi",
# Tell it to put the data in the project directory already defined
project = project,
# Specify the MODIS product
product = 'MOD13Q1.061',
layer = '_250m_16_days_NDVI',
# Set a start date and end data
# I'm using January for reasons described above
start_date = "01-01",
end_date = "01-30",
# Recurring meaning we want to download those dates over multiple years
recurring = True,
# Specify the range of years of interest
year_range = [2000, 2024],
# Specify the polygon we want to get NDVI data for
polygon = rmnp_gdf
)
# Download the prepared download
ndvi_downloader.download_files(cache=True)
<generator object Path.rglob at 0x72e32d737450>
# Get a sorted list of NDVI file paths
ndvi_paths = sorted(list(project.project_dir.rglob('*NDVI*.tif')))
# View the paths
ndvi_paths
[PosixPath('/workspaces/data/ndvi_project/rmnp_ndvi/MOD13Q1.061_1999351_to_2024017/MOD13Q1.061__250m_16_days_NDVI_doy2000353000000_aid0001.tif'),
PosixPath('/workspaces/data/ndvi_project/rmnp_ndvi/MOD13Q1.061_1999351_to_2024017/MOD13Q1.061__250m_16_days_NDVI_doy2001001000000_aid0001.tif'),
PosixPath('/workspaces/data/ndvi_project/rmnp_ndvi/MOD13Q1.061_1999351_to_2024017/MOD13Q1.061__250m_16_days_NDVI_doy2001017000000_aid0001.tif'),
PosixPath('/workspaces/data/ndvi_project/rmnp_ndvi/MOD13Q1.061_1999351_to_2024017/MOD13Q1.061__250m_16_days_NDVI_doy2001353000000_aid0001.tif'),
PosixPath('/workspaces/data/ndvi_project/rmnp_ndvi/MOD13Q1.061_1999351_to_2024017/MOD13Q1.061__250m_16_days_NDVI_doy2002001000000_aid0001.tif'),
PosixPath('/workspaces/data/ndvi_project/rmnp_ndvi/MOD13Q1.061_1999351_to_2024017/MOD13Q1.061__250m_16_days_NDVI_doy2002017000000_aid0001.tif'),
PosixPath('/workspaces/data/ndvi_project/rmnp_ndvi/MOD13Q1.061_1999351_to_2024017/MOD13Q1.061__250m_16_days_NDVI_doy2002353000000_aid0001.tif'),
PosixPath('/workspaces/data/ndvi_project/rmnp_ndvi/MOD13Q1.061_1999351_to_2024017/MOD13Q1.061__250m_16_days_NDVI_doy2003001000000_aid0001.tif'),
PosixPath('/workspaces/data/ndvi_project/rmnp_ndvi/MOD13Q1.061_1999351_to_2024017/MOD13Q1.061__250m_16_days_NDVI_doy2003017000000_aid0001.tif'),
PosixPath('/workspaces/data/ndvi_project/rmnp_ndvi/MOD13Q1.061_1999351_to_2024017/MOD13Q1.061__250m_16_days_NDVI_doy2003353000000_aid0001.tif'),
PosixPath('/workspaces/data/ndvi_project/rmnp_ndvi/MOD13Q1.061_1999351_to_2024017/MOD13Q1.061__250m_16_days_NDVI_doy2004001000000_aid0001.tif'),
PosixPath('/workspaces/data/ndvi_project/rmnp_ndvi/MOD13Q1.061_1999351_to_2024017/MOD13Q1.061__250m_16_days_NDVI_doy2004017000000_aid0001.tif'),
PosixPath('/workspaces/data/ndvi_project/rmnp_ndvi/MOD13Q1.061_1999351_to_2024017/MOD13Q1.061__250m_16_days_NDVI_doy2004353000000_aid0001.tif'),
PosixPath('/workspaces/data/ndvi_project/rmnp_ndvi/MOD13Q1.061_1999351_to_2024017/MOD13Q1.061__250m_16_days_NDVI_doy2005001000000_aid0001.tif'),
PosixPath('/workspaces/data/ndvi_project/rmnp_ndvi/MOD13Q1.061_1999351_to_2024017/MOD13Q1.061__250m_16_days_NDVI_doy2005017000000_aid0001.tif'),
PosixPath('/workspaces/data/ndvi_project/rmnp_ndvi/MOD13Q1.061_1999351_to_2024017/MOD13Q1.061__250m_16_days_NDVI_doy2005353000000_aid0001.tif'),
PosixPath('/workspaces/data/ndvi_project/rmnp_ndvi/MOD13Q1.061_1999351_to_2024017/MOD13Q1.061__250m_16_days_NDVI_doy2006001000000_aid0001.tif'),
PosixPath('/workspaces/data/ndvi_project/rmnp_ndvi/MOD13Q1.061_1999351_to_2024017/MOD13Q1.061__250m_16_days_NDVI_doy2006017000000_aid0001.tif'),
PosixPath('/workspaces/data/ndvi_project/rmnp_ndvi/MOD13Q1.061_1999351_to_2024017/MOD13Q1.061__250m_16_days_NDVI_doy2006353000000_aid0001.tif'),
PosixPath('/workspaces/data/ndvi_project/rmnp_ndvi/MOD13Q1.061_1999351_to_2024017/MOD13Q1.061__250m_16_days_NDVI_doy2007001000000_aid0001.tif'),
PosixPath('/workspaces/data/ndvi_project/rmnp_ndvi/MOD13Q1.061_1999351_to_2024017/MOD13Q1.061__250m_16_days_NDVI_doy2007017000000_aid0001.tif'),
PosixPath('/workspaces/data/ndvi_project/rmnp_ndvi/MOD13Q1.061_1999351_to_2024017/MOD13Q1.061__250m_16_days_NDVI_doy2007353000000_aid0001.tif'),
PosixPath('/workspaces/data/ndvi_project/rmnp_ndvi/MOD13Q1.061_1999351_to_2024017/MOD13Q1.061__250m_16_days_NDVI_doy2008001000000_aid0001.tif'),
PosixPath('/workspaces/data/ndvi_project/rmnp_ndvi/MOD13Q1.061_1999351_to_2024017/MOD13Q1.061__250m_16_days_NDVI_doy2008017000000_aid0001.tif'),
PosixPath('/workspaces/data/ndvi_project/rmnp_ndvi/MOD13Q1.061_1999351_to_2024017/MOD13Q1.061__250m_16_days_NDVI_doy2008353000000_aid0001.tif'),
PosixPath('/workspaces/data/ndvi_project/rmnp_ndvi/MOD13Q1.061_1999351_to_2024017/MOD13Q1.061__250m_16_days_NDVI_doy2009001000000_aid0001.tif'),
PosixPath('/workspaces/data/ndvi_project/rmnp_ndvi/MOD13Q1.061_1999351_to_2024017/MOD13Q1.061__250m_16_days_NDVI_doy2009017000000_aid0001.tif'),
PosixPath('/workspaces/data/ndvi_project/rmnp_ndvi/MOD13Q1.061_1999351_to_2024017/MOD13Q1.061__250m_16_days_NDVI_doy2009353000000_aid0001.tif'),
PosixPath('/workspaces/data/ndvi_project/rmnp_ndvi/MOD13Q1.061_1999351_to_2024017/MOD13Q1.061__250m_16_days_NDVI_doy2010001000000_aid0001.tif'),
PosixPath('/workspaces/data/ndvi_project/rmnp_ndvi/MOD13Q1.061_1999351_to_2024017/MOD13Q1.061__250m_16_days_NDVI_doy2010017000000_aid0001.tif'),
PosixPath('/workspaces/data/ndvi_project/rmnp_ndvi/MOD13Q1.061_1999351_to_2024017/MOD13Q1.061__250m_16_days_NDVI_doy2010353000000_aid0001.tif'),
PosixPath('/workspaces/data/ndvi_project/rmnp_ndvi/MOD13Q1.061_1999351_to_2024017/MOD13Q1.061__250m_16_days_NDVI_doy2011001000000_aid0001.tif'),
PosixPath('/workspaces/data/ndvi_project/rmnp_ndvi/MOD13Q1.061_1999351_to_2024017/MOD13Q1.061__250m_16_days_NDVI_doy2011017000000_aid0001.tif'),
PosixPath('/workspaces/data/ndvi_project/rmnp_ndvi/MOD13Q1.061_1999351_to_2024017/MOD13Q1.061__250m_16_days_NDVI_doy2011353000000_aid0001.tif'),
PosixPath('/workspaces/data/ndvi_project/rmnp_ndvi/MOD13Q1.061_1999351_to_2024017/MOD13Q1.061__250m_16_days_NDVI_doy2012001000000_aid0001.tif'),
PosixPath('/workspaces/data/ndvi_project/rmnp_ndvi/MOD13Q1.061_1999351_to_2024017/MOD13Q1.061__250m_16_days_NDVI_doy2012017000000_aid0001.tif'),
PosixPath('/workspaces/data/ndvi_project/rmnp_ndvi/MOD13Q1.061_1999351_to_2024017/MOD13Q1.061__250m_16_days_NDVI_doy2012353000000_aid0001.tif'),
PosixPath('/workspaces/data/ndvi_project/rmnp_ndvi/MOD13Q1.061_1999351_to_2024017/MOD13Q1.061__250m_16_days_NDVI_doy2013001000000_aid0001.tif'),
PosixPath('/workspaces/data/ndvi_project/rmnp_ndvi/MOD13Q1.061_1999351_to_2024017/MOD13Q1.061__250m_16_days_NDVI_doy2013017000000_aid0001.tif'),
PosixPath('/workspaces/data/ndvi_project/rmnp_ndvi/MOD13Q1.061_1999351_to_2024017/MOD13Q1.061__250m_16_days_NDVI_doy2013353000000_aid0001.tif'),
PosixPath('/workspaces/data/ndvi_project/rmnp_ndvi/MOD13Q1.061_1999351_to_2024017/MOD13Q1.061__250m_16_days_NDVI_doy2014001000000_aid0001.tif'),
PosixPath('/workspaces/data/ndvi_project/rmnp_ndvi/MOD13Q1.061_1999351_to_2024017/MOD13Q1.061__250m_16_days_NDVI_doy2014017000000_aid0001.tif'),
PosixPath('/workspaces/data/ndvi_project/rmnp_ndvi/MOD13Q1.061_1999351_to_2024017/MOD13Q1.061__250m_16_days_NDVI_doy2014353000000_aid0001.tif'),
PosixPath('/workspaces/data/ndvi_project/rmnp_ndvi/MOD13Q1.061_1999351_to_2024017/MOD13Q1.061__250m_16_days_NDVI_doy2015001000000_aid0001.tif'),
PosixPath('/workspaces/data/ndvi_project/rmnp_ndvi/MOD13Q1.061_1999351_to_2024017/MOD13Q1.061__250m_16_days_NDVI_doy2015017000000_aid0001.tif'),
PosixPath('/workspaces/data/ndvi_project/rmnp_ndvi/MOD13Q1.061_1999351_to_2024017/MOD13Q1.061__250m_16_days_NDVI_doy2015353000000_aid0001.tif'),
PosixPath('/workspaces/data/ndvi_project/rmnp_ndvi/MOD13Q1.061_1999351_to_2024017/MOD13Q1.061__250m_16_days_NDVI_doy2016001000000_aid0001.tif'),
PosixPath('/workspaces/data/ndvi_project/rmnp_ndvi/MOD13Q1.061_1999351_to_2024017/MOD13Q1.061__250m_16_days_NDVI_doy2016017000000_aid0001.tif'),
PosixPath('/workspaces/data/ndvi_project/rmnp_ndvi/MOD13Q1.061_1999351_to_2024017/MOD13Q1.061__250m_16_days_NDVI_doy2016353000000_aid0001.tif'),
PosixPath('/workspaces/data/ndvi_project/rmnp_ndvi/MOD13Q1.061_1999351_to_2024017/MOD13Q1.061__250m_16_days_NDVI_doy2017001000000_aid0001.tif'),
PosixPath('/workspaces/data/ndvi_project/rmnp_ndvi/MOD13Q1.061_1999351_to_2024017/MOD13Q1.061__250m_16_days_NDVI_doy2017017000000_aid0001.tif'),
PosixPath('/workspaces/data/ndvi_project/rmnp_ndvi/MOD13Q1.061_1999351_to_2024017/MOD13Q1.061__250m_16_days_NDVI_doy2017353000000_aid0001.tif'),
PosixPath('/workspaces/data/ndvi_project/rmnp_ndvi/MOD13Q1.061_1999351_to_2024017/MOD13Q1.061__250m_16_days_NDVI_doy2018001000000_aid0001.tif'),
PosixPath('/workspaces/data/ndvi_project/rmnp_ndvi/MOD13Q1.061_1999351_to_2024017/MOD13Q1.061__250m_16_days_NDVI_doy2018017000000_aid0001.tif'),
PosixPath('/workspaces/data/ndvi_project/rmnp_ndvi/MOD13Q1.061_1999351_to_2024017/MOD13Q1.061__250m_16_days_NDVI_doy2018353000000_aid0001.tif'),
PosixPath('/workspaces/data/ndvi_project/rmnp_ndvi/MOD13Q1.061_1999351_to_2024017/MOD13Q1.061__250m_16_days_NDVI_doy2019001000000_aid0001.tif'),
PosixPath('/workspaces/data/ndvi_project/rmnp_ndvi/MOD13Q1.061_1999351_to_2024017/MOD13Q1.061__250m_16_days_NDVI_doy2019017000000_aid0001.tif'),
PosixPath('/workspaces/data/ndvi_project/rmnp_ndvi/MOD13Q1.061_1999351_to_2024017/MOD13Q1.061__250m_16_days_NDVI_doy2019353000000_aid0001.tif'),
PosixPath('/workspaces/data/ndvi_project/rmnp_ndvi/MOD13Q1.061_1999351_to_2024017/MOD13Q1.061__250m_16_days_NDVI_doy2020001000000_aid0001.tif'),
PosixPath('/workspaces/data/ndvi_project/rmnp_ndvi/MOD13Q1.061_1999351_to_2024017/MOD13Q1.061__250m_16_days_NDVI_doy2020017000000_aid0001.tif'),
PosixPath('/workspaces/data/ndvi_project/rmnp_ndvi/MOD13Q1.061_1999351_to_2024017/MOD13Q1.061__250m_16_days_NDVI_doy2020353000000_aid0001.tif'),
PosixPath('/workspaces/data/ndvi_project/rmnp_ndvi/MOD13Q1.061_1999351_to_2024017/MOD13Q1.061__250m_16_days_NDVI_doy2021001000000_aid0001.tif'),
PosixPath('/workspaces/data/ndvi_project/rmnp_ndvi/MOD13Q1.061_1999351_to_2024017/MOD13Q1.061__250m_16_days_NDVI_doy2021017000000_aid0001.tif'),
PosixPath('/workspaces/data/ndvi_project/rmnp_ndvi/MOD13Q1.061_1999351_to_2024017/MOD13Q1.061__250m_16_days_NDVI_doy2021353000000_aid0001.tif'),
PosixPath('/workspaces/data/ndvi_project/rmnp_ndvi/MOD13Q1.061_1999351_to_2024017/MOD13Q1.061__250m_16_days_NDVI_doy2022001000000_aid0001.tif'),
PosixPath('/workspaces/data/ndvi_project/rmnp_ndvi/MOD13Q1.061_1999351_to_2024017/MOD13Q1.061__250m_16_days_NDVI_doy2022017000000_aid0001.tif'),
PosixPath('/workspaces/data/ndvi_project/rmnp_ndvi/MOD13Q1.061_1999351_to_2024017/MOD13Q1.061__250m_16_days_NDVI_doy2022353000000_aid0001.tif'),
PosixPath('/workspaces/data/ndvi_project/rmnp_ndvi/MOD13Q1.061_1999351_to_2024017/MOD13Q1.061__250m_16_days_NDVI_doy2023001000000_aid0001.tif'),
PosixPath('/workspaces/data/ndvi_project/rmnp_ndvi/MOD13Q1.061_1999351_to_2024017/MOD13Q1.061__250m_16_days_NDVI_doy2023017000000_aid0001.tif'),
PosixPath('/workspaces/data/ndvi_project/rmnp_ndvi/MOD13Q1.061_1999351_to_2024017/MOD13Q1.061__250m_16_days_NDVI_doy2023353000000_aid0001.tif'),
PosixPath('/workspaces/data/ndvi_project/rmnp_ndvi/MOD13Q1.061_1999351_to_2024017/MOD13Q1.061__250m_16_days_NDVI_doy2024001000000_aid0001.tif'),
PosixPath('/workspaces/data/ndvi_project/rmnp_ndvi/MOD13Q1.061_1999351_to_2024017/MOD13Q1.061__250m_16_days_NDVI_doy2024017000000_aid0001.tif')]
# Plot the results with ESRI web tile images
rmnp_gdf.hvplot(
geo=True,
tiles='EsriImagery',
fill_color=None,
line_color='red',
line_width=2,
title='Rocky Mountain National Park Boundary',
frame_width=500)
Figure 2: This map shows the park boundary without fill color, and with an ESRI satellite imagery background layer to show the general landscape and area of the park.
# Extracting the string of date values from the file names
# Including the .tif in character count
doy_start = -25 # Counting backwards to the first character of the date string
doy_end = -18 # Counting backwards to the last character of the date string
# Loop through each NDVI image
ndvi_das = []
for ndvi_path in ndvi_paths:
# Get date from file name
doy = ndvi_path.name[doy_start:doy_end] #Extract the date
date = pd.to_datetime(doy, format='%Y%j') #Turn it into pandas date-time object
# Open the NDVI dataset (da = data array)
da = rxr.open_rasterio(ndvi_path, mask_and_scale=True).squeeze()
# Add date dimension and clean up metadata
da = da.assign_coords({'date': date})
da = da.expand_dims({'date': 1})
da.name = 'NDVI'
# Prepare for concatenation
ndvi_das.append(da)
len(ndvi_das) # Check the number of NDVI files in ndvi_das
72
This code chunk above is adding a pandas date time object to each NDVI file, since the raw data just has the date in the file name.
# Combine NDVI images from all dates
# This gives us two arrays (within ndvi_da) with cooresponding dates and NDVI values
ndvi_da = xr.combine_by_coords(ndvi_das, coords=['date'])
ndvi_da
/tmp/ipykernel_1496/1185387115.py:3: FutureWarning: In a future version of xarray the default value for compat will change from compat='no_conflicts' to compat='override'. This is likely to lead to different results when combining overlapping variables with the same name. To opt in to new defaults and get rid of these warnings now use `set_options(use_new_combine_kwarg_defaults=True) or set compat explicitly. ndvi_da = xr.combine_by_coords(ndvi_das, coords=['date']) /tmp/ipykernel_1496/1185387115.py:3: FutureWarning: In a future version of xarray the default value for compat will change from compat='no_conflicts' to compat='override'. This is likely to lead to different results when combining overlapping variables with the same name. To opt in to new defaults and get rid of these warnings now use `set_options(use_new_combine_kwarg_defaults=True) or set compat explicitly. ndvi_da = xr.combine_by_coords(ndvi_das, coords=['date'])
<xarray.Dataset> Size: 11MB
Dimensions: (date: 72, y: 191, x: 203)
Coordinates:
band int64 8B 1
* x (x) float64 2kB -105.9 -105.9 -105.9 ... -105.5 -105.5 -105.5
* y (y) float64 2kB 40.55 40.55 40.55 40.55 ... 40.16 40.16 40.16
spatial_ref int64 8B 0
* date (date) datetime64[ns] 576B 2000-12-18 2001-01-01 ... 2024-01-17
Data variables:
NDVI (date, y, x) float32 11MB 0.1285 0.1285 ... 0.4633 0.4343# Plot the 1st tif file from ndvi_paths
first_image = rxr.open_rasterio(ndvi_paths[0], mask_and_scale=True).squeeze().plot()
first_image.axes.set_title("RMNP NDVI image 1 from March, 2000")
Text(0.5, 1.0, 'RMNP NDVI image 1 from March, 2000')
Figure 3: This is a test visualization of the first NDVI image downloaded. Higher values (dark red) represent healthier/denser vegetation. Lower values (white, light pink, or blue) represent less healthy/dense vegetation.
# Plot 1st and last NDVI images side-by-side for comparison
old_ndvi = rxr.open_rasterio(ndvi_paths[0], mask_and_scale=True).squeeze()
recent_ndvi = rxr.open_rasterio(ndvi_paths[-1], mask_and_scale=True).squeeze()
# Create side-by-side subplots
# This part was created thanks to the video demonstration
fig, axes = plt.subplots(1 , 2, figsize=[12, 4])
# Plot each in their own axes. PiYG is a value for the green/pink color scheme
old_ndvi.plot(ax=axes[0], cmap=plt.cm.PiYG, vmin=-1, vmax=1)
rmnp_gdf.plot(ax=axes[0], edgecolor='black', facecolor='none', linewidth=1)
axes[0].set_title("NDVI - RMNP 2000")
# Ploting the recent ndvi image using axis 1
recent_ndvi.plot(ax=axes[1], cmap=plt.cm.PiYG, vmin=-1, vmax=1)
rmnp_gdf.plot(ax=axes[1], edgecolor='black', facecolor='none', linewidth=1)
axes[1].set_title("NDVI - RMNP 2024")
Text(0.5, 1.0, 'NDVI - RMNP 2024')
Figure 4: These side-by-side maps compare NDVI in 2000 vs. 2024. Darker greens represent healthier vegetation, whereas whites/pinks represent less healthy vegetation.
Calculating a Difference Map¶
# Compute the difference in NDVI from the first and second half of the study period
ndvi_diff = (
ndvi_da
.sel(date=slice('2012', '2024'))
.mean('date')
.NDVI
- ndvi_da
.sel(date=slice('2000','2012'))
.mean('date')
.NDVI
)
ndvi_diff
# Here we are subtracting the mean NDVI of 2000-2012 from the mean NDVI of 2012-2024.
<xarray.DataArray 'NDVI' (y: 191, x: 203)> Size: 155kB
array([[ 7.3079765e-03, 1.5445661e-02, 1.1941388e-02, ...,
-1.8431789e-01, -1.8790740e-01, -1.9800520e-01],
[-1.2013502e-03, 7.6419674e-03, 8.8908169e-03, ...,
-2.0381254e-01, -1.6738772e-01, -1.8971199e-01],
[-1.4498785e-02, 1.5648007e-03, 1.5648007e-03, ...,
-2.4673939e-01, -1.9475049e-01, -1.9475049e-01],
...,
[-1.6470850e-03, 6.3335001e-03, 7.8702755e-03, ...,
5.3730220e-02, 4.2035460e-02, 2.3916066e-02],
[ 1.5147030e-04, 2.7697831e-03, 8.8216420e-03, ...,
1.9586712e-02, 5.4537535e-02, 5.4537535e-02],
[-3.3464432e-03, 1.5298715e-03, 1.5298715e-03, ...,
1.7007262e-02, 1.7007262e-02, 2.9352933e-02]],
shape=(191, 203), dtype=float32)
Coordinates:
band int64 8B 1
* x (x) float64 2kB -105.9 -105.9 -105.9 ... -105.5 -105.5 -105.5
* y (y) float64 2kB 40.55 40.55 40.55 40.55 ... 40.16 40.16 40.16
spatial_ref int64 8B 0# Plot the difference
(
ndvi_diff.hvplot(
x='x',
y='y',
cmap='PiYG',
geo=True,
title="NDVI Changes in Rocky Mountain National Park \n2012-2024 vs. 2000-2012")
* # Adding the boundary plot from earlier
rmnp_gdf.hvplot(
geo=True,
fill_color=None,
line_color='black',
frame_width=500)
)
Figure 5: This map shows change in NDVI values for the two date slices. I subtracted the earlier date slice (2000-2012) from the more recent date slice (2012-2024). The result shows large areas of pink which indicates a decrease in vegetation health over the time period. There are also some green areas which indicate an increase in vegetation health over the time period. Overall, it's not surprising to see a decrease in NDVI values due to the impacts from pine beetles and wildfires, as described above.
NDVI Outside vs. Inside the Park¶
To compare how NDVI is changing outside vs. inside the park, I will be quantifying NDVI in a rectangle around the park boundary.
# Compute the area outside the RMNP
outside_rmnp_gdf = (
gpd.GeoDataFrame(geometry=rmnp_gdf.envelope)
.overlay(rmnp_gdf, how='difference')
)
# Taking the 'difference' allows us to select the area outside the boundary
# Plot the outside RMNP gdf
outside_rmnp_gdf.plot()
<Axes: >
Figure 6: This map shows the footprint of the park boundary, except it is inverted to select the area in a rectangle outside the park boundary. The rectangle is defined by the North/South/East/West extents of the park boundary.
# Clip NDVI data to both inside and outside the boundary
ndvi_inside = ndvi_da.rio.clip(rmnp_gdf.geometry, from_disk=True)
print(ndvi_inside)
ndvi_outside = ndvi_da.rio.clip(outside_rmnp_gdf.geometry, from_disk=True)
print(ndvi_outside)
<xarray.Dataset> Size: 11MB
Dimensions: (x: 201, y: 190, date: 72)
Coordinates:
band int64 8B 1
* x (x) float64 2kB -105.9 -105.9 -105.9 ... -105.5 -105.5 -105.5
* y (y) float64 2kB 40.55 40.55 40.55 40.55 ... 40.16 40.16 40.16
* date (date) datetime64[ns] 576B 2000-12-18 2001-01-01 ... 2024-01-17
spatial_ref int64 8B 0
Data variables:
NDVI (date, y, x) float32 11MB nan nan nan nan ... nan nan nan nan
<xarray.Dataset> Size: 11MB
Dimensions: (x: 202, y: 190, date: 72)
Coordinates:
band int64 8B 1
* x (x) float64 2kB -105.9 -105.9 -105.9 ... -105.5 -105.5 -105.5
* y (y) float64 2kB 40.55 40.55 40.55 40.55 ... 40.16 40.16 40.16
* date (date) datetime64[ns] 576B 2000-12-18 2001-01-01 ... 2024-01-17
spatial_ref int64 8B 0
Data variables:
NDVI (date, y, x) float32 11MB 0.1285 0.1285 0.0432 ... 0.341 0.4484
# Compute mean annual January NDVI
jan_ndvi_inside_df = (
ndvi_inside
.groupby(ndvi_inside.date.dt.year)
.mean(...)
.NDVI
.to_dataframe()
)
print("Mean NDVI inside RMNP Boundary", jan_ndvi_inside_df.head())
jan_ndvi_outside_df = (
ndvi_outside
.groupby(ndvi_outside.date.dt.year)
.mean(...)
.NDVI
.to_dataframe()
)
print("Mean NDVI outside RMNP Boundary", jan_ndvi_outside_df.head())
# Produces a NDVI Column that represents the average January NDVI for each year
Mean NDVI inside RMNP Boundary band spatial_ref NDVI year 2000 1 0 0.202534 2001 1 0 0.262729 2002 1 0 0.245032 2003 1 0 0.242526 2004 1 0 0.239577 Mean NDVI outside RMNP Boundary band spatial_ref NDVI year 2000 1 0 0.276744 2001 1 0 0.334754 2002 1 0 0.314590 2003 1 0 0.325537 2004 1 0 0.306118
# Join inside and outside RMNP df and rename columns
jan_ndvi_df = (
jan_ndvi_inside_df[['NDVI']]
.join(jan_ndvi_outside_df[['NDVI']],
lsuffix='_Inside_RMNP',
rsuffix='_Outside_RMNP')
)
jan_ndvi_df
| NDVI_Inside_RMNP | NDVI_Outside_RMNP | |
|---|---|---|
| year | ||
| 2000 | 0.202534 | 0.276744 |
| 2001 | 0.262729 | 0.334754 |
| 2002 | 0.245032 | 0.314590 |
| 2003 | 0.242526 | 0.325537 |
| 2004 | 0.239577 | 0.306118 |
| 2005 | 0.224774 | 0.305300 |
| 2006 | 0.203688 | 0.262113 |
| 2007 | 0.171129 | 0.214271 |
| 2008 | 0.187286 | 0.243634 |
| 2009 | 0.191575 | 0.255900 |
| 2010 | 0.172240 | 0.224974 |
| 2011 | 0.155885 | 0.217404 |
| 2012 | 0.195166 | 0.246179 |
| 2013 | 0.201369 | 0.259090 |
| 2014 | 0.164304 | 0.229065 |
| 2015 | 0.159257 | 0.214458 |
| 2016 | 0.166110 | 0.222576 |
| 2017 | 0.164789 | 0.219892 |
| 2018 | 0.205836 | 0.284953 |
| 2019 | 0.175625 | 0.243770 |
| 2020 | 0.174723 | 0.240372 |
| 2021 | 0.194546 | 0.244862 |
| 2022 | 0.131220 | 0.153884 |
| 2023 | 0.142583 | 0.177607 |
| 2024 | 0.195786 | 0.239010 |
# Plot mean NDVI inside and outside RMNP boundary
jan_ndvi_df.hvplot(title='Mean January NDVI Inside and Outside RMNP \n2000 - 2024')
Figure 7: This plot shows the mean January NDVI values in RMNP from 2000 to 2024. The plot shows that NDVI values outside the park within the rectangle boundary are higher than within the park, but annual trends for increases/decreases track very closely.
# Plot difference inside and outside the boundary
jan_ndvi_df['Difference'] = (
jan_ndvi_df['NDVI_Inside_RMNP']
- jan_ndvi_df['NDVI_Outside_RMNP']
)
jan_ndvi_df
| NDVI_Inside_RMNP | NDVI_Outside_RMNP | Difference | |
|---|---|---|---|
| year | |||
| 2000 | 0.202534 | 0.276744 | -0.074210 |
| 2001 | 0.262729 | 0.334754 | -0.072026 |
| 2002 | 0.245032 | 0.314590 | -0.069557 |
| 2003 | 0.242526 | 0.325537 | -0.083011 |
| 2004 | 0.239577 | 0.306118 | -0.066541 |
| 2005 | 0.224774 | 0.305300 | -0.080526 |
| 2006 | 0.203688 | 0.262113 | -0.058425 |
| 2007 | 0.171129 | 0.214271 | -0.043143 |
| 2008 | 0.187286 | 0.243634 | -0.056348 |
| 2009 | 0.191575 | 0.255900 | -0.064325 |
| 2010 | 0.172240 | 0.224974 | -0.052734 |
| 2011 | 0.155885 | 0.217404 | -0.061520 |
| 2012 | 0.195166 | 0.246179 | -0.051013 |
| 2013 | 0.201369 | 0.259090 | -0.057721 |
| 2014 | 0.164304 | 0.229065 | -0.064761 |
| 2015 | 0.159257 | 0.214458 | -0.055201 |
| 2016 | 0.166110 | 0.222576 | -0.056466 |
| 2017 | 0.164789 | 0.219892 | -0.055102 |
| 2018 | 0.205836 | 0.284953 | -0.079117 |
| 2019 | 0.175625 | 0.243770 | -0.068145 |
| 2020 | 0.174723 | 0.240372 | -0.065649 |
| 2021 | 0.194546 | 0.244862 | -0.050315 |
| 2022 | 0.131220 | 0.153884 | -0.022663 |
| 2023 | 0.142583 | 0.177607 | -0.035023 |
| 2024 | 0.195786 | 0.239010 | -0.043224 |
jan_ndvi_df.Difference.hvplot(
title='Difference in NDVI within and outside RMNP \n 2000 to 2020'
)
Figure 8: This plot represents the difference between NDVI values within and outside the RMNP boundary over the time period. Values closer to 0 indicate more similarity in NDVI values. Over time, the values are increasing towards 0, meaning the two different zones are becoming more similar.
Adding a Trend Line to NDVI Data Inside RMNP¶
To quantify the decrease in NDVI values that is observable in Figures 5 and 7, I will be adding a trendline (linear regression) to the data for NDVI inside the park. I will be using the linear regression method I learned in an earlier EarthLab assignment.
# Define the X and y values for the model
# I needed to reshape the index (year) to be 2-dimensional
X = jan_ndvi_df.index.values.reshape(-1, 1)
y = jan_ndvi_df["NDVI_Inside_RMNP"].values
# Method for running the linear regression using sklearn
model = LinearRegression()
model.fit(X, y)
# finding the slope with .coef_
slope = model.coef_[0]
# finding the intercept with .intercept
intercept = model.intercept_
# Create evenly spaced x-range for trendline
x_vals = np.linspace(
jan_ndvi_df.index.min(),
jan_ndvi_df.index.max(),
200)
# Applying the y=mx+b formula for lines
y_trend = slope * x_vals + intercept
#
fig, axes = plt.subplots(figsize=(10, 6))
axes.scatter(jan_ndvi_df.index, y, label="NDVI values")
axes.plot(x_vals, y_trend, linewidth=2, label="Trend line", color="red")
axes.set(
xlabel="Year",
ylabel="NDVI Inside RMNP",
title="NDVI Inside RMNP from 2000 to 2024",
)
axes.legend()
plt.show()
print("Slope:",slope)
Slope: -0.003081726351609596
Figure 9: Plot showing NDVI values from 2000 to 2024 inside RMNP with a linear regression trendline. The slope of the trendline is -0.00308, and has a clear negative trend on the plot, meaning that NDVI values have decreased in the park over the study period.
Conclusion¶
The results of Figures 5, 7, and 9 all indicate a significant decrease in NDVI values in Rocky Mountain National Park over the 24 year study period. In my data preparation and methods, I was careful to download data that would give results for pine trees specifically, by downloading data from January when other vegetation is dormant. This should remove potential errors that could have otherwise occurred if I used data from the summer for example.
The decrease in NDVI values can be explained by the pine beetle infestations starting in the early 2000s through 2009 (Mountain Pine Beetle. NPS). Additionally, significant wildfires in the park during 2012 and 2020 are another factor that is likely causing the decrease in NDVI values (Fire History. NPS)
This analysis shows the usefulness of using NDVI data to quantify changes in vegetation health. This analysis could be taken further by using other datasets, such as wildfire maps, or beetle kill maps (if they exist) to cross check the results found above. Additionally, knowing vegetation health trends is an important tool for conservation efforts and climate change understanding and mitigation.
Works Cited:¶
Colorado State Forest Service. Mountain Pine Beetle. https://csfs.colostate.edu/forest-management/common-forest-insects-diseases/mountain-pine-beetle/
Colorado State Forest Service. Forests and Trees. https://csfs.colostate.edu/forests-trees/
NASA EarthData. Normalized Difference Vegetation Index (NDVI). https://www.earthdata.nasa.gov/topics/land-surface/normalized-difference-vegetation-index-ndvi
National Park Service. Mountain Pine Beetle. Rocky Mountain National Park. https://www.nps.gov/romo/learn/nature/mtn_pine_beetle_background.htm
National Park Service. Fire History. Rocky Mountain National Park. https://www.nps.gov/romo/learn/fire-history.htm