Extracting Time Series at Multiple Polygons with Xarray

Extracting Time Series at Multiple Polygons with Xarray#

Introduction#

Zonal Statistics is used to summarize the values of a raster dataset within the zones of a vector dataset. We can also extend this process over time to extract data for multiple time-steps. The xvec library works with the vector geometries directly and handles a raster’s extra time dimension automatically, giving us a vector data cube of summary values for each zone at each time-step.

This is a large computation that is enabled by cloud-hosted NightTime Lights data in the Cloud-Optimized GeoTIFF (COG) format and parallel computing on a local dask cluster.

Overview of the Task#

We select all Admin1 units of a country and calculate a sum of nighttime light pixel intensities over multiple years.

Input Layers:

ne_10m_admin_1_states_provinces.zip: Admin-1 (States/Provinces) shapefile from Natural Earth
2015.tif, 2016.tif, … 2020.tif: Yearly global rasters of nighttime light data. These files were download from Harvard Dataverse and converted to Cloud-Optimized GeoTIFFs using GDAL.

Example command:
```
gdalwarp -of COG 2021_HasMask/LongNTL_2021.tif 2021.tif \
  -te -180 -90 180 90 -dstnodata 0 \
  -co COMPRESS=DEFLATE -co PREDICTOR=2 -co NUM_THREADS=ALL_CPUS
```

Output Layers:

output.csv : A CSV file with the extracted sum of nighttime lights for each admin-1 region for each year.

Data Credit:

Made with Natural Earth. Free vector and raster map data @ naturalearthdata.com.
Chen, Zuoqi; Yu, Bailang; Yang, Chengshu; Zhou, Yuyu; Yao, Shenjun; Qian, Xingjian; Wang, Congxiao; Wu, Bin; Wu, Jianping; Liao, Lingxing; Shi, Kaifang, 2020, “The global NPP-VIIRS-like nighttime light data (Version 2) for 1992-2024”, https://doi.org/10.7910/DVN/YGIVCD, Harvard Dataverse, V8

Running the Notebook: The preferred way to run this notebook is on Google Colab.

Setup and Data Download#

The following blocks of code will install the required packages and download the datasets to your Colab environment.

%%capture
if 'google.colab' in str(get_ipython()):
    !pip install rioxarray xvec exactextract

import os
import glob
import pandas as pd
import geopandas as gpd
import numpy as np
import xarray as xr
import rioxarray as rxr
import matplotlib.pyplot as plt
from datetime import datetime
import dask
import xvec

data_folder = 'data'
output_folder = 'output'

if not os.path.exists(data_folder):
    os.mkdir(data_folder)
if not os.path.exists(output_folder):
    os.mkdir(output_folder)

Setup a local Dask cluster. This distributes the computation across multiple workers on your computer.

from dask.distributed import Client, progress
client = Client()  # set up local cluster on the machine
client

If you are running this notebook in Colab, you will need to create and use a proxy URL to see the dashboard running on the local server.

if 'google.colab' in str(get_ipython()):
    from google.colab import output
    port_to_expose = 8787  # This is the default port for Dask dashboard
    print(output.eval_js(f'google.colab.kernel.proxyPort({port_to_expose})'))

def download(url):
    filename = os.path.join(data_folder, os.path.basename(url))
    if not os.path.exists(filename):
        from urllib.request import urlretrieve
        local, _ = urlretrieve(url, filename)
        print('Downloaded ' + local)

admin1_zipfile = 'ne_10m_admin_1_states_provinces.zip'
admin1_url = 'https://naciscdn.org/naturalearth/10m/cultural/'

download(admin1_url + admin1_zipfile)

Data Pre-Processing#

First we will read the Admin-1 shapefile and filter to a country. The iso_a2 column contains the 2-digit ISO code for the country. Here we are seelcting all provinces in Sri Lanka.

country_code = 'LK'
admin1_file_path = os.path.join(data_folder, admin1_zipfile)

admin1_df = gpd.read_file(admin1_file_path)

zones = admin1_df[admin1_df['iso_a2'] == country_code]
zones = zones[['adm1_code', 'name', 'iso_a2', 'geometry']].copy()
zones['id'] = zones.reset_index().index + 1

Plot the selected zones.

zones.explore()

Make this Notebook Trusted to load map: File -> Trust Notebook

Next we read the Yearly Night-time Lights Rasters and create an Xarray Dataset.

start_year = 2015
end_year = 2020

data_folder = 'https://github.com/spatialthoughts/geopython-tutorials/releases/download/data/'

bbox = zones.total_bounds

da_list = []

for year in range(start_year, end_year + 1):
  print(f'Processing year: {year}')
  cog_url = f'{data_folder}/{year}.tif'
  da = rxr.open_rasterio(
      cog_url,
      chunks=True).rio.clip_box(*bbox)
  dt = pd.to_datetime(year, format='%Y')
  da = da.assign_coords(time = dt)
  da = da.expand_dims(dim="time")
  da_list.append(da)

ntl_datacube = xr.concat(da_list, dim='time').chunk('auto')
ntl_datacube

Since all the rasters have a single band of information, we select the first band and drop that dimension.

ntl_datacube = ntl_datacube.sel(band=1, drop=True)
ntl_datacube

<xarray.DataArray (time: 6, y: 870, x: 498)> Size: 10MB
dask.array<getitem, shape=(6, 870, 498), dtype=float32, chunksize=(6, 870, 498), chunktype=numpy.ndarray>
Coordinates:
  * time         (time) datetime64[us] 48B 2015-01-01 2016-01-01 ... 2020-01-01
  * y            (y) float64 7kB 9.828 9.823 9.819 9.814 ... 5.933 5.929 5.924
  * x            (x) float64 4kB 79.66 79.66 79.66 79.67 ... 81.88 81.88 81.89
    spatial_ref  int64 8B 0
Attributes:
    DataType:            Generic
    AREA_OR_POINT:       Area
    BandName:            Band_1
    RepresentationType:  ATHEMATIC
    scale_factor:        1.0
    add_offset:          0.0
    long_name:           Band_1
    _FillValue:          0.0

Zonal Stats#

Now we will extract the sum of the raster pixel values for every admin-1 region in the selected country.

We will use xvec.zonal_stats() to summarize the raster for each zone. Unlike a rasterization-based approach, xvec works with the polygon geometries directly and handles the extra time dimension of the data cube automatically.

When doing Zonal Stats it is preferred to keep the raster in its original projection to minimize distortions, so we reproject the zone polygons to match the CRS of the raster.

zones_reprojected = zones.to_crs(ntl_datacube.rio.crs)

Next we call zonal_stats() on the data cube, passing the zone geometries. We use the exactextract method, which weights each pixel by the fraction of its area that falls within the polygon. Since the raster has a time dimension, the result is a DataArray with both time and geometry dimensions holding the sum of nighttime lights for each region for each year.

Dask will now distribute the computation across all chunks using the available workers.

It will be useful to carry over some attributes from the original GeoDataFrame. As we have an XArray vector data cube, we add them as coordinate variables along the geometry dimension.

for column in ['adm1_code', 'name', 'iso_a2']:
    result[column] = ('geometry', zones_reprojected[column].values)
    result = result.assign_coords({column: result[column]})
result

Convert the results to a Pandas DataFrame. We now have the total of nighttime lights values for each region for each year.

stats = result.drop_vars('spatial_ref', errors='ignore').to_dataframe('sum').reset_index()
stats['year'] = stats['time'].dt.year
stats.head()

	geometry	time	index	adm1_code	name	iso_a2	sum	year
0	POLYGON ((80.92292 8.98892, 80.92319 8.97321, ...	2015-01-01	1859	LKA-2454	Trikuṇāmalaya	LK	450.266036	2015
1	POLYGON ((80.92292 8.98892, 80.92319 8.97321, ...	2016-01-01	1859	LKA-2454	Trikuṇāmalaya	LK	453.701265	2016
2	POLYGON ((80.92292 8.98892, 80.92319 8.97321, ...	2017-01-01	1859	LKA-2454	Trikuṇāmalaya	LK	641.043588	2017
3	POLYGON ((80.92292 8.98892, 80.92319 8.97321, ...	2018-01-01	1859	LKA-2454	Trikuṇāmalaya	LK	657.221597	2018
4	POLYGON ((80.92292 8.98892, 80.92319 8.97321, ...	2019-01-01	1859	LKA-2454	Trikuṇāmalaya	LK	856.091326	2019

Select the columns we need for the output.

output = stats[['adm1_code', 'name', 'iso_a2', 'year', 'sum']]
output.head()

	adm1_code	name	iso_a2	year	sum
0	LKA-2454	Trikuṇāmalaya	LK	2015	450.266036
1	LKA-2454	Trikuṇāmalaya	LK	2016	453.701265
2	LKA-2454	Trikuṇāmalaya	LK	2017	641.043588
3	LKA-2454	Trikuṇāmalaya	LK	2018	657.221597
4	LKA-2454	Trikuṇāmalaya	LK	2019	856.091326

Finally, we save the result to disk.

output_file = 'output.csv'
output_path = os.path.join(output_folder, output_file)

output.to_csv(output_path, index=False)
print('Successfully written output file at {}'.format(output_path))

If you want to give feedback or share your experience with this tutorial, please comment below. (requires GitHub account)