Rioxarray

13. Rioxarray#

13.1. Overview#

rioxarray is an extension of the powerful Python library Xarray that focuses on geospatial raster data. It provides easy access to georeferencing information and geospatial transforms using Xarray’s labeled, multi-dimensional data structures, which makes it an ideal tool for working with geospatial data like satellite imagery or climate data.

The key feature of rioxarray is its seamless integration of rasterio’s geospatial data handling capabilities (such as CRS and affine transforms) with Xarray’s efficient multi-dimensional array handling. This allows you to manipulate, analyze, and visualize raster data with ease.

13.2. Learning Objectives#

By the end of this lecture, you should be able to:

Understand how rioxarray extends Xarray for geospatial data handling.
Load and inspect georeferenced raster datasets using rioxarray.
Perform basic geospatial operations, such as clipping, reprojection, and masking, using rioxarray.
Use rioxarray to manage CRS and spatial dimensions in raster datasets.
Export and visualize geospatial raster datasets.

13.3. Installing rioxarray#

To use rioxarray, you’ll need to install it along with rasterio and its dependencies. You can install it via pip by uncommenting the following line:

# %pip install rioxarray rasterio

13.4. Importing rioxarray#

You can start by importing rioxarray and other necessary libraries:

import rioxarray
import numpy as np
import xarray as xr
import matplotlib.pyplot as plt

xr.set_options(keep_attrs=True, display_expand_data=False)

<xarray.core.options.set_options at 0x7fc3ba9458b0>

13.6. Basic Geospatial Operations#

13.6.1. Reprojecting a Dataset#

Reprojecting raster data to another CRS is common in geospatial analysis. For example, you may want to reproject the dataset from its native projection to the WGS84 geographic coordinate system (EPSG:4326):

# Reproject the dataset to WGS84 (EPSG:4326)
data_reprojected = data.rio.reproject("EPSG:4326")
print(data_reprojected.rio.crs)

EPSG:4326

13.6.2. Clipping a Raster#

Clipping a raster dataset is useful when you only want to focus on a specific geographic area. You can clip a dataset using a bounding box in the same CRS as the data:

# Define a bounding box (in the same CRS as the dataset)
bbox = [-115.391, 35.982, -114.988, 36.425]

# Clip the raster to the bounding box
clipped_data = data_reprojected.rio.clip_box(*bbox)

clipped_data.shape

(7, 492, 447)

Alternatively, you can clip the raster using a vector dataset containing polygon geometries:

import geopandas as gpd

# Load a geojson with regions of interest
geojson_path = "https://github.com/opengeos/datasets/releases/download/places/las_vegas_bounds_utm.geojson"
bounds = gpd.read_file(geojson_path)

# Clip the raster to the shape
clipped_data2 = data.rio.clip(bounds.geometry, bounds.crs)

clipped_data2.shape

(7, 522, 514)

13.7. Working with Spatial Dimensions#

rioxarray supports operations on spatial dimensions (latitude/longitude or x/y coordinates) like resampling, reducing, or slicing.

13.7.1. Resampling#

To resample the raster dataset to a different resolution (e.g., 1 km), use the rio.resample method:

# Resample to 1km resolution (using an average resampling method)
resampled_data = data.rio.reproject(data.rio.crs, resolution=(1000, 1000))

resampled_data.shape

(7, 231, 224)

13.7.2. Extracting Spatial Subsets#

You can extract spatial subsets of the dataset by selecting specific coordinate ranges:

# Select a subset of the data within a lat/lon range
min_x, max_x = -115.391, -114.988
min_y, max_y = 35.982, 36.425
subset = data_reprojected.sel(
    x=slice(min_x, max_x), y=slice(max_y, min_y)
)  # Slice y in reverse order
subset.shape

(7, 491, 447)

13.8. Visualization of Georeferenced Data#

Once you have performed operations on the data, you can visualize it using matplotlib. For example, to plot a multi-band image using bands 4, 3, and 2:

# Plot the raster data
plt.figure(figsize=(8, 8))
data_reprojected.sel(band=[4, 3, 2]).plot.imshow(vmin=0, vmax=0.3)
plt.title("Landsat Image covering Las Vegas")
plt.xlabel("Longitude")
plt.ylabel("Latitude")
plt.show()

../../_images/125a2a9dea5fb975919a040430486f08d185773fba5b05cf11a6410bf11fee10.png

You can also visualize clipped or masked data in the same way:

# Plot the raster data
plt.figure(figsize=(8, 8))
clipped_data.sel(band=[4, 3, 2]).plot.imshow(vmin=0, vmax=0.3)
plt.title("Landsat Image covering Las Vegas")
plt.xlabel("Longitude")
plt.ylabel("Latitude")
plt.show()

../../_images/0a464cfab016e81ebb7ca9ec75c30938bfcd67a84a9137a43d262e4ac0594dec.png

For more advanced plots, such as overlaying a vector dataset on the raster data, you can combine rioxarray with geopandas and matplotlib:

# Plot raster with GeoJSON overlay
fig, ax = plt.subplots(figsize=(8, 8))
data.attrs["long_name"] = "Surface Reflectance"  # Update the long_name attribute
data.sel(band=4).plot.imshow(ax=ax, vmin=0, vmax=0.4, cmap="gray")
bounds.boundary.plot(ax=ax, color="red")
plt.title("Raster with Vector Overlay")
plt.show()

../../_images/79ee821c008a74947587075c76bb05c0941724b0ee59513192a2f477b6f56d7d.png

13.9. Saving Data#

Just like loading data, you can export rioxarray datasets to disk. For example, you can save the modified or processed raster data as a GeoTIFF file:

# Save the DataArray as a GeoTIFF file
data.rio.to_raster("output_raster.tif")

13.10. Handling NoData Values#

If your dataset contains NoData values, you can manage them using the following functions:

# Assign NoData value
data2 = data.rio.set_nodata(-9999)

# Remove NoData values (mask them)
data_clean = data2.rio.write_nodata(-9999, inplace=True)

13.11. Reproject to Multiple CRS#

You can reproject the dataset to multiple CRS and compare them. For instance:

# Reproject to WGS 84 (EPSG:4326)
data = data.rio.reproject("EPSG:4326")
print(data.rio.crs)

EPSG:4326

# Reproject to EPSG:3857 (Web Mercator)
mercator_data = data.rio.reproject("EPSG:3857")
print(mercator_data.rio.crs)

EPSG:3857

# Plot the raster data in WGS84
plt.figure(figsize=(6, 6))
data.sel(band=[4, 3, 2]).plot.imshow(vmin=0, vmax=0.3)
plt.title("EPSG:4326")
plt.xlabel("Longitude")
plt.ylabel("Latitude")
plt.show()

../../_images/63513d7ca8b786c6205bb6cc4d888bb92bb65ac4bb91a87cc2d35e3c5a8374de.png

# Plot the raster data in Web Mercator
plt.figure(figsize=(6, 6))
mercator_data.sel(band=[4, 3, 2]).plot.imshow(vmin=0, vmax=0.3)
plt.title("EPSG:3857")
plt.xlabel("X")
plt.ylabel("Y")
plt.show()

../../_images/d1057f824fd9b3507e7fe748c054ffaf2df2bbb1336cb4a6381076cebffea705.png

13.12. Basic Band Math (NDVI Calculation)#

Band math enables us to perform computations across different bands. A common application is calculating the Normalized Difference Vegetation Index (NDVI), which is an indicator of vegetation health.

NDVI is calculated as:

NDVI = (NIR - Red) / (NIR + Red)

We can compute and plot the NDVI as follows:

# Select the red (band 4) and NIR (band 5) bands
red_band = data.sel(band=4)
nir_band = data.sel(band=5)

# Calculate NDVI
ndvi = (nir_band - red_band) / (nir_band + red_band)
ndvi = ndvi.clip(min=-1, max=1)  # Clip values to the range [-1, 1]
ndvi.attrs["long_name"] = "NDVI"

To visualize the NDVI, we can plot it using matplotlib:

# Plot the NDVI values
ndvi.plot(cmap="RdYlGn", vmin=-1, vmax=1)
plt.title("NDVI of the Landsat Image")
plt.xlabel("Longitude")
plt.ylabel("Latitude")
plt.show()

../../_images/5e5a0f008e666ca1a79f9468fea44f045914ce4ca13c3f1751bc11483fc9c093.png

You can also mask out non-vegetated areas or areas with invalid NDVI values (such as water or urban regions) by applying a threshold:

# Mask out non-vegetated areas (NDVI < 0.2)
ndvi_clean = ndvi.where(ndvi > 0.2)
ndvi_clean.plot(cmap="Greens", vmin=0.2, vmax=0.5)
plt.title("Cleaned NDVI (non-vegetated areas masked)")
plt.xlabel("Longitude")
plt.ylabel("Latitude")
plt.show()

../../_images/75f40f282a6770053bfb099fb03d015beb44951fbd22f9b9b26e0c99e1b58314.png

13.13. Exercises#

13.13.1. Sample Dataset#

For the exercises, we will use a sample GeoTIFF raster dataset of Libya, which is available at the following URL: opengeos/datasets

13.13.2. Exercise 1: Load and Inspect a Raster Dataset#

Use rioxarray to load the GeoTIFF raster file.
Inspect the dataset by printing its dimensions, coordinates, and attributes.
Check and print the CRS and affine transformation of the dataset.

13.13.3. Exercise 2: Reproject the Raster to a New CRS#

Reproject the loaded raster dataset from its original CRS to EPSG:4326 (WGS84).
Print the new CRS and check the dimensions and coordinates of the reprojected data.
Plot the original and reprojected datasets for comparison.

13.13.4. Exercise 3: Clip the Raster Using a Bounding Box#

Define a bounding box (e.g., xmin, ymin, xmax, ymax) that covers the land area of Libya.
Clip the raster dataset using this bounding box.
Plot the clipped data to visualize the result.

13.13.5. Exercise 4: Mask the Raster Using a Vector Dataset#

Load the GeoJSON file at opengeos/datasets using geopandas.
Use the GeoJSON to mask the raster dataset, keeping only the data within the GeoJSON boundaries.
Plot the masked raster data.

13.13.6. Exercise 5: Resample the Raster to a Different Resolution#

Resample the raster dataset to a 3m resolution, using an average resampling method.
Check the new dimensions and coordinates after resampling.
Save the resampled raster dataset as a new GeoTIFF file.

13.14. Summary#

In this lecture, we have explored the basic functionality of rioxarray, a powerful extension of Xarray designed for geospatial raster data. Key points include:

Loading and inspecting geospatial raster data with CRS and transform metadata.
Performing essential geospatial operations, such as reprojection, clipping, and masking.
Visualizing and exporting raster data.
Working with spatial dimensions (x, y) using slicing, resampling, and other operations.

By integrating Xarray’s multi-dimensional data handling capabilities with rasterio’s geospatial features, rioxarray makes it easier to manage and analyze geospatial raster datasets. It is a versatile tool for anyone working with geospatial data in scientific computing, environmental analysis, or remote sensing.