In [1]:
# Install the development version of the earthpy package
# !pip install git+https://github.com/earthlab/earthpy@apppears
In [2]:
import getpass
import json
import os
import pathlib
from glob import glob
import matplotlib.pyplot as plt
# 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
In [3]:
abuja_dir = os.path.join(pathlib.Path.home(), 'abuja-data')
# Make the data directory
os.makedirs(abuja_dir, exist_ok=True)
abuja_dir
Out[3]:
'/home/jovyan/abuja-data'
In [4]:
# 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[4]:
| STATE | LGA | AREA | PERIMETER | LONGITUDE | LATITUDE | FULL_NAME | geometry | |
|---|---|---|---|---|---|---|---|---|
| 0 | Sokoto | Gada | 1193.977 | 170.095 | None | None | None | POLYGON ((5.53632 13.88793, 5.53480 13.88488, ... |
| 1 | Sokoto | Illela | 1298.423 | 174.726 | None | None | None | POLYGON ((5.53632 13.88793, 5.54517 13.88419, ... |
| 2 | Sokoto | Tangaza | 2460.715 | 209.702 | None | None | None | POLYGON ((4.85548 13.76724, 4.86189 13.78085, ... |
| 3 | Borno | Abadam | 2430.515 | 288.957 | None | None | None | POLYGON ((12.83189 13.39871, 12.83397 13.40439... |
| 4 | Lake | Lake chad | 5225.912 | 497.039 | None | None | None | POLYGON ((13.48608 13.30821, 13.48296 13.31344... |
| ... | ... | ... | ... | ... | ... | ... | ... | ... |
| 770 | Delta | Isoko North | 485.467 | 169.369 | None | None | None | MULTIPOLYGON (((6.31996 5.63341, 6.32003 5.633... |
| 771 | Niger | Lavun | 3951.431 | 424.153 | None | None | None | MULTIPOLYGON (((6.12188 9.09441, 6.12001 9.094... |
| 772 | Yobe | Bade | 817.260 | 216.207 | None | None | None | MULTIPOLYGON (((11.01052 12.80457, 11.00747 12... |
| 773 | Zamfara | Maru | 7795.261 | 536.500 | None | None | None | MULTIPOLYGON (((6.43894 12.41104, 6.43609 12.4... |
| 774 | Akwa Ibom | Oron | 81.472 | 57.846 | None | None | None | POLYGON ((8.22063 4.84473, 8.23405 4.82974, 8.... |
775 rows × 8 columns
In [8]:
print(nigeria_gdf.columns)
Index(['STATE', 'LGA', 'AREA', 'PERIMETER', 'LONGITUDE', 'LATITUDE',
'FULL_NAME', 'geometry'],
dtype='object')
In [14]:
# SEARCH FOR STATE IN DATA FRAME.
search_state = 'Abuja' # Note: Use capital letters if the column is in capital letters
# Check if the state exists in the specified column
if search_state in nigeria_gdf['STATE'].values:
print(f"{search_state} exists in the GeoDataFrame.")
else:
print(f"{search_state} does not exist in the GeoDataFrame.")
Abuja exists in the GeoDataFrame.
In [16]:
# Select out Abuja
Abuja = nigeria_gdf[nigeria_gdf["STATE"]=="Abuja"]
Abuja
Out[16]:
| STATE | LGA | AREA | PERIMETER | LONGITUDE | LATITUDE | FULL_NAME | geometry | |
|---|---|---|---|---|---|---|---|---|
| 300 | Abuja | Bwari | 818.909 | 153.605 | None | None | None | POLYGON ((7.62143 9.30766, 7.62241 9.30840, 7.... |
| 307 | Abuja | Gwagwalada | 967.762 | 129.429 | None | None | None | POLYGON ((7.12020 9.12689, 7.12679 9.11870, 7.... |
| 308 | Abuja | Abaji | 1002.453 | 257.578 | None | None | None | POLYGON ((6.96319 9.23240, 6.96163 9.23123, 6.... |
| 314 | Abuja | Abuja Municipal | 1585.320 | 206.501 | None | None | None | POLYGON ((7.17027 9.09925, 7.17036 9.09829, 7.... |
| 329 | Abuja | Kwali | 1299.111 | 151.853 | None | None | None | POLYGON ((7.05974 8.89200, 7.06064 8.89094, 7.... |
| 331 | Abuja | Kuje | 1708.466 | 217.281 | None | None | None | POLYGON ((7.07285 8.90097, 7.08316 8.90024, 7.... |
In [17]:
# Display Interactive map of state, In this case (Abuja)
Abuja.explore()
Out[17]:
Make this Notebook Trusted to load map: File -> Trust Notebook
In [21]:
# Select Gwagwalada
Gwagwalada = Abuja[Abuja["LGA"] == 'Gwagwalada']
Gwagwalada
Out[21]:
| STATE | LGA | AREA | PERIMETER | LONGITUDE | LATITUDE | FULL_NAME | geometry | |
|---|---|---|---|---|---|---|---|---|
| 307 | Abuja | Gwagwalada | 967.762 | 129.429 | None | None | None | POLYGON ((7.12020 9.12689, 7.12679 9.11870, 7.... |
In [22]:
# Display Gwagwalada Map
Gwagwalada.explore()
Out[22]:
Make this Notebook Trusted to load map: File -> Trust Notebook
In [23]:
# Initialize AppeearsDownloader for MODIS NDVI data
ndvi_downloader = eaapp.AppeearsDownloader(
download_key='abuja-ndvi',
ea_dir=abuja_dir,
product='MOD13Q1.061',
layer='_250m_16_days_NDVI',
start_date="07-01",
end_date="07-31",
recurring=True,
year_range=[2015, 2020],
polygon= Gwagwalada
)
ndvi_downloader.download_files(cache=True)
In [24]:
# Get a list of NDVI tif file paths
abuja_paths = sorted(glob(os.path.join(abuja_dir, 'abuja-ndvi', '*', '*NDVI*.tif')))
len(abuja_paths)
Out[24]:
18
In [27]:
scale_factor = 10000
doy_start = -19
doy_end = -12
In [28]:
abuja_das = []
for abuja_path in abuja_paths:
# Get date from file name
doy = abuja_path [doy_start:doy_end]
date = pd.to_datetime(doy, format='%Y%j')
# Open dataset
da = rxr.open_rasterio(abuja_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
abuja_das.append(da)
len(abuja_das)
Out[28]:
18
In [29]:
abuja_da = xr.combine_by_coords(abuja_das, coords=['date'])
abuja_da
Out[29]:
<xarray.Dataset>
Dimensions: (x: 165, y: 167, date: 18)
Coordinates:
band int64 1
* x (x) float64 6.824 6.826 6.828 6.83 ... 7.159 7.161 7.164 7.166
* y (y) float64 9.236 9.234 9.232 9.23 ... 8.897 8.895 8.893 8.891
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.5796 0.5933 0.5727 ... 0.6463 0.6531In [30]:
# Compute the difference in NDVI before and after the fire
abuja_diff = (
abuja_da
.sel(date=slice('2017', '2020'))
.mean('date')
.NDVI
- abuja_da
.sel(date=slice('2015', '2017'))
.mean('date')
.NDVI
)
In [32]:
# Plot the difference
fig, ax = plt.subplots()
abuja_diff.plot(x='x', y='y', cmap='PiYG', robust=True, ax=ax)
Gwagwalada.plot(edgecolor='black', color='none', ax=ax)
plt.show()
In [33]:
Gwagwalada_out_gdf = (
gpd.GeoDataFrame(geometry=Gwagwalada.envelope)
.overlay(Gwagwalada, how='difference'))
Gwagwalada_out_gdf
Out[33]:
| geometry | |
|---|---|
| 0 | MULTIPOLYGON (((6.82315 8.88978, 6.82315 9.109... |
In [34]:
# Clip data to both inside and outside the boundary
ndvi_abuja_da = abuja_da.rio.clip(Gwagwalada.geometry, from_disk=True)
ndvi_abuja_out_da = abuja_da.rio.clip(Gwagwalada_out_gdf.geometry, from_disk=True)
In [35]:
# Compute mean annual July NDVI
jul_ndvi_abuja_df = (
ndvi_abuja_da
.groupby(ndvi_abuja_da.date.dt.year)
.mean(...)
.NDVI.to_dataframe())
jul_ndvi_abuja_out_df = (
ndvi_abuja_out_da
.groupby(ndvi_abuja_out_da.date.dt.year)
.mean(...)
.NDVI.to_dataframe())
# Plot inside and outside the reservation
abuja_jul_ndvi_df = (
jul_ndvi_abuja_df[['NDVI']]
.join(
jul_ndvi_abuja_out_df[['NDVI']],
lsuffix=' Inside Gwagwalada', rsuffix=' Outside Gwagwalada')
)
In [36]:
# Plot difference inside and outside the reservation
abuja_jul_ndvi_df.plot(
title='NDVI before and after 2017; Gwagwalada, Abuja, Nigeria')
Out[36]:
<Axes: title={'center': 'NDVI before and after 2017; Gwagwalada, Abuja, Nigeria'}, xlabel='year'>
