The Cameron Peak Fire, Colorado, USA¶

This analysis looks at Cameron Peak and other fires

No description has been provided for this image

It’s another ESIIL Earth Data Science Workflow¶

This notebook contains your next environmental data science (EDS) coding challenge! Before we get started, make sure to read or review the guidelines below. These will help make sure that your code is readable and reproducible.

Don’t get caught by these interactive coding notebook gotchas¶

No description has been provided for this image

Image source: https://alaskausfws.medium.com/whats-big-and-brown-and-loves-salmon-e1803579ee36

These are the most common issues that will keep you from getting started and delay your code review:

Run your code in the right environment to avoid import errors¶

We’ve created a coding environment for you to use that already has all the software and libraries you will need! When you try to run some code, you may be prompted to select a kernel. The kernel refers to the version of Python you are using. You should use the base kernel, which should be the default option.

Always run your code start to finish before submitting¶

Before you commit your work, make sure it runs reproducibly by clicking:

  1. Restart (this button won’t appear until you’ve run some code), then
  2. Run All

Check your code to make sure it’s clean and easy to read¶

No description has been provided for this image

  • Format all cells prior to submitting (right click on your code).

  • Use expressive names for variables so you or the reader knows what they are.

  • Use comments to explain your code – e.g.

    # This is a comment, it starts with a hash sign

Label and describe your plots¶

Source: https://xkcd.com/833
Source: https://xkcd.com/833

Make sure each plot has:

  • A title that explains where and when the data are from
  • x- and y- axis labels with units where appropriate
  • A legend where appropriate

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.

Healthy leaves reflect a lot of NIR radiation compared to dead or stressed leaves > 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.

In [ ]:
# Install the development version of the earthpy package
!pip install git+https://github.com/earthlab/earthpy@apppears
In [1]:
# Import python libraries
import getpass
import json
import os
import pathlib
from glob import glob

# Library to work with tabular data
import pandas as pd

# Library to work with 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 # uses numpy
import matplotlib.pyplot as plt
In [3]:
fire_path = "/workspaces/shortcourse-02-vegetation-nquarder/data/fire_archive_M-C61_464561.shp"
fire_gdf = gpd.read_file(fire_path)
fire_gdf
Out[3]:
LATITUDE LONGITUDE BRIGHTNESS SCAN TRACK ACQ_DATE ACQ_TIME SATELLITE INSTRUMENT CONFIDENCE VERSION BRIGHT_T31 FRP DAYNIGHT TYPE geometry
0 46.7993 -113.3773 315.1 1.4 1.2 2005-07-21 2052 Aqua MODIS 30 6.03 304.7 8.2 D 0 POINT (-113.37730 46.79930)
1 46.7975 -113.3961 320.3 1.4 1.2 2005-07-21 2052 Aqua MODIS 65 6.03 308.4 16.3 D 0 POINT (-113.39610 46.79750)
2 46.8114 -111.8638 323.4 1.7 1.3 2005-07-23 1859 Terra MODIS 76 6.03 304.9 26.8 D 0 POINT (-111.86380 46.81140)
3 46.8159 -111.8784 325.4 1.2 1.1 2005-07-23 2040 Aqua MODIS 62 6.03 312.9 13.0 D 0 POINT (-111.87840 46.81590)
4 49.1815 -113.2053 313.2 1.1 1.0 2005-07-27 2016 Aqua MODIS 38 6.03 299.4 5.4 D 0 POINT (-113.20530 49.18150)
... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ...
23296 47.3181 -114.4514 324.0 1.0 1.0 2020-07-30 0538 Terra MODIS 100 6.03 289.8 21.5 N 0 POINT (-114.45140 47.31810)
23297 47.3163 -114.4651 337.7 1.0 1.0 2020-07-30 0538 Terra MODIS 100 6.03 293.6 40.3 N 0 POINT (-114.46510 47.31630)
23298 47.3146 -114.4786 313.3 1.0 1.0 2020-07-30 0538 Terra MODIS 87 6.03 290.4 12.0 N 0 POINT (-114.47860 47.31460)
23299 47.3171 -114.4624 310.9 1.0 1.0 2020-07-30 0947 Aqua MODIS 44 6.03 293.5 9.7 N 0 POINT (-114.46240 47.31710)
23300 47.3152 -114.4494 316.6 1.0 1.0 2020-07-30 0947 Aqua MODIS 94 6.03 295.0 14.2 N 0 POINT (-114.44940 47.31520)

23301 rows × 16 columns

In [4]:
fire_gdf.explore()
Out[4]:
Make this Notebook Trusted to load map: File -> Trust Notebook

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.

In [6]:
data_dir = os.path.join(pathlib.Path.home(), 'cameron-peak-data')
# Make the data directory
os.makedirs(data_dir, exist_ok=True)

data_dir
Out[6]:
'/home/jovyan/cameron-peak-data'

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.
In [5]:
# Download the Cameron Peak fire boundary
cameron_peak_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")
cameron_peak_url

# Open data with geopandas
cameron_peak_gdf = gpd.read_file(cameron_peak_url)
cameron_peak_gdf
Out[5]:
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 NaN 1678806822810 1602019800000 {53741A13-D269-4CD5-AF91-02E094B944DA} ... 1638820092647 FODR NaN NaN NaN {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

In [10]:
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
In [4]:
# Plot the Cameron Peak Fire boundary
cameron_peak_gdf.explore()
Out[4]:
Make this Notebook Trusted to load map: File -> Trust Notebook

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?

In [12]:
# 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=cameron_peak_gdf
)

ndvi_downloader.download_files(cache=True)

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.

In [8]:
# Get a list of NDVI tif file paths
ndvi_paths = sorted(glob(os.path.join(data_dir, 'cp-ndvi', '*', '*NDVI*.tif')))
len(ndvi_paths)
Out[8]:
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
In [3]:
scale_factor = 10000
doy_start = -19
doy_end = -12
In [10]:
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[10]:
18
In [11]:
ndvi_das
Out[11]:
[<xarray.DataArray 'NDVI' (date: 1, y: 153, x: 339)>
 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]]],
       dtype=float32)
 Coordinates:
     band         int64 1
   * x            (x) float64 -105.9 -105.9 -105.9 ... -105.2 -105.2 -105.2
   * y            (y) float64 40.78 40.78 40.77 40.77 ... 40.47 40.47 40.46 40.46
     spatial_ref  int64 0
   * date         (date) datetime64[ns] 2018-06-26,
 <xarray.DataArray 'NDVI' (date: 1, y: 153, x: 339)>
 array([[[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]]],
       dtype=float32)
 Coordinates:
     band         int64 1
   * x            (x) float64 -105.9 -105.9 -105.9 ... -105.2 -105.2 -105.2
   * y            (y) float64 40.78 40.78 40.77 40.77 ... 40.47 40.47 40.46 40.46
     spatial_ref  int64 0
   * date         (date) datetime64[ns] 2018-07-12,
 <xarray.DataArray 'NDVI' (date: 1, y: 153, x: 339)>
 array([[[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.2342, 0.2342, 0.2641, ..., 0.6167, 0.5051, 0.4424],
         [0.3243, 0.3079, 0.2732, ..., 0.5653, 0.5653, 0.5468],
         [0.4528, 0.2801, 0.2443, ..., 0.5428, 0.5418, 0.5418]]],
       dtype=float32)
 Coordinates:
     band         int64 1
   * x            (x) float64 -105.9 -105.9 -105.9 ... -105.2 -105.2 -105.2
   * y            (y) float64 40.78 40.78 40.77 40.77 ... 40.47 40.47 40.46 40.46
     spatial_ref  int64 0
   * date         (date) datetime64[ns] 2018-07-28,
 <xarray.DataArray 'NDVI' (date: 1, y: 153, x: 339)>
 array([[[0.6198, 0.6046, 0.6046, ..., 0.6461, 0.6461, 0.5353],
         [0.6056, 0.6964, 0.6964, ..., 0.5134, 0.6104, 0.5998],
         [0.6625, 0.6625, 0.6964, ..., 0.5429, 0.5478, 0.5478],
         ...,
         [0.0373, 0.0373, 0.0338, ..., 0.5641, 0.5623, 0.5623],
         [0.0735, 0.0511, 0.0428, ..., 0.5022, 0.5729, 0.5217],
         [0.0958, 0.0572, 0.0457, ..., 0.602 , 0.5729, 0.5729]]],
       dtype=float32)
 Coordinates:
     band         int64 1
   * x            (x) float64 -105.9 -105.9 -105.9 ... -105.2 -105.2 -105.2
   * y            (y) float64 40.78 40.78 40.77 40.77 ... 40.47 40.47 40.46 40.46
     spatial_ref  int64 0
   * date         (date) datetime64[ns] 2019-06-26,
 <xarray.DataArray 'NDVI' (date: 1, y: 153, x: 339)>
 array([[[0.6042, 0.6158, 0.5882, ..., 0.4353, 0.4353, 0.4208],
         [0.6042, 0.671 , 0.671 , ..., 0.4236, 0.4202, 0.4208],
         [0.6147, 0.6147, 0.667 , ..., 0.4149, 0.4087, 0.4694],
         ...,
         [0.1867, 0.1867, 0.1722, ..., 0.4773, 0.4938, 0.4302],
         [0.2452, 0.2452, 0.1384, ..., 0.4151, 0.5035, 0.4979],
         [0.3021, 0.2488, 0.1415, ..., 0.4289, 0.5248, 0.5248]]],
       dtype=float32)
 Coordinates:
     band         int64 1
   * x            (x) float64 -105.9 -105.9 -105.9 ... -105.2 -105.2 -105.2
   * y            (y) float64 40.78 40.78 40.77 40.77 ... 40.47 40.47 40.46 40.46
     spatial_ref  int64 0
   * date         (date) datetime64[ns] 2019-07-12,
 <xarray.DataArray 'NDVI' (date: 1, y: 153, x: 339)>
 array([[[0.6287, 0.6287, 0.6362, ..., 0.4046, 0.4046, 0.4021],
         [0.6292, 0.697 , 0.697 , ..., 0.4679, 0.3959, 0.4021],
         [0.625 , 0.625 , 0.6793, ..., 0.4253, 0.3919, 0.4128],
         ...,
         [0.2745, 0.2745, 0.1563, ..., 0.4104, 0.4541, 0.3733],
         [0.3348, 0.3348, 0.2914, ..., 0.4412, 0.4291, 0.4518],
         [0.436 , 0.347 , 0.2914, ..., 0.387 , 0.4627, 0.4627]]],
       dtype=float32)
 Coordinates:
     band         int64 1
   * x            (x) float64 -105.9 -105.9 -105.9 ... -105.2 -105.2 -105.2
   * y            (y) float64 40.78 40.78 40.77 40.77 ... 40.47 40.47 40.46 40.46
     spatial_ref  int64 0
   * date         (date) datetime64[ns] 2019-07-28,
 <xarray.DataArray 'NDVI' (date: 1, y: 153, x: 339)>
 array([[[0.5801, 0.589 , 0.5523, ..., 0.3445, 0.3445, 0.3624],
         [0.5898, 0.6868, 0.6868, ..., 0.4048, 0.3638, 0.3586],
         [0.6143, 0.6143, 0.6868, ..., 0.3604, 0.3733, 0.3814],
         ...,
         [0.2375, 0.2375, 0.2375, ..., 0.4297, 0.4803, 0.3789],
         [0.3265, 0.2375, 0.2375, ..., 0.3969, 0.4803, 0.443 ],
         [0.3419, 0.2943, 0.2736, ..., 0.4618, 0.4788, 0.4788]]],
       dtype=float32)
 Coordinates:
     band         int64 1
   * x            (x) float64 -105.9 -105.9 -105.9 ... -105.2 -105.2 -105.2
   * y            (y) float64 40.78 40.78 40.77 40.77 ... 40.47 40.47 40.46 40.46
     spatial_ref  int64 0
   * date         (date) datetime64[ns] 2020-06-25,
 <xarray.DataArray 'NDVI' (date: 1, y: 153, x: 339)>
 array([[[0.5915, 0.6164, 0.5988, ..., 0.3275, 0.3275, 0.301 ],
         [0.6295, 0.6295, 0.6295, ..., 0.3328, 0.3275, 0.3192],
         [0.6354, 0.6354, 0.6781, ..., 0.3259, 0.3328, 0.3265],
         ...,
         [0.2495, 0.2495, 0.2156, ..., 0.4628, 0.4628, 0.367 ],
         [0.311 , 0.297 , 0.297 , ..., 0.4088, 0.4088, 0.3885],
         [0.4352, 0.3389, 0.297 , ..., 0.4088, 0.4088, 0.4088]]],
       dtype=float32)
 Coordinates:
     band         int64 1
   * x            (x) float64 -105.9 -105.9 -105.9 ... -105.2 -105.2 -105.2
   * y            (y) float64 40.78 40.78 40.77 40.77 ... 40.47 40.47 40.46 40.46
     spatial_ref  int64 0
   * date         (date) datetime64[ns] 2020-07-11,
 <xarray.DataArray 'NDVI' (date: 1, y: 153, x: 339)>
 array([[[0.6042, 0.6054, 0.5813, ..., 0.3002, 0.3002, 0.3313],
         [0.6042, 0.637 , 0.637 , ..., 0.3337, 0.3161, 0.3058],
         [0.6605, 0.6605, 0.6756, ..., 0.3207, 0.3207, 0.3155],
         ...,
         [0.2475, 0.2475, 0.2296, ..., 0.4566, 0.4367, 0.318 ],
         [0.3478, 0.2492, 0.2296, ..., 0.3201, 0.4268, 0.3801],
         [0.3973, 0.2501, 0.2501, ..., 0.4786, 0.4645, 0.4645]]],
       dtype=float32)
 Coordinates:
     band         int64 1
   * x            (x) float64 -105.9 -105.9 -105.9 ... -105.2 -105.2 -105.2
   * y            (y) float64 40.78 40.78 40.77 40.77 ... 40.47 40.47 40.46 40.46
     spatial_ref  int64 0
   * date         (date) datetime64[ns] 2020-07-27,
 <xarray.DataArray 'NDVI' (date: 1, y: 153, x: 339)>
 array([[[0.6083, 0.6195, 0.6195, ..., 0.4828, 0.4828, 0.4954],
         [0.6012, 0.6837, 0.6837, ..., 0.4801, 0.4672, 0.4688],
         [0.6354, 0.6354, 0.6837, ..., 0.4705, 0.4701, 0.5122],
         ...,
         [0.2752, 0.2752, 0.2015, ..., 0.6013, 0.6057, 0.5   ],
         [0.2606, 0.2918, 0.2739, ..., 0.5711, 0.5975, 0.5975],
         [0.4402, 0.3109, 0.2506, ..., 0.5632, 0.6123, 0.6123]]],
       dtype=float32)
 Coordinates:
     band         int64 1
   * x            (x) float64 -105.9 -105.9 -105.9 ... -105.2 -105.2 -105.2
   * y            (y) float64 40.78 40.78 40.77 40.77 ... 40.47 40.47 40.46 40.46
     spatial_ref  int64 0
   * date         (date) datetime64[ns] 2021-06-26,
 <xarray.DataArray 'NDVI' (date: 1, y: 153, x: 339)>
 array([[[0.6103, 0.6302, 0.6302, ..., 0.3999, 0.3999, 0.4207],
         [0.6277, 0.6677, 0.6677, ..., 0.447 , 0.391 , 0.4075],
         [0.6643, 0.6643, 0.6758, ..., 0.3848, 0.4462, 0.391 ],
         ...,
         [0.2796, 0.2796, 0.226 , ..., 0.5718, 0.587 , 0.475 ],
         [0.2692, 0.2744, 0.2602, ..., 0.5335, 0.5589, 0.475 ],
         [0.3598, 0.2473, 0.2005, ..., 0.5754, 0.5624, 0.5624]]],
       dtype=float32)
 Coordinates:
     band         int64 1
   * x            (x) float64 -105.9 -105.9 -105.9 ... -105.2 -105.2 -105.2
   * y            (y) float64 40.78 40.78 40.77 40.77 ... 40.47 40.47 40.46 40.46
     spatial_ref  int64 0
   * date         (date) datetime64[ns] 2021-07-12,
 <xarray.DataArray 'NDVI' (date: 1, y: 153, x: 339)>
 array([[[0.6004, 0.6111, 0.577 , ..., 0.4086, 0.4086, 0.4266],
         [0.5732, 0.6728, 0.6728, ..., 0.4042, 0.4266, 0.4266],
         [0.5732, 0.5732, 0.6728, ..., 0.4042, 0.4217, 0.4217],
         ...,
         [0.2494, 0.2494, 0.2494, ..., 0.5226, 0.4771, 0.4771],
         [0.3017, 0.2898, 0.171 , ..., 0.5065, 0.4771, 0.4771],
         [0.373 , 0.373 , 0.226 , ..., 0.5248, 0.4838, 0.4838]]],
       dtype=float32)
 Coordinates:
     band         int64 1
   * x            (x) float64 -105.9 -105.9 -105.9 ... -105.2 -105.2 -105.2
   * y            (y) float64 40.78 40.78 40.77 40.77 ... 40.47 40.47 40.46 40.46
     spatial_ref  int64 0
   * date         (date) datetime64[ns] 2021-07-28,
 <xarray.DataArray 'NDVI' (date: 1, y: 153, x: 339)>
 array([[[0.5762, 0.581 , 0.5797, ..., 0.2721, 0.2721, 0.2937],
         [0.5721, 0.6287, 0.6287, ..., 0.3098, 0.3098, 0.2807],
         [0.6161, 0.6161, 0.6561, ..., 0.2896, 0.2894, 0.3165],
         ...,
         [0.2301, 0.2301, 0.1457, ..., 0.5113, 0.4463, 0.3555],
         [0.559 , 0.559 , 0.1457, ..., 0.3438, 0.4463, 0.3555],
         [0.3742, 0.2145, 0.2145, ..., 0.4332, 0.3889, 0.3889]]],
       dtype=float32)
 Coordinates:
     band         int64 1
   * x            (x) float64 -105.9 -105.9 -105.9 ... -105.2 -105.2 -105.2
   * y            (y) float64 40.78 40.78 40.77 40.77 ... 40.47 40.47 40.46 40.46
     spatial_ref  int64 0
   * date         (date) datetime64[ns] 2022-06-26,
 <xarray.DataArray 'NDVI' (date: 1, y: 153, x: 339)>
 array([[[0.6245, 0.664 , 0.6318, ..., 0.297 , 0.297 , 0.3184],
         [0.6245, 0.6975, 0.6975, ..., 0.3189, 0.3189, 0.3027],
         [0.6408, 0.6408, 0.6751, ..., 0.3086, 0.3373, 0.3027],
         ...,
         [0.2819, 0.2819, 0.103 , ..., 0.5089, 0.4961, 0.3709],
         [0.3659, 0.3659, 0.2201, ..., 0.3797, 0.4028, 0.4028],
         [0.4405, 0.3574, 0.2201, ..., 0.43  , 0.4036, 0.4036]]],
       dtype=float32)
 Coordinates:
     band         int64 1
   * x            (x) float64 -105.9 -105.9 -105.9 ... -105.2 -105.2 -105.2
   * y            (y) float64 40.78 40.78 40.77 40.77 ... 40.47 40.47 40.46 40.46
     spatial_ref  int64 0
   * date         (date) datetime64[ns] 2022-07-12,
 <xarray.DataArray 'NDVI' (date: 1, y: 153, x: 339)>
 array([[[0.5689, 0.5789, 0.5478, ..., 0.3454, 0.3454, 0.3319],
         [0.5781, 0.5912, 0.5912, ..., 0.3536, 0.3364, 0.3176],
         [0.6664, 0.6664, 0.6961, ..., 0.3359, 0.326 , 0.326 ],
         ...,
         [0.3177, 0.3177, 0.3177, ..., 0.5063, 0.5054, 0.349 ],
         [0.4454, 0.2504, 0.2833, ..., 0.4733, 0.4538, 0.4029],
         [0.4454, 0.275 , 0.262 , ..., 0.5324, 0.4879, 0.4879]]],
       dtype=float32)
 Coordinates:
     band         int64 1
   * x            (x) float64 -105.9 -105.9 -105.9 ... -105.2 -105.2 -105.2
   * y            (y) float64 40.78 40.78 40.77 40.77 ... 40.47 40.47 40.46 40.46
     spatial_ref  int64 0
   * date         (date) datetime64[ns] 2022-07-28,
 <xarray.DataArray 'NDVI' (date: 1, y: 153, x: 339)>
 array([[[0.5411, 0.5622, 0.5526, ..., 0.4669, 0.4669, 0.5193],
         [0.5322, 0.6519, 0.6519, ..., 0.4944, 0.4669, 0.5075],
         [0.566 , 0.566 , 0.6079, ..., 0.4613, 0.4944, 0.4638],
         ...,
         [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]]],
       dtype=float32)
 Coordinates:
     band         int64 1
   * x            (x) float64 -105.9 -105.9 -105.9 ... -105.2 -105.2 -105.2
   * y            (y) float64 40.78 40.78 40.77 40.77 ... 40.47 40.47 40.46 40.46
     spatial_ref  int64 0
   * date         (date) datetime64[ns] 2023-06-26,
 <xarray.DataArray 'NDVI' (date: 1, y: 153, x: 339)>
 array([[[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]]],
       dtype=float32)
 Coordinates:
     band         int64 1
   * x            (x) float64 -105.9 -105.9 -105.9 ... -105.2 -105.2 -105.2
   * y            (y) float64 40.78 40.78 40.77 40.77 ... 40.47 40.47 40.46 40.46
     spatial_ref  int64 0
   * date         (date) datetime64[ns] 2023-07-12,
 <xarray.DataArray 'NDVI' (date: 1, y: 153, x: 339)>
 array([[[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)
 Coordinates:
     band         int64 1
   * x            (x) float64 -105.9 -105.9 -105.9 ... -105.2 -105.2 -105.2
   * y            (y) float64 40.78 40.78 40.77 40.77 ... 40.47 40.47 40.46 40.46
     spatial_ref  int64 0
   * date         (date) datetime64[ns] 2023-07-28]

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

In [12]:
ndvi_da = xr.combine_by_coords(ndvi_das, coords=['date'])
ndvi_da
Out[12]:
<xarray.Dataset>
Dimensions:      (x: 339, y: 153, date: 18)
Coordinates:
    band         int64 1
  * x            (x) float64 -105.9 -105.9 -105.9 ... -105.2 -105.2 -105.2
  * y            (y) float64 40.78 40.78 40.77 40.77 ... 40.47 40.47 40.46 40.46
    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.5672 0.5711 0.6027 ... 0.6327 0.6327
xarray.Dataset
    • x: 339
    • y: 153
    • date: 18
    • 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]')
    • 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)
    • 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))
In [14]:
# 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
)
In [16]:
# Plot the difference
fig, ax = plt.subplots()

ndvi_diff.plot(x='x', y='y', cmap='PiYG', robust=True, ax=ax)
cameron_peak_gdf.plot(edgecolor='black', color='none', ax=ax)
plt.title("NDVI - After vs. Before Cameron Peak Fire, 2020")
plt.show()
No description has been provided for this image

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 [17]:
out_gdf = (
    gpd.GeoDataFrame(geometry=cameron_peak_gdf.envelope)
    .overlay(cameron_peak_gdf, how='difference'))
out_gdf
Out[17]:
geometry
0 MULTIPOLYGON (((-105.90371 40.46134, -105.9037...

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

In [18]:
# Clip data to both inside and outside the boundary
ndvi_cp_da = ndvi_da.rio.clip(cameron_peak_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.

:::

In [19]:
# 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')
)

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 [20]:
# Plot difference inside and outside the reservation
jul_ndvi_df.plot(
    title='NDVI before and after the Cameron Peak Fire')
Out[20]:
<Axes: title={'center': 'NDVI before and after the Cameron Peak Fire'}, xlabel='year'>
No description has been provided for this image

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

Looking at Lagos, Nigeria¶

In [21]:
lagos_dir = os.path.join(pathlib.Path.home(), 'lagos-data')
# Make the data directory
os.makedirs(lagos_dir, exist_ok=True)

lagos_dir
Out[21]:
'/home/jovyan/lagos-data'
In [2]:
# Define url to Nigeria .shp
nigeria_url = "https://open.africa/dataset/83582021-d8f7-4bfd-928f-e9c6c8cb1247/resource/372a616a-66cc-41f7-ac91-d8af8f23bc2b/download/nigeria-lgas.zip"

# Open Nigeria .shp using geopandas
nigeria_gdf = gpd.read_file(nigeria_url)
nigeria_gdf
Out[2]:
STATE LGA AREA PERIMETER LONGITUDE LATITUDE FULL_NAME geometry
0 Sokoto Gada 1193.977 170.095 NaN NaN NaN POLYGON ((5.53632 13.88793, 5.53480 13.88488, ...
1 Sokoto Illela 1298.423 174.726 NaN NaN NaN POLYGON ((5.53632 13.88793, 5.54517 13.88419, ...
2 Sokoto Tangaza 2460.715 209.702 NaN NaN NaN POLYGON ((4.85548 13.76724, 4.86189 13.78085, ...
3 Borno Abadam 2430.515 288.957 NaN NaN NaN POLYGON ((12.83189 13.39871, 12.83397 13.40439...
4 Lake Lake chad 5225.912 497.039 NaN NaN NaN POLYGON ((13.48608 13.30821, 13.48296 13.31344...
... ... ... ... ... ... ... ... ...
770 Delta Isoko North 485.467 169.369 NaN NaN NaN MULTIPOLYGON (((6.31996 5.63341, 6.32003 5.633...
771 Niger Lavun 3951.431 424.153 NaN NaN NaN MULTIPOLYGON (((6.12188 9.09441, 6.12001 9.094...
772 Yobe Bade 817.260 216.207 NaN NaN NaN MULTIPOLYGON (((11.01052 12.80457, 11.00747 12...
773 Zamfara Maru 7795.261 536.500 NaN NaN NaN MULTIPOLYGON (((6.43894 12.41104, 6.43609 12.4...
774 Akwa Ibom Oron 81.472 57.846 NaN NaN NaN POLYGON ((8.22063 4.84473, 8.23405 4.82974, 8....

775 rows × 8 columns

In [23]:
# Select out Lagos
lagos = nigeria_gdf[nigeria_gdf["STATE"]=="Lagos"]
lagos
Out[23]:
STATE LGA AREA PERIMETER LONGITUDE LATITUDE FULL_NAME geometry
542 Lagos Ifako/Ijaye 33.615 28.115 NaN NaN NaN POLYGON ((3.26599 6.70696, 3.26787 6.70693, 3....
543 Lagos Alimosho 148.062 55.194 NaN NaN NaN POLYGON ((3.28961 6.67267, 3.28807 6.66979, 3....
546 Lagos Ikorodu 369.494 88.090 NaN NaN NaN POLYGON ((3.69866 6.68905, 3.69502 6.68524, 3....
549 Lagos Agege 16.575 18.782 NaN NaN NaN POLYGON ((3.33228 6.64879, 3.33108 6.64763, 3....
551 Lagos Ikeja 45.744 32.736 NaN NaN NaN POLYGON ((3.35382 6.66920, 3.35453 6.66877, 3....
554 Lagos Kosofe 74.118 38.369 NaN NaN NaN POLYGON ((3.43208 6.64850, 3.43255 6.64591, 3....
558 Lagos Lagos Island 238.990 94.966 NaN NaN NaN POLYGON ((3.46255 6.60946, 3.46357 6.60944, 3....
560 Lagos Oshodi/Isolo 38.595 28.349 NaN NaN NaN POLYGON ((3.33961 6.58872, 3.34026 6.58733, 3....
561 Lagos Shomolu 18.304 17.122 NaN NaN NaN POLYGON ((3.38974 6.57181, 3.39050 6.57035, 3....
562 Lagos Mushin 17.387 18.422 NaN NaN NaN POLYGON ((3.36006 6.57612, 3.36123 6.57177, 3....
564 Lagos Ibeju Lekki 488.075 157.064 NaN NaN NaN POLYGON ((3.61947 6.41924, 3.62283 6.42978, 3....
565 Lagos Ojo 196.583 68.195 NaN NaN NaN POLYGON ((3.21582 6.52357, 3.21650 6.52330, 3....
567 Lagos Lagos Mainland 22.781 22.449 NaN NaN NaN POLYGON ((3.39173 6.53537, 3.39223 6.53461, 3....
568 Lagos Surulere 30.717 23.004 NaN NaN NaN POLYGON ((3.35727 6.53463, 3.35851 6.53242, 3....
569 Lagos Amuwo Odofin 190.955 67.320 NaN NaN NaN POLYGON ((3.27456 6.52441, 3.27533 6.52369, 3....
570 Lagos Badagry 516.531 115.339 NaN NaN NaN POLYGON ((3.08298 6.52197, 3.08746 6.51811, 3....
574 Lagos Ajeromi/ Ifelodun 11.645 13.097 NaN NaN NaN POLYGON ((3.33895 6.47589, 3.33891 6.47583, 3....
749 Lagos Epe 1205.291 247.437 NaN NaN NaN POLYGON ((3.90320 6.68705, 4.06281 6.68550, 4....
764 Lagos Eti-Osa 201.279 85.069 NaN NaN NaN POLYGON ((3.66212 6.49440, 3.66465 6.49059, 3....
765 Lagos Apapa 43.227 34.006 NaN NaN NaN POLYGON ((3.38973 6.40838, 3.39593 6.39834, 3....
In [38]:
# Interactive map of Lagos
lagos.explore()
Out[38]:
Make this Notebook Trusted to load map: File -> Trust Notebook
In [24]:
# Select Alimosho
alimosho = lagos[lagos["LGA"]=='Alimosho']
alimosho
Out[24]:
STATE LGA AREA PERIMETER LONGITUDE LATITUDE FULL_NAME geometry
543 Lagos Alimosho 148.062 55.194 NaN NaN NaN POLYGON ((3.28961 6.67267, 3.28807 6.66979, 3....
In [40]:
alimosho.explore()
Out[40]:
Make this Notebook Trusted to load map: File -> Trust Notebook
In [42]:
# Initialize AppeearsDownloader for MODIS NDVI data
ndvi_downloader = eaapp.AppeearsDownloader(
    download_key='lagos-ndvi',
    ea_dir=lagos_dir,
    product='MOD13Q1.061',
    layer='_250m_16_days_NDVI',
    start_date="07-01",
    end_date="07-31",
    recurring=True,
    year_range=[2015, 2020],
    polygon=alimosho
)

ndvi_downloader.download_files(cache=True)
In [25]:
# Get a list of NDVI tif file paths
lagos_paths = sorted(glob(os.path.join(lagos_dir, 'lagos-ndvi', '*', '*NDVI*.tif')))
len(lagos_paths)
Out[25]:
18
In [26]:
lagos_das = []
for lagos_path in lagos_paths:
    # Get date from file name
    doy = lagos_path[doy_start:doy_end]
    date = pd.to_datetime(doy, format='%Y%j')

    # Open dataset
    da = rxr.open_rasterio(lagos_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
    lagos_das.append(da)

len(lagos_das)
Out[26]:
18
In [27]:
lagos_da = xr.combine_by_coords(lagos_das, coords=['date'])
lagos_da
Out[27]:
<xarray.Dataset>
Dimensions:      (x: 56, y: 89, date: 18)
Coordinates:
    band         int64 1
  * x            (x) float64 3.186 3.189 3.191 3.193 ... 3.295 3.297 3.299 3.301
  * y            (y) float64 6.705 6.703 6.701 6.699 ... 6.528 6.526 6.524 6.522
    spatial_ref  int64 0
  * date         (date) datetime64[ns] 2015-06-26 2015-07-12 ... 2020-07-27
Data variables:
    NDVI         (date, y, x) float32 0.6021 0.6098 0.5921 ... 0.2961 0.297
xarray.Dataset
    • x: 56
    • y: 89
    • date: 18
    • band
      ()
      int64
      1
      array(1)
    • x
      (x)
      float64
      3.186 3.189 3.191 ... 3.299 3.301
      array([3.186458, 3.188542, 3.190625, 3.192708, 3.194792, 3.196875, 3.198958,
       3.201042, 3.203125, 3.205208, 3.207292, 3.209375, 3.211458, 3.213542,
       3.215625, 3.217708, 3.219792, 3.221875, 3.223958, 3.226042, 3.228125,
       3.230208, 3.232292, 3.234375, 3.236458, 3.238542, 3.240625, 3.242708,
       3.244792, 3.246875, 3.248958, 3.251042, 3.253125, 3.255208, 3.257292,
       3.259375, 3.261458, 3.263542, 3.265625, 3.267708, 3.269792, 3.271875,
       3.273958, 3.276042, 3.278125, 3.280208, 3.282292, 3.284375, 3.286458,
       3.288542, 3.290625, 3.292708, 3.294792, 3.296875, 3.298958, 3.301042])
    • y
      (y)
      float64
      6.705 6.703 6.701 ... 6.524 6.522
      array([6.705208, 6.703125, 6.701042, 6.698958, 6.696875, 6.694792, 6.692708,
       6.690625, 6.688542, 6.686458, 6.684375, 6.682292, 6.680208, 6.678125,
       6.676042, 6.673958, 6.671875, 6.669792, 6.667708, 6.665625, 6.663542,
       6.661458, 6.659375, 6.657292, 6.655208, 6.653125, 6.651042, 6.648958,
       6.646875, 6.644792, 6.642708, 6.640625, 6.638542, 6.636458, 6.634375,
       6.632292, 6.630208, 6.628125, 6.626042, 6.623958, 6.621875, 6.619792,
       6.617708, 6.615625, 6.613542, 6.611458, 6.609375, 6.607292, 6.605208,
       6.603125, 6.601042, 6.598958, 6.596875, 6.594792, 6.592708, 6.590625,
       6.588542, 6.586458, 6.584375, 6.582292, 6.580208, 6.578125, 6.576042,
       6.573958, 6.571875, 6.569792, 6.567708, 6.565625, 6.563542, 6.561458,
       6.559375, 6.557292, 6.555208, 6.553125, 6.551042, 6.548958, 6.546875,
       6.544792, 6.542708, 6.540625, 6.538542, 6.536458, 6.534375, 6.532292,
       6.530208, 6.528125, 6.526042, 6.523958, 6.521875])
    • 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 :
      3.1854166663813004 0.0020833333331466974 0.0 6.706249999399219 0.0 -0.0020833333331466974
      array(0)
    • date
      (date)
      datetime64[ns]
      2015-06-26 ... 2020-07-27
      array(['2015-06-26T00:00:00.000000000', '2015-07-12T00:00:00.000000000',
       '2015-07-28T00:00:00.000000000', '2016-06-25T00:00:00.000000000',
       '2016-07-11T00:00:00.000000000', '2016-07-27T00:00:00.000000000',
       '2017-06-26T00:00:00.000000000', '2017-07-12T00:00:00.000000000',
       '2017-07-28T00:00:00.000000000', '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'],
      dtype='datetime64[ns]')
    • NDVI
      (date, y, x)
      float32
      0.6021 0.6098 ... 0.2961 0.297
      array([[[0.6021, 0.6098, 0.5921, ..., 0.1796, 0.198 , 0.202 ],
        [0.6021, 0.6021, 0.4927, ..., 0.0988, 0.0988, 0.0988],
        [0.483 , 0.3454, 0.3454, ..., 0.0645, 0.0725, 0.0725],
        ...,
        [0.3179, 0.4178, 0.4178, ..., 0.2027, 0.1739, 0.1728],
        [0.3179, 0.3079, 0.4311, ..., 0.4688, 0.1997, 0.1997],
        [0.2625, 0.3079, 0.3191, ..., 0.633 , 0.3045, 0.1636]],

       [[0.5312, 0.5312, 0.6211, ..., 0.3192, 0.328 , 0.328 ],
        [0.5707, 0.5707, 0.4551, ..., 0.3192, 0.2925, 0.1879],
        [0.4171, 0.2903, 0.2903, ..., 0.332 , 0.2105, 0.1741],
        ...,
        [0.8989, 0.8989, 0.8994, ..., 0.3523, 0.2193, 0.0826],
        [0.592 , 0.4235, 0.3674, ..., 0.2957, 0.2206, 0.0938],
        [0.592 , 0.4155, 0.3724, ..., 0.2957, 0.2206, 0.1781]],

       [[0.0872, 0.098 , 0.0885, ..., 0.0595, 0.0877, 0.0877],
        [0.0545, 0.0569, 0.0608, ..., 0.0786, 0.1001, 0.0986],
        [0.059 , 0.059 , 0.0699, ..., 0.0832, 0.0893, 0.1054],
        ...,
...
        [0.2843, 0.2843, 0.2747, ..., 0.0335, 0.0265, 0.0854],
        [0.3272, 0.2843, 0.2747, ..., 0.0837, 0.0837, 0.0924],
        [0.3272, 0.3189, 0.3189, ..., 0.1042, 0.1042, 0.0809]],

       [[0.3649, 0.3649, 0.184 , ..., 0.0684, 0.0674, 0.0433],
        [0.0634, 0.1497, 0.103 , ..., 0.0579, 0.0684, 0.0591],
        [0.0571, 0.1222, 0.1198, ..., 0.0683, 0.0708, 0.0615],
        ...,
        [0.1463, 0.0741, 0.0715, ..., 0.0739, 0.0704, 0.0625],
        [0.2168, 0.2168, 0.2168, ..., 0.1174, 0.0796, 0.0461],
        [0.2091, 0.2091, 0.2091, ..., 0.0716, 0.078 , 0.0461]],

       [[0.5378, 0.6008, 0.6325, ..., 0.2449, 0.1424, 0.2253],
        [0.5378, 0.6008, 0.6008, ..., 0.2263, 0.1478, 0.2187],
        [0.2133, 0.2643, 0.2643, ..., 0.1357, 0.1357, 0.1593],
        ...,
        [0.9117, 0.6923, 0.2849, ..., 0.458 , 0.4329, 0.4329],
        [0.8087, 0.8237, 0.7588, ..., 0.2097, 0.2961, 0.3422],
        [0.8301, 0.8237, 0.8328, ..., 0.2015, 0.2961, 0.297 ]]],
      dtype=float32)
    • x
      PandasIndex
      PandasIndex(Index([3.1864583330478737, 3.1885416663810204,  3.190624999714167,
        3.192708333047314, 3.1947916663804605, 3.1968749997136072,
        3.198958333046754, 3.2010416663799006, 3.2031249997130473,
        3.205208333046194, 3.2072916663793407, 3.2093749997124874,
        3.211458333045634,  3.213541666378781, 3.2156249997119275,
        3.217708333045074,  3.219791666378221, 3.2218749997113676,
       3.2239583330445143,  3.226041666377661, 3.2281249997108077,
       3.2302083330439544,  3.232291666377101, 3.2343749997102478,
       3.2364583330433945,  3.238541666376541,  3.240624999709688,
       3.2427083330428346, 3.2447916663759813,  3.246874999709128,
       3.2489583330422747, 3.2510416663754214,  3.253124999708568,
       3.2552083330417148, 3.2572916663748614,  3.259374999708008,
        3.261458333041155, 3.2635416663743015, 3.2656249997074482,
        3.267708333040595, 3.2697916663737416, 3.2718749997068883,
        3.273958333040035, 3.2760416663731817, 3.2781249997063284,
        3.280208333039475,  3.282291666372622, 3.2843749997057685,
        3.286458333038915,  3.288541666372062, 3.2906249997052086,
       3.2927083330383553,  3.294791666371502, 3.2968749997046487,
       3.2989583330377954,  3.301041666370942],
      dtype='float64', name='x'))
    • y
      PandasIndex
      PandasIndex(Index([ 6.705208332732646,  6.703124999399499,  6.701041666066352,
        6.698958332733206,  6.696874999400059,  6.694791666066912,
       6.6927083327337655,  6.690624999400619,  6.688541666067472,
        6.686458332734325,  6.684374999401179,  6.682291666068032,
        6.680208332734885,  6.678124999401739,  6.676041666068592,
        6.673958332735445,  6.671874999402299,  6.669791666069152,
        6.667708332736005, 6.6656249994028585,  6.663541666069712,
        6.661458332736565,  6.659374999403418,  6.657291666070272,
        6.655208332737125,  6.653124999403978,  6.651041666070832,
        6.648958332737685,  6.646874999404538, 6.6447916660713915,
        6.642708332738245,  6.640624999405098,  6.638541666071951,
        6.636458332738805,  6.634374999405658,  6.632291666072511,
        6.630208332739365,  6.628124999406218,  6.626041666073071,
       6.6239583327399245,  6.621874999406778,  6.619791666073631,
        6.617708332740484,  6.615624999407338,  6.613541666074191,
        6.611458332741044,  6.609374999407898,  6.607291666074751,
        6.605208332741604, 6.6031249994084575,  6.601041666075311,
        6.598958332742164, 6.5968749994090174,  6.594791666075871,
        6.592708332742724,  6.590624999409577,  6.588541666076431,
        6.586458332743284,  6.584374999410137,  6.582291666076991,
        6.580208332743844,  6.578124999410697, 6.5760416660775505,
        6.573958332744404,  6.571874999411257,   6.56979166607811,
        6.567708332744964,  6.565624999411817,   6.56354166607867,
        6.561458332745524,  6.559374999412377,   6.55729166607923,
       6.5552083327460835,  6.553124999412937,   6.55104166607979,
        6.548958332746643,  6.546874999413497,   6.54479166608035,
        6.542708332747203,  6.540624999414057,   6.53854166608091,
        6.536458332747763, 6.5343749994146165,   6.53229166608147,
        6.530208332748323,  6.528124999415176,   6.52604166608203,
        6.523958332748883,  6.521874999415736],
      dtype='float64', name='y'))
    • date
      PandasIndex
      PandasIndex(DatetimeIndex(['2015-06-26', '2015-07-12', '2015-07-28', '2016-06-25',
               '2016-07-11', '2016-07-27', '2017-06-26', '2017-07-12',
               '2017-07-28', '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'],
              dtype='datetime64[ns]', name='date', freq=None))
In [29]:
# Compute the difference in NDVI before and after the fire
lagos_diff = (
    lagos_da
        .sel(date=slice('2017', '2020'))
        .mean('date')
        .NDVI 
   - lagos_da
        .sel(date=slice('2015', '2017'))
        .mean('date')
        .NDVI
)
In [31]:
# Plot the difference
fig, ax = plt.subplots()

lagos_diff.plot(x='x', y='y', cmap='PiYG', robust=True, ax=ax)
alimosho.plot(edgecolor='black', color='none', ax=ax)
plt.show()
No description has been provided for this image
In [32]:
alimosho_out_gdf = (
    gpd.GeoDataFrame(geometry=alimosho.envelope)
    .overlay(alimosho, how='difference'))
alimosho_out_gdf
Out[32]:
geometry
0 MULTIPOLYGON (((3.18644 6.52276, 3.18644 6.549...
In [33]:
# Clip data to both inside and outside the boundary
ndvi_lagos_da = lagos_da.rio.clip(alimosho.geometry, from_disk=True)
ndvi_lagos_out_da = lagos_da.rio.clip(alimosho_out_gdf.geometry, from_disk=True)
In [36]:
# Compute mean annual July NDVI
jul_ndvi_lagos_df = (
    ndvi_lagos_da
    .groupby(ndvi_lagos_da.date.dt.year)
    .mean(...)
    .NDVI.to_dataframe())
jul_ndvi_lagos_out_df = (
    ndvi_lagos_out_da
    .groupby(ndvi_lagos_out_da.date.dt.year)
    .mean(...)
    .NDVI.to_dataframe())

# Plot inside and outside the reservation
lagos_jul_ndvi_df = (
    jul_ndvi_lagos_df[['NDVI']]
    .join(
        jul_ndvi_lagos_out_df[['NDVI']], 
        lsuffix=' Inside Alimosho', rsuffix=' Outside Alimosho')
)
In [37]:
# Plot difference inside and outside the reservation
lagos_jul_ndvi_df.plot(
    title='NDVI before and after 2017; Alimosho, Lagos, Nigeria')
Out[37]:
<Axes: title={'center': 'NDVI before and after 2017; Alimosho, Lagos, Nigeria'}, xlabel='year'>
No description has been provided for this image

Rivers, Nigeria¶

In [3]:
# Select out Rivers
rivers = nigeria_gdf[nigeria_gdf["STATE"]=="Rivers"]
rivers
Out[3]:
STATE LGA AREA PERIMETER LONGITUDE LATITUDE FULL_NAME geometry
652 Rivers Ogba/ Egbema/ Ndoni 700.062 184.440 NaN NaN NaN POLYGON ((6.66876 5.69434, 6.66190 5.68314, 6....
687 Rivers Emuoha 762.646 161.086 NaN NaN NaN POLYGON ((6.80440 5.21118, 6.80077 5.20519, 6....
690 Rivers Ikwerre 673.248 117.630 NaN NaN NaN POLYGON ((6.91888 5.19793, 6.91591 5.19311, 6....
691 Rivers Ahoada West 440.862 115.787 NaN NaN NaN POLYGON ((6.58052 5.20643, 6.58048 5.20639, 6....
693 Rivers Ahoada East 354.096 114.092 NaN NaN NaN POLYGON ((6.70663 5.18277, 6.70695 5.18098, 6....
694 Rivers Etche 568.879 121.100 NaN NaN NaN POLYGON ((7.00468 5.18696, 7.00668 5.18656, 7....
695 Rivers Omumma 302.128 72.275 NaN NaN NaN POLYGON ((7.25150 5.19474, 7.26031 5.19401, 7....
714 Rivers Abua /Odual 673.838 140.800 NaN NaN NaN POLYGON ((6.65027 4.90263, 6.65031 4.90244, 6....
718 Rivers Obio/Akpor 280.085 83.299 NaN NaN NaN POLYGON ((7.08107 4.89117, 7.08948 4.88486, 7....
720 Rivers Oyigbo 342.323 107.718 NaN NaN NaN POLYGON ((7.15481 4.91673, 7.15315 4.90538, 7....
723 Rivers Degema 1023.248 223.875 NaN NaN NaN POLYGON ((6.77979 4.79679, 6.78072 4.79332, 6....
724 Rivers Eleme 125.840 46.559 NaN NaN NaN POLYGON ((7.14285 4.83191, 7.14649 4.82998, 7....
726 Rivers Khana 520.661 103.493 NaN NaN NaN POLYGON ((7.49638 4.81743, 7.49414 4.81443, 7....
727 Rivers Port Harcourt 194.759 64.004 NaN NaN NaN POLYGON ((7.09815 4.80543, 7.09758 4.80108, 7....
728 Rivers Asari-Toru 99.479 55.917 NaN NaN NaN POLYGON ((6.87581 4.77949, 6.87508 4.77951, 6....
729 Rivers Tai 101.697 45.200 NaN NaN NaN POLYGON ((7.30234 4.76531, 7.30180 4.76260, 7....
732 Rivers Okrika 226.901 89.791 NaN NaN NaN POLYGON ((7.19800 4.69843, 7.19764 4.69777, 7....
735 Rivers Akuku-Toru 533.696 125.944 NaN NaN NaN POLYGON ((6.73779 4.73185, 6.74073 4.73109, 6....
736 Rivers Ogu /Bolo 67.644 48.186 NaN NaN NaN POLYGON ((7.21447 4.71585, 7.21419 4.70832, 7....
737 Rivers Gokana 93.218 43.954 NaN NaN NaN POLYGON ((7.33225 4.70824, 7.33472 4.70121, 7....
743 Rivers Bonny 647.422 139.810 NaN NaN NaN POLYGON ((7.31188 4.61225, 7.30808 4.60857, 7....
744 Rivers Andoni 235.360 79.493 NaN NaN NaN POLYGON ((7.42821 4.60147, 7.42838 4.60124, 7....
745 Rivers Opobo /Nkoro 167.605 62.395 NaN NaN NaN POLYGON ((7.52345 4.56744, 7.52331 4.56502, 7....
In [ ]:
rivers.explore()
In [5]:
# Dissolve into single polygon
rivers_dissolve = rivers.dissolve(by='STATE')
rivers_dissolve
Out[5]:
geometry LGA AREA PERIMETER LONGITUDE LATITUDE FULL_NAME
STATE
Rivers POLYGON ((6.75750 4.34611, 6.73802 4.34505, 6.... Ogba/ Egbema/ Ndoni 700.062 184.44 NaN NaN NaN
In [6]:
rivers_dir = os.path.join(pathlib.Path.home(), 'rivers-data')
# Make the data directory
os.makedirs(rivers_dir, exist_ok=True)

rivers_dir
Out[6]:
'/home/jovyan/rivers-data'
In [7]:
# Initialize AppeearsDownloader for MODIS NDVI data
ndvi_downloader = eaapp.AppeearsDownloader(
    download_key='rivers-ndvi',
    ea_dir=rivers_dir,
    product='MOD13Q1.061',
    layer='_250m_16_days_NDVI',
    start_date="07-01",
    end_date="07-31",
    recurring=True,
    year_range=[2016, 2020],
    polygon=rivers_dissolve
)

ndvi_downloader.download_files(cache=True)
In [10]:
# Get a list of NDVI tif file paths
rivers_paths = sorted(glob(os.path.join(rivers_dir, 'rivers-ndvi', '*', '*NDVI*.tif')))
len(rivers_paths)
Out[10]:
15
In [11]:
rivers_das = []
for rivers_path in rivers_paths:
    # Get date from file name
    doy = rivers_path[doy_start:doy_end]
    date = pd.to_datetime(doy, format='%Y%j')

    # Open dataset
    da = rxr.open_rasterio(rivers_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
    rivers_das.append(da)

len(rivers_das)
Out[11]:
15
In [12]:
rivers_da = xr.combine_by_coords(rivers_das, coords=['date'])
rivers_da
Out[12]:
<xarray.Dataset>
Dimensions:      (x: 577, y: 661, date: 15)
Coordinates:
    band         int64 1
  * x            (x) float64 6.384 6.386 6.389 6.391 ... 7.578 7.58 7.582 7.584
  * y            (y) float64 5.72 5.718 5.716 5.714 ... 4.351 4.349 4.347 4.345
    spatial_ref  int64 0
  * date         (date) datetime64[ns] 2016-06-25 2016-07-11 ... 2020-07-27
Data variables:
    NDVI         (date, y, x) float32 0.8904 0.465 0.6473 0.8786 ... nan nan nan
xarray.Dataset
    • x: 577
    • y: 661
    • date: 15
    • band
      ()
      int64
      1
      array(1)
    • x
      (x)
      float64
      6.384 6.386 6.389 ... 7.582 7.584
      array([6.384375, 6.386458, 6.388542, ..., 7.580208, 7.582292, 7.584375])
    • y
      (y)
      float64
      5.72 5.718 5.716 ... 4.347 4.345
      array([5.719792, 5.717708, 5.715625, ..., 4.348958, 4.346875, 4.344792])
    • 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 :
      6.383333332761481 0.0020833333331466974 0.0 5.720833332820831 0.0 -0.0020833333331466974
      array(0)
    • date
      (date)
      datetime64[ns]
      2016-06-25 ... 2020-07-27
      array(['2016-06-25T00:00:00.000000000', '2016-07-11T00:00:00.000000000',
       '2016-07-27T00:00:00.000000000', '2017-06-26T00:00:00.000000000',
       '2017-07-12T00:00:00.000000000', '2017-07-28T00:00:00.000000000',
       '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'], dtype='datetime64[ns]')
    • NDVI
      (date, y, x)
      float32
      0.8904 0.465 0.6473 ... nan nan nan
      array([[[ 0.8904,  0.465 ,  0.6473, ...,  0.8777,  0.8834,  0.8834],
        [ 0.8901,  0.5729,  0.5296, ...,  0.8881,  0.8875,  0.8837],
        [ 0.9437,  0.5729,  0.8546, ...,  0.8881,  0.8875,  0.8816],
        ...,
        [ 0.1625,     nan,     nan, ...,     nan,     nan,     nan],
        [ 0.1322,     nan,     nan, ...,     nan,     nan,     nan],
        [ 0.15  ,  0.2023,  0.2023, ...,     nan,     nan,     nan]],

       [[ 0.2261,  0.2287,  0.2287, ...,  0.9381,  0.7785,  0.9479],
        [ 0.2096,  0.1555,  0.1555, ...,  0.8059,  0.8059,  0.9444],
        [ 0.3129,  0.1355,  0.154 , ...,  0.9446,  0.9446,  0.9507],
        ...,
        [ 0.1146,  0.0942,  0.0948, ...,     nan,     nan,     nan],
        [ 0.1163,  0.1032,  0.0948, ...,     nan,     nan,     nan],
        [ 0.1163,  0.0347,  0.028 , ...,     nan,     nan,     nan]],

       [[ 0.5496,  0.5179,  0.1984, ...,  0.2043,  0.2279,  0.5842],
        [ 0.674 ,  0.5179,  0.2581, ...,  0.1644,  0.187 ,  0.3623],
        [ 0.2718,  0.6089,  0.5185, ...,  0.2138,  0.3592,  0.2947],
        ...,
...
        [ 0.0736, -0.1373, -0.1107, ...,     nan,     nan,     nan],
        [ 0.3357,  0.2139,  0.017 , ...,     nan,     nan,     nan],
        [ 0.3861,  0.6458,  0.2188, ...,     nan,     nan,     nan]],

       [[ 0.2493,  0.2044,  0.2044, ...,  0.8304,  0.9198,  0.8775],
        [ 0.2493,  0.2428,  0.2428, ...,  0.8738,  0.8738,  0.9311],
        [ 0.3582,  0.4052,  0.4052, ...,  0.8614,  0.9454,  0.9203],
        ...,
        [ 0.1338,     nan,     nan, ...,     nan,     nan,     nan],
        [ 0.1338,     nan,     nan, ...,     nan,     nan,     nan],
        [ 0.2489,     nan,     nan, ...,     nan,     nan,     nan]],

       [[ 0.39  ,  0.4789,  0.545 , ...,  0.8239,  0.8106,  0.8173],
        [ 0.2882,  0.4789,  0.3315, ...,  0.8145,  0.8368,  0.7928],
        [ 0.2888,  0.1862,  0.1734, ...,  0.756 ,  0.756 ,  0.1663],
        ...,
        [ 0.2513,     nan,     nan, ...,     nan,     nan,     nan],
        [ 0.1007,     nan,     nan, ...,     nan,     nan,     nan],
        [    nan,     nan,     nan, ...,     nan,     nan,     nan]]],
      dtype=float32)
    • x
      PandasIndex
      PandasIndex(Index([ 6.384374999428054,  6.386458332761201,  6.388541666094348,
        6.390624999427494,  6.392708332760641,  6.394791666093788,
       6.3968749994269345,  6.398958332760081,  6.401041666093228,
        6.403124999426375,
       ...
        7.565624999322232, 7.5677083326553785,  7.569791665988525,
        7.571874999321672, 7.5739583326548185,  7.576041665987965,
        7.578124999321112,  7.580208332654259,  7.582291665987405,
        7.584374999320552],
      dtype='float64', name='x', length=577))
    • y
      PandasIndex
      PandasIndex(Index([ 5.719791666154258,  5.717708332821111,  5.715624999487964,
        5.713541666154818,  5.711458332821671,  5.709374999488524,
        5.707291666155378,  5.705208332822231,  5.703124999489084,
       5.7010416661559375,
       ...
        4.363541666275758,  4.361458332942611,  4.359374999609464,
        4.357291666276318,  4.355208332943171,  4.353124999610024,
        4.351041666276878,  4.348958332943731,  4.346874999610584,
       4.3447916662774375],
      dtype='float64', name='y', length=661))
    • date
      PandasIndex
      PandasIndex(DatetimeIndex(['2016-06-25', '2016-07-11', '2016-07-27', '2017-06-26',
               '2017-07-12', '2017-07-28', '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'],
              dtype='datetime64[ns]', name='date', freq=None))
In [13]:
# Compute the difference in NDVI before and after the fire
rivers_diff = (
    rivers_da
        .sel(date=slice('2018', '2020'))
        .mean('date')
        .NDVI 
   - rivers_da
        .sel(date=slice('2016', '2018'))
        .mean('date')
        .NDVI
)
In [17]:
# Plot the difference
fig, ax = plt.subplots(figsize=(8,10))

rivers_diff.plot(x='x', y='y', cmap='PiYG', robust=True, ax=ax)
rivers_dissolve.plot(edgecolor='black', color='none', ax=ax, linewidth=3)
plt.show()
No description has been provided for this image
In [20]:
rivers_out_gdf = (
    gpd.GeoDataFrame(geometry=rivers_dissolve.envelope)
    .overlay(rivers, how='difference'))
rivers_out_gdf
Out[20]:
geometry
0 MULTIPOLYGON (((6.38367 4.82114, 6.38367 5.718...
In [21]:
# Clip data to both inside and outside the boundary
ndvi_rivers_da = rivers_da.rio.clip(rivers_dissolve.geometry, from_disk=True)
ndvi_rivers_out_da = rivers_da.rio.clip(rivers_out_gdf.geometry, from_disk=True)
In [22]:
# Compute mean annual July NDVI
jul_ndvi_rivers_df = (
    ndvi_rivers_da
    .groupby(ndvi_rivers_da.date.dt.year)
    .mean(...)
    .NDVI.to_dataframe())
jul_ndvi_rivers_out_df = (
    ndvi_rivers_out_da
    .groupby(ndvi_rivers_out_da.date.dt.year)
    .mean(...)
    .NDVI.to_dataframe())

# Plot inside and outside the reservation
rivers_jul_ndvi_df = (
    jul_ndvi_rivers_df[['NDVI']]
    .join(
        jul_ndvi_rivers_out_df[['NDVI']], 
        lsuffix=' Inside Rivers', rsuffix=' Outside Rivers')
)
In [23]:
# Plot difference inside and outside the reservation
rivers_jul_ndvi_df.plot(
    title='NDVI before and after 2018; Rivers, Nigeria')
Out[23]:
<Axes: title={'center': 'NDVI before and after 2018; Rivers, Nigeria'}, xlabel='year'>
No description has been provided for this image

Looking at Cedar Rapids, IA¶

In [4]:
cr_url = "https://sos.iowa.gov/shapefiles/City%20Precincts/Cedar%20Rapids%20Precinct%20Boundaries.zip"
cr_gdf = gpd.read_file(cr_url)
cr_gdf.explore()
Out[4]:
Make this Notebook Trusted to load map: File -> Trust Notebook
In [5]:
# Dissolve into single polygon
cr_gdf['new_column'] = 0
cr_gdf_dissolve = cr_gdf.dissolve(by='new_column')
cr_gdf_dissolve
Out[5]:
geometry LONGNAME SHORTNAME DISTRICT COLOR TOTAL TARGET_DEV
new_column
0 POLYGON ((-91.64995 41.90342, -91.65247 41.903... Cedar Rapids 01 CR01 1 -211298291 3466 -34
In [6]:
cr_dir = os.path.join(pathlib.Path.home(), 'cr-data')
# Make the data directory
os.makedirs(cr_dir, exist_ok=True)

cr_dir
Out[6]:
'/home/jovyan/cr-data'
In [57]:
# Initialize AppeearsDownloader for MODIS NDVI data
ndvi_downloader = eaapp.AppeearsDownloader(
    download_key='cr-ndvi',
    ea_dir=cr_dir,
    product='MOD13Q1.061',
    layer='_250m_16_days_NDVI',
    start_date="07-01",
    end_date="07-31",
    recurring=True,
    year_range=[2018, 2023],
    polygon=cr_gdf_dissolve
)

ndvi_downloader.download_files(cache=True)
In [7]:
# Get a list of NDVI tif file paths
cr_paths = sorted(glob(os.path.join(cr_dir, 'cr-ndvi', '*', '*NDVI*.tif')))
len(cr_paths)
Out[7]:
18
In [11]:
cr_das = []
for cr_path in cr_paths:
    # Get date from file name
    doy = cr_path[doy_start:doy_end]
    date = pd.to_datetime(doy, format='%Y%j')

    # Open dataset
    da = rxr.open_rasterio(cr_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
    cr_das.append(da)

len(cr_das)
Out[11]:
18
In [12]:
cr_da = xr.combine_by_coords(cr_das, coords=['date'])
cr_da
Out[12]:
<xarray.Dataset>
Dimensions:      (x: 108, y: 100, date: 18)
Coordinates:
    band         int64 1
  * x            (x) float64 -91.77 -91.77 -91.77 ... -91.56 -91.55 -91.55
  * y            (y) float64 42.07 42.07 42.06 42.06 ... 41.87 41.87 41.86 41.86
    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.7344 0.7344 0.719 ... 0.7285 0.7285
xarray.Dataset
    • x: 108
    • y: 100
    • date: 18
    • band
      ()
      int64
      1
      array(1)
    • x
      (x)
      float64
      -91.77 -91.77 ... -91.55 -91.55
      array([-91.773958, -91.771875, -91.769792, -91.767708, -91.765625, -91.763542,
       -91.761458, -91.759375, -91.757292, -91.755208, -91.753125, -91.751042,
       -91.748958, -91.746875, -91.744792, -91.742708, -91.740625, -91.738542,
       -91.736458, -91.734375, -91.732292, -91.730208, -91.728125, -91.726042,
       -91.723958, -91.721875, -91.719792, -91.717708, -91.715625, -91.713542,
       -91.711458, -91.709375, -91.707292, -91.705208, -91.703125, -91.701042,
       -91.698958, -91.696875, -91.694792, -91.692708, -91.690625, -91.688542,
       -91.686458, -91.684375, -91.682292, -91.680208, -91.678125, -91.676042,
       -91.673958, -91.671875, -91.669792, -91.667708, -91.665625, -91.663542,
       -91.661458, -91.659375, -91.657292, -91.655208, -91.653125, -91.651042,
       -91.648958, -91.646875, -91.644792, -91.642708, -91.640625, -91.638542,
       -91.636458, -91.634375, -91.632292, -91.630208, -91.628125, -91.626042,
       -91.623958, -91.621875, -91.619792, -91.617708, -91.615625, -91.613542,
       -91.611458, -91.609375, -91.607292, -91.605208, -91.603125, -91.601042,
       -91.598958, -91.596875, -91.594792, -91.592708, -91.590625, -91.588542,
       -91.586458, -91.584375, -91.582292, -91.580208, -91.578125, -91.576042,
       -91.573958, -91.571875, -91.569792, -91.567708, -91.565625, -91.563542,
       -91.561458, -91.559375, -91.557292, -91.555208, -91.553125, -91.551042])
    • y
      (y)
      float64
      42.07 42.07 42.06 ... 41.86 41.86
      array([42.067708, 42.065625, 42.063542, 42.061458, 42.059375, 42.057292,
       42.055208, 42.053125, 42.051042, 42.048958, 42.046875, 42.044792,
       42.042708, 42.040625, 42.038542, 42.036458, 42.034375, 42.032292,
       42.030208, 42.028125, 42.026042, 42.023958, 42.021875, 42.019792,
       42.017708, 42.015625, 42.013542, 42.011458, 42.009375, 42.007292,
       42.005208, 42.003125, 42.001042, 41.998958, 41.996875, 41.994792,
       41.992708, 41.990625, 41.988542, 41.986458, 41.984375, 41.982292,
       41.980208, 41.978125, 41.976042, 41.973958, 41.971875, 41.969792,
       41.967708, 41.965625, 41.963542, 41.961458, 41.959375, 41.957292,
       41.955208, 41.953125, 41.951042, 41.948958, 41.946875, 41.944792,
       41.942708, 41.940625, 41.938542, 41.936458, 41.934375, 41.932292,
       41.930208, 41.928125, 41.926042, 41.923958, 41.921875, 41.919792,
       41.917708, 41.915625, 41.913542, 41.911458, 41.909375, 41.907292,
       41.905208, 41.903125, 41.901042, 41.898958, 41.896875, 41.894792,
       41.892708, 41.890625, 41.888542, 41.886458, 41.884375, 41.882292,
       41.880208, 41.878125, 41.876042, 41.873958, 41.871875, 41.869792,
       41.867708, 41.865625, 41.863542, 41.861458])
    • 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 :
      -91.77499999146531 0.002083333333139592 0.0 42.06874999608778 0.0 -0.002083333333139592
      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]')
    • NDVI
      (date, y, x)
      float32
      0.7344 0.7344 ... 0.7285 0.7285
      array([[[0.7344, 0.7344, 0.719 , ..., 0.7983, 0.7983, 0.7383],
        [0.8021, 0.8021, 0.8197, ..., 0.7721, 0.7721, 0.7555],
        [0.8526, 0.8526, 0.8416, ..., 0.8549, 0.8549, 0.7907],
        ...,
        [0.8978, 0.901 , 0.901 , ..., 0.718 , 0.6464, 0.6464],
        [0.8978, 0.8978, 0.9028, ..., 0.633 , 0.718 , 0.718 ],
        [0.8513, 0.8513, 0.831 , ..., 0.7464, 0.5944, 0.5944]],

       [[0.7968, 0.7968, 0.8048, ..., 0.8797, 0.8797, 0.8948],
        [0.7804, 0.7804, 0.8205, ..., 0.8964, 0.8964, 0.8964],
        [0.8693, 0.8693, 0.873 , ..., 0.909 , 0.909 , 0.8964],
        ...,
        [0.9098, 0.9351, 0.9351, ..., 0.8329, 0.6793, 0.6793],
        [0.9123, 0.9123, 0.9385, ..., 0.7642, 0.6793, 0.6793],
        [0.8743, 0.8743, 0.8781, ..., 0.7642, 0.6374, 0.6374]],

       [[0.7546, 0.7546, 0.7418, ..., 0.9215, 0.9215, 0.8677],
        [0.7947, 0.7947, 0.7947, ..., 0.9296, 0.9296, 0.8768],
        [0.8141, 0.8141, 0.8356, ..., 0.9296, 0.9296, 0.8768],
        ...,
...
        [0.8336, 0.8476, 0.8476, ..., 0.8438, 0.8111, 0.8111],
        [0.8369, 0.8369, 0.8583, ..., 0.8434, 0.7873, 0.7873],
        [0.832 , 0.832 , 0.8359, ..., 0.7822, 0.6777, 0.6777]],

       [[0.7369, 0.7369, 0.7158, ..., 0.8765, 0.8765, 0.8314],
        [0.8348, 0.8348, 0.7158, ..., 0.8933, 0.8933, 0.8167],
        [0.8171, 0.8171, 0.8348, ..., 0.9142, 0.9142, 0.8167],
        ...,
        [0.8865, 0.9266, 0.9266, ..., 0.8661, 0.7776, 0.7776],
        [0.8918, 0.8918, 0.9266, ..., 0.8647, 0.7407, 0.7407],
        [0.931 , 0.931 , 0.8729, ..., 0.7705, 0.7407, 0.7407]],

       [[0.7315, 0.7315, 0.7294, ..., 0.8792, 0.8792, 0.8201],
        [0.7721, 0.7721, 0.7082, ..., 0.8535, 0.8535, 0.8293],
        [0.8094, 0.8094, 0.7961, ..., 0.9014, 0.9014, 0.8095],
        ...,
        [0.922 , 0.9417, 0.9417, ..., 0.8874, 0.8326, 0.8326],
        [0.9053, 0.9053, 0.9053, ..., 0.8874, 0.7962, 0.7962],
        [0.9297, 0.9297, 0.9053, ..., 0.7298, 0.7285, 0.7285]]],
      dtype=float32)
    • x
      PandasIndex
      PandasIndex(Index([-91.77395832479874, -91.77187499146561, -91.76979165813246,
       -91.76770832479932, -91.76562499146618, -91.76354165813305,
        -91.7614583247999, -91.75937499146676, -91.75729165813362,
       -91.75520832480049,
       ...
       -91.56979165815106, -91.56770832481791, -91.56562499148478,
       -91.56354165815165,  -91.5614583248185, -91.55937499148536,
       -91.55729165815222, -91.55520832481909, -91.55312499148594,
        -91.5510416581528],
      dtype='float64', name='x', length=108))
    • y
      PandasIndex
      PandasIndex(Index([ 42.06770832942121,  42.06562499608807,  42.06354166275493,
        42.06145832942179,  42.05937499608865,  42.05729166275551,
        42.05520832942237,  42.05312499608923,  42.05104166275609,
        42.04895832942295,  42.04687499608981,  42.04479166275667,
        42.04270832942353,  42.04062499609039, 42.038541662757254,
       42.036458329424114, 42.034374996090975, 42.032291662757835,
       42.030208329424696, 42.028124996091556, 42.026041662758416,
        42.02395832942528,  42.02187499609214,    42.019791662759,
        42.01770832942586,  42.01562499609272,  42.01354166275958,
        42.01145832942644,   42.0093749960933,  42.00729166276016,
        42.00520832942702,  42.00312499609388,  42.00104166276074,
         41.9989583294276,  41.99687499609446,  41.99479166276132,
        41.99270832942818,  41.99062499609504, 41.988541662761904,
       41.986458329428764, 41.984374996095625, 41.982291662762485,
       41.980208329429345, 41.978124996096206, 41.976041662763066,
        41.97395832942993,  41.97187499609679,  41.96979166276365,
        41.96770832943051,  41.96562499609737,  41.96354166276423,
        41.96145832943109,  41.95937499609795,  41.95729166276481,
        41.95520832943167,  41.95312499609853,  41.95104166276539,
        41.94895832943225,  41.94687499609911,  41.94479166276597,
        41.94270832943283,  41.94062499609969,  41.93854166276655,
       41.936458329433414, 41.934374996100274, 41.932291662767135,
       41.930208329433995, 41.928124996100856, 41.926041662767716,
       41.923958329434576,  41.92187499610144,   41.9197916627683,
        41.91770832943516,  41.91562499610202,  41.91354166276888,
        41.91145832943574,   41.9093749961026,  41.90729166276946,
        41.90520832943632,  41.90312499610318,  41.90104166277004,
         41.8989583294369,  41.89687499610376,  41.89479166277062,
        41.89270832943748,  41.89062499610434,   41.8885416627712,
       41.886458329438064, 41.884374996104924, 41.882291662771785,
       41.880208329438645, 41.878124996105505, 41.876041662772366,
       41.873958329439226,  41.87187499610609,  41.86979166277295,
        41.86770832943981,  41.86562499610667,  41.86354166277353,
        41.86145832944039],
      dtype='float64', name='y'))
    • 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))
In [13]:
# Compute the difference in NDVI before and after the fire
cr_diff = (
    cr_da
        .sel(date=slice('2020', '2023'))
        .mean('date')
        .NDVI 
   - cr_da
        .sel(date=slice('2018', '2020'))
        .mean('date')
        .NDVI
)
In [25]:
# Plot the difference
fig, ax = plt.subplots()

cr_diff.plot(x='x', y='y', cmap='PiYG', robust=True, ax=ax)
cr_gdf_dissolve.plot(edgecolor='black', color='none', ax=ax)
plt.title('NDVI Difference (July) @ Cedar Rapids, IA\nPost/Pre 2020 Derecho')
plt.show()
No description has been provided for this image
In [21]:
cr_out_gdf = (
    gpd.GeoDataFrame(geometry=cr_gdf_dissolve.envelope)
    .overlay(cr_gdf_dissolve, how='difference'))
cr_out_gdf
Out[21]:
geometry
0 MULTIPOLYGON (((-91.77459 41.86159, -91.77459 ...
In [22]:
# Clip data to both inside and outside the boundary
ndvi_cr_da = cr_da.rio.clip(cr_gdf_dissolve.geometry, from_disk=True)
ndvi_cr_out_da = cr_da.rio.clip(cr_out_gdf.geometry, from_disk=True)
In [23]:
# Compute mean annual July NDVI
jul_ndvi_cr_df = (
    ndvi_cr_da
    .groupby(ndvi_cr_da.date.dt.year)
    .mean(...)
    .NDVI.to_dataframe())
jul_ndvi_cr_out_df = (
    ndvi_cr_out_da
    .groupby(ndvi_cr_out_da.date.dt.year)
    .mean(...)
    .NDVI.to_dataframe())

# Plot inside and outside the reservation
derecho_ndvi_df = (
    jul_ndvi_cr_df[['NDVI']]
    .join(
        jul_ndvi_cr_out_df[['NDVI']], 
        lsuffix=' Inside Cedar Rapids', rsuffix=' Outside Cedar Rapids')
)
In [24]:
# Plot difference inside and outside the reservation
derecho_ndvi_df.plot(
    title='NDVI before and after the 2020 Derecho')
Out[24]:
<Axes: title={'center': 'NDVI before and after the 2020 Derecho'}, xlabel='year'>
No description has been provided for this image

Glacier National Park¶

Looking at NDVI changes at Glacier National Park (Montana, USA) 2015-2020

In [2]:
# Read in all National Park Boundaries (US)
nps_boundary = gpd.read_file("http://gstore.unm.edu/apps/rgisarchive/"
                             "datasets/7bbe8af5-029b-4adf-b06c-134f0dd57226/"
                             "nps_boundary.original.zip")
nps_boundary
Out[2]:
GIS_LOC_ID UNIT_CODE GROUP_CODE UNIT_NAME UNIT_TYPE META_MIDF LANDS_CODE DATE_EDIT GIS_NOTES geometry
0 NaN UPDE NaN Upper Delaware Scenic and Recreational River NaN NaN 2006-06-06 POC for this update: leslie_morlock@nps.gov -... POLYGON ((-74.92112 41.49030, -74.91918 41.491...
1 NaN BAND NaN Bandelier National Monument NaN NaN NaN Park Created Boundary - contact Kay Beeley MULTIPOLYGON (((-106.19585 35.85378, -106.2214...
2 NaN ACAD NaN Acadia National Park NaN NaN NaN Shift 0.03 mi. MULTIPOLYGON (((-68.22595 44.39584, -68.22766 ...
3 NaN ROCR NaN Achbold Parkway Park NaN NaN NaN NaN MULTIPOLYGON (((-77.08430 38.91673, -77.08482 ...
4 NaN ADAM NaN Adams National Historical Park NaN NaN NaN Shift 0.03 mi. POLYGON ((-71.00948 42.25625, -71.01133 42.255...
... ... ... ... ... ... ... ... ... ... ...
505 NaN CLBA NaN Clara Barton National Historic Site NaN NaN 2008-11-17 Lands - http://landsnet.nps.gov/tractsnet/docu... POLYGON ((-77.14072 38.96699, -77.14075 38.967...
506 NaN FOWA NaN Fort Washington Park Park NaN NaN 2008-11-17 Lands - http://landsnet.nps.gov/tractsnet/docu... POLYGON ((-77.02299 38.72054, -77.02032 38.718...
507 NaN KICA NaN Kings Canyon National Park NaN NaN 2003-11-07 Lands MULTIPOLYGON (((-118.77919 37.23609, -118.7791...
508 NaN SEQU NaN Sequoia National Park NaN NaN 2004-01-09 Lands POLYGON ((-118.70637 36.66475, -118.70636 36.6...
509 NaN CABR NaN Cabrillo National Monument NaN NaN 2008-08-28 Lands - http://landsnet.nps.gov/tractsnet/docu... POLYGON ((-117.23862 32.66936, -117.24036 32.6...

510 rows × 10 columns

In [3]:
# Select out Glacier National Park
glacier = nps_boundary[nps_boundary["UNIT_NAME"] == "Glacier"]
glacier
Out[3]:
GIS_LOC_ID UNIT_CODE GROUP_CODE UNIT_NAME UNIT_TYPE META_MIDF LANDS_CODE DATE_EDIT GIS_NOTES geometry
361 NaN GLAC NaN Glacier National Park NaN NaN 2006-03-09 POC for this update: richard_menicke@nps.gov ... POLYGON ((-114.47552 49.00091, -114.46101 49.0...
In [4]:
# Plot the boundary to inspect
glacier.explore()
Out[4]:
Make this Notebook Trusted to load map: File -> Trust Notebook
In [5]:
glacier_dir = os.path.join(pathlib.Path.home(), 'glacier-data')
# Make the data directory
os.makedirs(glacier_dir, exist_ok=True)

glacier_dir
Out[5]:
'/home/jovyan/glacier-data'
In [8]:
# Initialize AppeearsDownloader for MODIS NDVI data
ndvi_downloader = eaapp.AppeearsDownloader(
    download_key='glacier-ndvi',
    ea_dir=glacier_dir,
    product='MOD13Q1.061',
    layer='_250m_16_days_NDVI',
    start_date="07-01",
    end_date="07-31",
    recurring=True,
    year_range=[2015, 2020],
    polygon=glacier
)

ndvi_downloader.download_files(cache=True)
In [9]:
# Get a list of NDVI tif file paths
glacier_paths = sorted(glob(os.path.join(glacier_dir, 'glacier-ndvi', '*', '*NDVI*.tif')))
len(glacier_paths)
Out[9]:
18
In [11]:
# Build data list of data arrays
glacier_das = []
for glacier_path in glacier_paths:
    # Get date from file name
    doy = glacier_path[doy_start:doy_end]
    date = pd.to_datetime(doy, format='%Y%j')

    # Open dataset
    da = rxr.open_rasterio(glacier_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
    glacier_das.append(da)

len(glacier_das)
Out[11]:
18
In [12]:
# Combine into single array by date
glacier_da = xr.combine_by_coords(glacier_das, coords=['date'])
glacier_da
Out[12]:
<xarray.Dataset>
Dimensions:      (x: 593, y: 369, date: 18)
Coordinates:
    band         int64 1
  * x            (x) float64 -114.5 -114.5 -114.5 ... -113.2 -113.2 -113.2
  * y            (y) float64 49.0 49.0 49.0 48.99 ... 48.24 48.24 48.24 48.23
    spatial_ref  int64 0
  * date         (date) datetime64[ns] 2015-06-26 2015-07-12 ... 2020-07-27
Data variables:
    NDVI         (date, y, x) float32 0.5447 0.5447 0.6585 ... 0.7097 0.7097
xarray.Dataset
    • x: 593
    • y: 369
    • date: 18
    • band
      ()
      int64
      1
      array(1)
    • x
      (x)
      float64
      -114.5 -114.5 ... -113.2 -113.2
      array([-114.476042, -114.473958, -114.471875, ..., -113.246875, -113.244792,
       -113.242708])
    • y
      (y)
      float64
      49.0 49.0 49.0 ... 48.24 48.23
      array([49.001042, 48.998958, 48.996875, ..., 48.238542, 48.236458, 48.234375])
    • 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 :
      -114.47708332307788 0.0020833333331466974 0.0 49.00208332894347 0.0 -0.0020833333331466974
      array(0)
    • date
      (date)
      datetime64[ns]
      2015-06-26 ... 2020-07-27
      array(['2015-06-26T00:00:00.000000000', '2015-07-12T00:00:00.000000000',
       '2015-07-28T00:00:00.000000000', '2016-06-25T00:00:00.000000000',
       '2016-07-11T00:00:00.000000000', '2016-07-27T00:00:00.000000000',
       '2017-06-26T00:00:00.000000000', '2017-07-12T00:00:00.000000000',
       '2017-07-28T00:00:00.000000000', '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'],
      dtype='datetime64[ns]')
    • NDVI
      (date, y, x)
      float32
      0.5447 0.5447 ... 0.7097 0.7097
      array([[[0.5447, 0.5447, 0.6585, ..., 0.8173, 0.832 , 0.832 ],
        [0.5769, 0.5408, 0.561 , ..., 0.7763, 0.7725, 0.832 ],
        [0.6659, 0.6659, 0.5408, ..., 0.7578, 0.7725, 0.7725],
        ...,
        [0.7366, 0.7622, 0.7622, ..., 0.7493, 0.7206, 0.7206],
        [0.7533, 0.7533, 0.7777, ..., 0.6806, 0.6806, 0.7206],
        [0.7139, 0.7636, 0.7636, ..., 0.5371, 0.6806, 0.6806]],

       [[0.6447, 0.6447, 0.6695, ..., 0.7528, 0.7857, 0.7857],
        [0.5996, 0.5996, 0.669 , ..., 0.7421, 0.7623, 0.7857],
        [0.7458, 0.7458, 0.5996, ..., 0.7421, 0.7623, 0.7623],
        ...,
        [0.7816, 0.7809, 0.7809, ..., 0.7702, 0.7488, 0.7488],
        [0.7841, 0.7841, 0.7809, ..., 0.7335, 0.7335, 0.693 ],
        [0.7675, 0.7816, 0.7816, ..., 0.6231, 0.6546, 0.6546]],

       [[0.5584, 0.5584, 0.6954, ..., 0.7219, 0.692 , 0.692 ],
        [0.5573, 0.4936, 0.6954, ..., 0.6881, 0.7219, 0.7325],
        [0.6886, 0.6886, 0.5887, ..., 0.7167, 0.7322, 0.7322],
        ...,
...
        [0.771 , 0.7731, 0.7731, ..., 0.7174, 0.684 , 0.684 ],
        [0.771 , 0.771 , 0.7769, ..., 0.6597, 0.6597, 0.5493],
        [0.7613, 0.7769, 0.7769, ..., 0.704 , 0.6564, 0.6564]],

       [[0.6542, 0.6542, 0.76  , ..., 0.7848, 0.8278, 0.8278],
        [0.6807, 0.6666, 0.763 , ..., 0.8117, 0.78  , 0.8278],
        [0.8111, 0.8111, 0.7676, ..., 0.7813, 0.8058, 0.8058],
        ...,
        [0.7831, 0.8209, 0.8209, ..., 0.8113, 0.8042, 0.8042],
        [0.8254, 0.8254, 0.8209, ..., 0.7148, 0.7148, 0.7134],
        [0.8184, 0.8254, 0.8254, ..., 0.6909, 0.6972, 0.6972]],

       [[0.627 , 0.627 , 0.73  , ..., 0.7044, 0.7158, 0.7158],
        [0.6401, 0.6102, 0.7485, ..., 0.7206, 0.7366, 0.812 ],
        [0.682 , 0.682 , 0.6695, ..., 0.7449, 0.77  , 0.77  ],
        ...,
        [0.7775, 0.803 , 0.803 , ..., 0.7913, 0.7678, 0.7678],
        [0.8184, 0.8184, 0.8158, ..., 0.7507, 0.7507, 0.7616],
        [0.7887, 0.8097, 0.8097, ..., 0.5577, 0.7097, 0.7097]]],
      dtype=float32)
    • x
      PandasIndex
      PandasIndex(Index([ -114.4760416564113, -114.47395832307815,   -114.471874989745,
       -114.46979165641186, -114.46770832307871, -114.46562498974556,
       -114.46354165641242, -114.46145832307927, -114.45937498974612,
       -114.45729165641298,
       ...
       -113.26145832318677, -113.25937498985363, -113.25729165652048,
       -113.25520832318733, -113.25312498985419, -113.25104165652104,
       -113.24895832318789, -113.24687498985475,  -113.2447916565216,
       -113.24270832318845],
      dtype='float64', name='x', length=593))
    • y
      PandasIndex
      PandasIndex(Index([  49.0010416622769,  48.99895832894375, 48.996874995610604,
        48.99479166227746,  48.99270832894431, 48.990624995611164,
        48.98854166227802,  48.98645832894487, 48.984374995611724,
        48.98229166227858,
       ...
        48.25312499567723, 48.251041662344086,  48.24895832901094,
        48.24687499567779, 48.244791662344646,   48.2427083290115,
        48.24062499567835, 48.238541662345206,  48.23645832901206,
        48.23437499567891],
      dtype='float64', name='y', length=369))
    • date
      PandasIndex
      PandasIndex(DatetimeIndex(['2015-06-26', '2015-07-12', '2015-07-28', '2016-06-25',
               '2016-07-11', '2016-07-27', '2017-06-26', '2017-07-12',
               '2017-07-28', '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'],
              dtype='datetime64[ns]', name='date', freq=None))
In [14]:
# Compute the difference in NDVI before and 2017/2018
glacier_diff = (
    glacier_da
        .sel(date=slice('2018', '2020'))
        .mean('date')
        .NDVI 
   - glacier_da
        .sel(date=slice('2015', '2017'))
        .mean('date')
        .NDVI
)

glacier_diff
Out[14]:
<xarray.DataArray 'NDVI' (y: 369, x: 593)>
array([[ 0.00964445,  0.00964445, -0.00370002, ...,  0.07951105,
         0.09175545,  0.09175545],
       [ 0.03782219,  0.03242224,  0.01788896, ...,  0.09422225,
         0.07023335,  0.0790112 ],
       [ 0.03799999,  0.03799999,  0.04133326, ...,  0.06477785,
         0.06571102,  0.06571102],
       ...,
       [ 0.00105554,  0.00936663,  0.00936663, ..., -0.00990003,
        -0.02216661, -0.02216661],
       [ 0.01837778,  0.01837778,  0.00757778, ..., -0.03024441,
        -0.03024441, -0.03652221],
       [ 0.01888889,  0.01841098,  0.01841098, ..., -0.04922229,
        -0.01131105, -0.01131105]], dtype=float32)
Coordinates:
    band         int64 1
  * x            (x) float64 -114.5 -114.5 -114.5 ... -113.2 -113.2 -113.2
  * y            (y) float64 49.0 49.0 49.0 48.99 ... 48.24 48.24 48.24 48.23
    spatial_ref  int64 0
xarray.DataArray
'NDVI'
  • y: 369
  • x: 593
  • 0.009644 0.009644 -0.0037 -0.001745 ... -0.04922 -0.01131 -0.01131
    array([[ 0.00964445,  0.00964445, -0.00370002, ...,  0.07951105,
         0.09175545,  0.09175545],
       [ 0.03782219,  0.03242224,  0.01788896, ...,  0.09422225,
         0.07023335,  0.0790112 ],
       [ 0.03799999,  0.03799999,  0.04133326, ...,  0.06477785,
         0.06571102,  0.06571102],
       ...,
       [ 0.00105554,  0.00936663,  0.00936663, ..., -0.00990003,
        -0.02216661, -0.02216661],
       [ 0.01837778,  0.01837778,  0.00757778, ..., -0.03024441,
        -0.03024441, -0.03652221],
       [ 0.01888889,  0.01841098,  0.01841098, ..., -0.04922229,
        -0.01131105, -0.01131105]], dtype=float32)
    • band
      ()
      int64
      1
      array(1)
    • x
      (x)
      float64
      -114.5 -114.5 ... -113.2 -113.2
      array([-114.476042, -114.473958, -114.471875, ..., -113.246875, -113.244792,
       -113.242708])
    • y
      (y)
      float64
      49.0 49.0 49.0 ... 48.24 48.23
      array([49.001042, 48.998958, 48.996875, ..., 48.238542, 48.236458, 48.234375])
    • 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 :
      -114.47708332307788 0.0020833333331466974 0.0 49.00208332894347 0.0 -0.0020833333331466974
      array(0)
    • x
      PandasIndex
      PandasIndex(Index([ -114.4760416564113, -114.47395832307815,   -114.471874989745,
       -114.46979165641186, -114.46770832307871, -114.46562498974556,
       -114.46354165641242, -114.46145832307927, -114.45937498974612,
       -114.45729165641298,
       ...
       -113.26145832318677, -113.25937498985363, -113.25729165652048,
       -113.25520832318733, -113.25312498985419, -113.25104165652104,
       -113.24895832318789, -113.24687498985475,  -113.2447916565216,
       -113.24270832318845],
      dtype='float64', name='x', length=593))
    • y
      PandasIndex
      PandasIndex(Index([  49.0010416622769,  48.99895832894375, 48.996874995610604,
        48.99479166227746,  48.99270832894431, 48.990624995611164,
        48.98854166227802,  48.98645832894487, 48.984374995611724,
        48.98229166227858,
       ...
        48.25312499567723, 48.251041662344086,  48.24895832901094,
        48.24687499567779, 48.244791662344646,   48.2427083290115,
        48.24062499567835, 48.238541662345206,  48.23645832901206,
        48.23437499567891],
      dtype='float64', name='y', length=369))
In [ ]:
In [16]:
# Plot the difference
fig, ax = plt.subplots()

glacier_diff.plot(x='x', y='y', cmap='PiYG', robust=True, ax=ax)
glacier.plot(edgecolor='black', color='none', ax=ax)
plt.title('NDVI Difference (July) @ Glacier National Park\n2018-2020 vs. 2015-2017')
plt.show()
No description has been provided for this image
In [17]:
glacier_out_gdf = (
    gpd.GeoDataFrame(geometry=glacier.envelope)
    .overlay(glacier, how='difference'))
glacier_out_gdf
Out[17]:
geometry
0 MULTIPOLYGON (((-114.47552 48.23369, -114.4755...
In [18]:
# Clip data to both inside and outside the boundary
ndvi_glacier_da = glacier_da.rio.clip(glacier.geometry, from_disk=True)
ndvi_glacier_out_da = glacier_da.rio.clip(glacier_out_gdf.geometry, from_disk=True)
In [19]:
# Compute mean annual July NDVI
jul_ndvi_glacier_df = (
    ndvi_glacier_da
    .groupby(ndvi_glacier_da.date.dt.year)
    .mean(...)
    .NDVI.to_dataframe())
jul_ndvi_glacier_out_df = (
    ndvi_glacier_out_da
    .groupby(ndvi_glacier_out_da.date.dt.year)
    .mean(...)
    .NDVI.to_dataframe())

# Plot inside and outside the reservation
glacier_ndvi_df = (
    jul_ndvi_glacier_df[['NDVI']]
    .join(
        jul_ndvi_glacier_out_df[['NDVI']], 
        lsuffix=' Inside Glacier National Park', rsuffix=' Outside Glacier National Park')
)
In [20]:
# Plot difference inside and outside the reservation
glacier_ndvi_df.plot(
    title='NDVI Differences at Glacier National Park\n2018-2020 vs. 2015-2107')
Out[20]:
<Axes: title={'center': 'NDVI Differences at Glacier National Park\n2018-2020 vs. 2015-2107'}, xlabel='year'>
No description has been provided for this image
In [ ]:
%%capture
%%bash
jupyter nbconvert *.ipynb --to html