Reprojecting datacube and raster data

• Products used: gm_s2_annual

• Special requirements: This notebook loads data from an external raster file (madagascar_CHPclim_02.tif) from the Supplementary_data folder of this repository

Keywords data used; sentinel-2 geomedian, data used; CHIRPS, data methods; reprojecting, data methods; resampling, data format; GeoTIFF, data methods; geobox

Background

It is often valuable to combine data from the datacube with other external raster datasets. For example, we may want to use a rainfall raster to focus our analysis on satellite data to mask areas of higher or lower precipitation. However, it can be challenging to combine datasets if they have different extents, resolutions (e.g. 20 m vs. 250 m), or coordinate reference systems (e.g. WGS 84 vs. Africa Albers Equal Area Conic). To be able to combine these datasets, we need to be able to reproject them into identical spatial grids prior to analysis.

Datacube stores information about the location, resolution and coordinate reference system of a rectangular grid of data using an object called a GeoBox. This GeoBox object is added to all data loaded from the datacube, and all rasters loaded with xr.open_rasterio (provided that import datacube is run before the raster is loaded). Datacube functions can use this object to provide a template that can be used to reproject raster and datacube data - either applying this reprojection directly when new data is being loaded, or to reproject existing data that has already been loaded into memory. This makes it straightforward to reproject one dataset into the exact spatial grid of another dataset for further analysis.

Description

This notebook demonstrates how to perform three key reprojection workflows when working with datacube data and external rasters:

1. Loading and reprojecting datacube data to match a raster file using dc.load

2. Reprojecting existing datacube data to match a raster using xr_reproject

3. Loading and reprojecting a raster to match datacube data using rio_slurp_xarray

---

Getting started

To run this analysis, run all the cells in the notebook, starting with the “Load packages” cell.

Load packages

Import Python packages that are used for the analysis.

import datacube
import rioxarray
import matplotlib.pyplot as plt
from odc.algo import xr_reproject
from datacube.testutils.io import rio_slurp_xarray
from datacube.utils.masking import mask_invalid_data

from deafrica_tools.plotting import rgb

Connect to the datacube

Connect to the datacube so we can access DE Africa data. The app parameter is a unique name for the analysis which is based on the notebook file name.

dc = datacube.Datacube(app="Reprojecting_data")

Loading and reprojecting datacube data to match a raster

Load a raster file

First we load a GeoTIFF raster from file using rioxarray.open_rasterio. The example below loads a long-term mean rainfall data for the month of February as according to Climate Hazards Group InfraRed Precipitation with Station (CHIRPS) for a section of Madagascar. It has a spatial resolution of ~0.05 decimal degrees (approximately equal to 5 km pixels at the equator) in the WGS 84 (EPSG:4326) coordinate reference system.

# Path to raster file
raster_path = "../Supplementary_data/Reprojecting_data/madagascar_CHPclim_02.tif"

# Load raster, and remove redundant "band" dimension
raster = rioxarray.open_rasterio(raster_path).squeeze("band", drop=True)

print(raster)

<xarray.DataArray (y: 20, x: 20)>
[400 values with dtype=float32]
Coordinates:
  * x            (x) float64 46.03 46.08 46.13 46.18 ... 46.83 46.88 46.93 46.98
  * y            (y) float64 -18.03 -18.08 -18.13 ... -18.88 -18.93 -18.98
    spatial_ref  int64 0
Attributes:
    AREA_OR_POINT:  Area
    scale_factor:   1.0
    add_offset:     0.0

If we plot our loaded raster, we can see there are areas with higher rainfall, and areas with lower rainfall.

raster.plot(size=6)

<matplotlib.collections.QuadMesh at 0x7fa4bdd79e70>

GeoBox objects

Now we have loaded our raster dataset, we can inspect its GeoBox object that we will use to allow us to reproject data. The GeoBox can be accessed using the .geobox method. It includes a set of information that together completely define the spatial grid of our data: * The width (e.g. 20) and height (e.g. 20) of our data in pixels * An Affine object which defines the spatial resolution (e.g. -0.05000000074505806 and 0.05000000074505806) and spatial position (e.g. 46.00000336766243 and -18.00000160932541) of our data * The coordinate reference system of our data (e.g. +init=epsg:4326)

raster.geobox

GeoBox(20, 20, Affine(0.05000000074505806, 0.0, 46.00000336766243,
       0.0, -0.05000000074505806, -18.00000160932541), 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"]])

The GeoBox also has a number of useful methods that can be used to view information about the spatial grid of our data. For example, we can inspect the spatial resolution of our data:

raster.geobox.resolution

(-0.05000000074505806, 0.05000000074505806)

Or obtain information about the data’s spatial extent:

raster.geobox.extent.boundingbox

BoundingBox(left=46.00000336766243, bottom=-19.00000162422657, right=47.00000338256359, top=-18.00000160932541)

Note: For more information about GeoBox objects and a complete list of their methods, refer to the datacube documentation.

Load and reproject datacube data

We can now use datacube to load and reproject satellite data to exactly match the coordinate reference system and resolution of our raster data. By specifying like=raster.geobox, we can load datacube data that will be reprojected to match the spatial grid of our raster.

# Load and reproject data from datacube
ds = dc.load(product="gm_s2_annual",
             measurements=["red", "green", "blue"],
             time=("2018"),
             like=raster.geobox,
             resampling="nearest",
             group_by="solar_day")

When we plot the result, it should appear similarly pixelated to the raster we loaded above, which has a lower resolution.

rgb(ds, col='time')

We can also directly compare the geobox of the two datasets to verify they share the same spatial grid:

ds.geobox == raster.geobox

True

Now that our two datasets share the same spatial grid, we can use our raster as a mask. For example, we can mask out all satellite pixels except those with higher rainfall (e.g. more than 350 mm/month of rain):

# Rename raster dimensions to match datacube conventions for data
# with geographic coordinates
if raster.geobox.crs.geographic:
    raster = raster.rename({"x": "longitude", "y": "latitude"})

# Identify high rainfall areas
high_rainfall = raster > 350

# Apply mask to set lower-rainfall areas to `NaN`
ds_masked = ds.where(high_rainfall)

# Plot the masked data
rgb(ds_masked)

/usr/local/lib/python3.10/dist-packages/matplotlib/cm.py:478: RuntimeWarning: invalid value encountered in cast
  xx = (xx * 255).astype(np.uint8)

Reprojecting existing datacube data to match a raster

The example above demonstrated how to load new satellite data from the datacube to match the spatial grid of a raster. However, sometimes we may have already loaded satellite data with a coordinate reference system and resolution that is different from our raster. In this case, we need to reproject this existing dataset to match our raster.

For example, we have loaded the 2018 Sentinel-2 geomedian data from the datacube with 20 m-resolution pixels in the EPSG:6933 coordinate reference system. Note that in this example we manually specify the x, y, resolution and output_crs parameters, rather than taking them directly from our raster using like=raster.geobox in the previous example.

# Load data from datacube
ds = dc.load(product="gm_s2_annual",
             measurements=["red", "green", "blue"],
             time=("2018"),
             x=(46.00, 47.00),
             y=(-19, -18),
             resolution=(-20, 20),
             output_crs="EPSG:6933",
             group_by="solar_day")

# Load raster, and remove redundant "band" dimension
raster = rioxarray.open_rasterio(raster_path).squeeze("band", drop=True)

If we plot our satellite data, we can see that it is much higher resolution than our pixelated CHIRPS precipitation raster:

rgb(ds)

Reproject datacube data

We can now use the xr_reproject function to reproject our existing high resolution satellite dataset. We specify geobox=raster.geobox to request that the data gets reprojected to match the spatial grid of our low resolution raster.

To control how the data is reprojected, we can specify a custom resampling method that will control how our high resolution pixels will be transformed into lower resolution pixels. In this case, we specify "average", which will set the value of each larger pixel to the average of all smaller pixels that fall within its pixel boundary.

Note: Refer to the rasterio documentation for a full list of available resampling methods.

# Temporary workaround for bug in `xr_reproject`
if raster.geobox.crs.geographic:
    ds = ds.rename({"x": "longitude", "y": "latitude"})

# Reproject data
ds_reprojected = xr_reproject(src=ds,
                              geobox=raster.geobox,
                              resampling="average")

#Set nodata to `NaN`
ds_reprojected = mask_invalid_data(ds_reprojected)

Now if we plot our reprojected dataset, we can see that our satellite imagery now has a similar resolution to our low resolution raster (e.g. with a pixelated appearance). Note however, that this result will appear smoother than the previous example due to the "average" resampling method used here (compared to "nearest").

rgb(ds_reprojected)

Once again, we can also verify that the two datasets have identical spatial grids:

ds_reprojected.geobox == raster.geobox

True

Loading and reprojecting a raster to match datacube data

Rather than reprojecting satellite data to match the resolution and projection system of our raster, we may instead wish to reproject our raster to match the spatial grid of our satellite data. This can be particularly useful when we have a lower resolution raster file (e.g. like the ~5 km resolution CHIRPS data we are using in this example), but we don’t want to lose the better spatial resolution of our satellite data.

Load datacube data

As in the previous example, we can load in satellite data from the datacube at 20 m spatial resolution and EPSG:6933 coordinate reference system:

# Load data from datacube
ds = dc.load(product="gm_s2_annual",
             measurements=["red", "green", "blue"],
             time=("2018"),
             x=(46.00, 47.00),
             y=(-19, -18),
             resolution=(-20, 20),
             output_crs="EPSG:6933",
             group_by="solar_day")

# Plot 20 m resolution satellite data
rgb(ds)

Load and reproject raster data

We can now use the rio_slurp_xarray function to load and reproject our raster file to match our higher resolution satellite dataset. We specify gbox=ds.geobox to request that our raster gets reprojected to match the spatial grid of ds. As in the previous xr_reproject example, we can also specify a custom resampling method which will be used during the reprojection. In this case, we specify 'bilinear', which will produce a smooth output without obvious pixel boundaries.

Note: Refer to the rasterio documentation for a full list of available resampling methods.

# Load raster and reproject to match satellite dataset
raster_reprojected = rio_slurp_xarray(fname=raster_path,
                                      gbox=ds.geobox,
                                      resampling="bilinear")

# Set nodata to `NaN`
raster_reprojected = mask_invalid_data(raster_reprojected)

If we plot our resampled raster data, it should now appear less pixelated than the original ~250 m resolution raster.

raster_reprojected.plot(size=6, vmin=200);

The resampled raster should also match the spatial grid of our higher resolution satellite data:

raster_reprojected.geobox == ds.geobox

True

Now both of our datasets share the same spatial grid, we can use our resampled raster to mask our higher resolution satellite dataset as we did in the first section (e.g. mask out all areas with less than 350 mm of rainfall in February). Compared to the first example, this masked satellite dataset should appear much higher resolution, with far less obvious pixelation:

# Identify high rainfall areas
high_rainfall = raster_reprojected > 350

# Apply mask to set high rainfall areas to `NaN`
ds_masked = ds.where(high_rainfall)

# Plot the masked data
rgb(ds_masked)

/usr/local/lib/python3.10/dist-packages/matplotlib/cm.py:478: RuntimeWarning: invalid value encountered in cast
  xx = (xx * 255).astype(np.uint8)

---

Additional information

License: The code in this notebook is licensed under the Apache License, Version 2.0. Digital Earth Africa data is licensed under the Creative Commons by Attribution 4.0 license.

Contact: If you need assistance, please post a question on the Open Data Cube Slack channel or on the GIS Stack Exchange using the open-data-cube tag (you can view previously asked questions here). If you would like to report an issue with this notebook, you can file one on Github.

Compatible datacube version:

print(datacube.__version__)

1.8.15

Last Tested:

from datetime import datetime
datetime.today().strftime('%Y-%m-%d')

'2023-08-14'