How to Use Xarray in Python for GIS and Raster Data Processing
Modern GIS and remote sensing often require multidimensional rasters. Data from satellites, climate grids, weather conditions, elevation models, and time series can be multidimensional. This means that Xarray for Python can be useful for working with such kinds of rasters with labeled dimensions and coordinates.
GIS specialists and remote sensing experts may find Xarray beneficial for working with multidimensional rasters.
What is Xarray?
Xarray is an open-source Python library used for handling labeled multi-dimensional arrays. This library builds on top of the functionality provided by the NumPy library by adding labels for dimensions, coordinates, and other metadata.
As opposed to using just positions of arrays such as rows and columns, Xarray enables working with meaningful dimensions such as:
x and y - for geographic coordinates
latitude and longitude
time - for time series datasets
band - for multi-spectral or hyper-spectral images
These features are especially helpful when working with geospatial raster datasets since spatial and temporal metadata are linked to their numeric data.
Why Use Xarray for GIS and Raster Processing?
While conventional NumPy arrays work effectively for mathematical operations, they are devoid of any notion of geographical coordinates or dimensions. Xarray brings a concept of a data model that makes multidimensional raster analysis easier and more understandable.
The applications of Xarray are as follows:
Analysis of satellite images
Processing of multispectral raster
Time series analysis
Climate/weather data
Processing of digital elevation models
Land cover analysis
Calculation of vegetation index
Analysis of raster data
Filtering in space and time
Large multidimensional data sets
Installing Xarray
The installation of Xarray is done via pip:
pip install xarray
For geospatial raster processes, the installation of rioxarray may be useful:
pip install rioxarray
rioxarray is an extension to Xarray that introduces functionality related to raster data processing based on Rasterio, such as coordinate reference system management, spatial transformations, clipping, reprojections, and exporting raster data.
Understanding Xarray Data Structures
There are three main data structures used by Xarray: DataArray, Dataset, and coordinates.
DataArray
DataArray represents a labeled multi-dimensional array. For instance, a single raster layer might contain y and x dimensions.
import xarray as xr
import numpy as np
data = xrDataArray(
np.randomrand(100, 100),
dims=("y", "x"),
name="elevation"
)
print(data)
Here, there are 100 rows and 100 columns, while Xarray explicitly labels them as y and x.
Opening Raster Data with rioxarray
rioxarray comes in handy in GIS operations to open geospatial rasters.
import rioxarray
raster = rioxarray.open_rasterio("satellite_image.tif")
print(raster)
If you have a MultiBand GeoTIFF, it can be opened by defining dimensions like this:
band
y
x
Spatial data and other metadata can also be retained.
You can find out what coordinate reference system you have with:
print(raster.rio.crs)
It is necessary, as GIS operations rely on it.
Selecting Raster Bands
Suppose a satellite image contains several spectral bands. Xarray makes it straightforward to select individual bands.
red = raster.sel(band=3)
nir = raster.sel(band=4)The exact band numbers depend on the satellite sensor and the organization of the raster file.
You can also select a range of coordinates:
subset = rastersel(
x=slice(500000, 510000),
y=slice(4500000, 4490000)
)This allows you to create spatial subsets without manually calculating array indexes.
Performing Raster Calculations
One of Xarray's major advantages is its ability to perform element-wise calculations while retaining coordinate information.
For example, the Normalized Difference Vegetation Index (NDVI) can be calculated as:
ndvi = (nir - red) / (nir + red)Because red and nir are Xarray objects, the calculation is performed across corresponding cells.
You can give the resulting layer a descriptive name:
ndvi.name = "NDVI"For production workflows, you should also handle division by zero and invalid pixels appropriately.
Working with NoData Values
Geospatial rasters are usually characterized by the presence of NoData pixels denoting regions that have no data at all.
It is critical to be aware of how NoData values are defined and handled when working with geospatial rasters.
One can use the capabilities of Xarray for dealing with missing data as follows:
valid_data = rasterwhere(raster != raster.rio.nodata)
The actual way to do it depends on the source data.
Handling Large Raster Datasets
Large satellite and scientific datasets can exceed available system memory. Xarray can work with Dask to process data in chunks.
For example:
import xarray as xr
dataset = xr.open_dataset(
"large_dataset.nc",
chunks={"time": 10, "y": 512, "x": 512}
)Instead of loading the entire dataset into memory immediately, the data can be processed in smaller chunks.
This is particularly useful for large Earth observation archives and time-series analysis.
Xarray and Cloud-Optimized Geospatial Workflows
Another application of Xarray can involve the use of current cloud-based data architectures. Data can be saved in a format like Zarr in chunks that can be efficiently processed.
A standard workflow could involve:
Cloud storage → Zarr/COG → Xarray → Dask → geospatial processing
The workflow architecture can be useful for scaling the processing of big remote sensing and environmental data.
As for the cloud-native GIS workflow, the correct chunk size is an important factor since improper chunking may lead to increased expenses.
Best Practices of Xarray Application in GIS
While employing Xarray for raster operations, follow the guidelines below:
Coordinate information should be maintained throughout the processing flow.
The CRS should be checked before conducting any spatial processing.
Dimensions and variables should have meaningful names.
Invalid values and NoData should be managed appropriately.
Use Dask for datasets that cannot fit into memory.
Select suitable resampling techniques during reprojections.
Preserve metadata during export of the processed datasets.
Use chunking wisely for multidimensional large datasets.
Raster alignment validation should be conducted prior to band math.
Python Xarray is an extensible data model that is ideal when dealing with multidimensional raster datasets. Xarray has features such as labeled dimensions, coordinate handling capabilities, metadata, and Dask support that make it suitable for satellite imagery and other types of Earth Observation datasets.
By incorporating rioxarray into Xarray, a practical workflow for geospatial processing of data that involves reading rasters, dealing with the Coordinate Reference System, clipping, reprojecting, doing raster algebra, masking, and exporting to GeoTIFF is possible. In the case of large datasets, Dask and Zarr can further be used to achieve scalability in geospatial processing.
GIS and remote sensing analysts handling big and multidimensional datasets can greatly benefit from Xarray.
To learn more about Xarray and its geospatial capabilities, click here.
For more information or any questions regarding NDVI, NDWI, NDBI, and SAVI, please don't hesitate to contact us at
Email: info@geowgs84.com
USA (HQ): (720) 702–4849
(A GeoWGS84 Corp Company)
Comments