# Introduction to xarray¶

• Prerequisites: Users of this notebook should have a basic understanding of:

**Keywords** :index:beginner's guide; xarray, :index:python package; xarray

## Background¶

Xarray is a python library which simplifies working with labelled multi-dimension arrays. Xarray introduces labels in the forms of dimensions, coordinates and attributes on top of raw numpy arrays, allowing for more intitutive and concise development. More information about xarray data structures and functions can be found here.

## Description¶

This notebook is designed to introduce users to xarray using Python code in Jupyter Notebooks via JupyterLab.

Topics covered include:

• How to use xarray functions in a Jupyter Notebook cell

• How to access xarray dimensions and metadata

• Using indexing to explore multi-dimensional xarray data

• Appliction of built-in xarray functions such as sum, std and mean

## Getting started¶

To run this notebook, run all the cells in the notebook starting with the “Load packages” cell. For help with running notebook cells, refer back to the Jupyter Notebooks notebook.

[1]:

%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
import xarray as xr


### Introduction to xarray¶

DE Africa uses xarray as its core data model. To better understand what it is, let first do a simple experiment using a combination of plain numpy arrays and Python dictionaries.

Suppose we have a satellite image with three bands: Red, NIR and SWIR. These bands are represented as 2-dimensional numpy arrays and the latitude and longitude coordinates for each dimension are represented using 1-dimensional arrays. Finally, we also have some metadata that comes with this image. The code below creates fake satellite data and structures the data as a dictionary.

[2]:

#create fake satellite data
red = np.random.rand(250,250)
nir = np.random.rand(250,250)
swir = np.random.rand(250,250)

#create some lats and lons
lats = np.linspace(-23.5, -26.0, num=red.shape[0], endpoint=False)
lons = np.linspace(110.0, 112.5, num=red.shape[1], endpoint=False)

title = "Image of the desert"
date = "2019-11-10"

#stack into a dictionary
image = {"red": red,
"nir": nir,
"swir": swir,
"latitude": lats,
"longitude": lons,
"title": title,
"date": date}


All our data is conveniently packed in a dictionary. Now we can use this dictionary to work with the data it contains:

[3]:

#date of satellite image
print(image["date"])

#mean of red values
image["red"].mean()

2019-11-10

[3]:

0.49993184086419146


Still, to select data we have to use numpy indexes. Wouldn’t it be convenient to be able to select data from the images using the coordinates of the pixels instead of their relative positions? This is exactly what xarray solves! Let’s see how it works:

To explore xarray we have a file containing some surface reflectance data extracted from the DE Africa platform. The object that we get ds is a xarray Dataset, which in some ways is very similar to the dictionary that we created before, but with lots of convenient functionality available.

[4]:

ds = xr.open_dataset('../Supplementary_data/07_Intro_to_xarray/example_netcdf.nc')
ds

[4]:

<xarray.Dataset>
Dimensions:      (time: 12, y: 601, x: 483)
Coordinates:
* time         (time) datetime64[ns] 2018-01-03T08:31:05 ... 2018-02-27T08:...
* y            (y) float64 -2.519e+06 -2.519e+06 ... -2.507e+06 -2.507e+06
* x            (x) float64 2.378e+06 2.378e+06 ... 2.388e+06 2.388e+06
spatial_ref  int32 ...
Data variables:
red          (time, y, x) uint16 ...
green        (time, y, x) uint16 ...
blue         (time, y, x) uint16 ...
Attributes:
crs:           EPSG:6933
grid_mapping:  spatial_ref

### Xarray dataset structure¶

A Dataset can be seen as a dictionary structure packing up the data, dimensions and attributes. Variables in a Dataset object are called DataArrays and they share dimensions with the higher level Dataset. The figure below provides an illustrative example:

To access a variable we can access as if it were a Python dictionary, or using the . notation, which is more convenient.

[5]:

ds["green"]

#or alternatively

ds.green

[5]:

<xarray.DataArray 'green' (time: 12, y: 601, x: 483)>
[3483396 values with dtype=uint16]
Coordinates:
* time         (time) datetime64[ns] 2018-01-03T08:31:05 ... 2018-02-27T08:...
* y            (y) float64 -2.519e+06 -2.519e+06 ... -2.507e+06 -2.507e+06
* x            (x) float64 2.378e+06 2.378e+06 ... 2.388e+06 2.388e+06
spatial_ref  int32 ...
Attributes:
units:         1
nodata:        0
crs:           EPSG:6933
grid_mapping:  spatial_ref

Dimensions are also stored as numeric arrays that we can easily access

[6]:

ds['time']

#or alternatively

ds.time

[6]:

<xarray.DataArray 'time' (time: 12)>
array(['2018-01-03T08:31:05.000000000', '2018-01-08T08:34:01.000000000',
'2018-01-13T08:30:41.000000000', '2018-01-18T08:30:42.000000000',
'2018-01-23T08:33:58.000000000', '2018-01-28T08:30:20.000000000',
'2018-02-07T08:30:53.000000000', '2018-02-12T08:31:43.000000000',
'2018-02-17T08:23:09.000000000', '2018-02-17T08:35:40.000000000',
'2018-02-22T08:34:52.000000000', '2018-02-27T08:31:36.000000000'],
dtype='datetime64[ns]')
Coordinates:
* time         (time) datetime64[ns] 2018-01-03T08:31:05 ... 2018-02-27T08:...
spatial_ref  int32 ...

Metadata is referred to as attributes and is internally stored under .attrs, but the same convenient . notation applies to them.

[7]:

ds.attrs['crs']

#or alternatively

ds.crs

[7]:

'EPSG:6933'


DataArrays store their data internally as multidimensional numpy arrays. But these arrays contain dimensions or labels that make it easier to handle the data. To access the underlaying numpy array of a DataArray we can use the .values notation.

[8]:

arr = ds.green.values

type(arr), arr.shape

[8]:

(numpy.ndarray, (12, 601, 483))


### Indexing¶

Xarray offers two different ways of selecting data. This includes the isel() approach, where data can be selected based on its index (like numpy).

[9]:

print(ds.time.values)

ss = ds.green.isel(time=0)
ss

['2018-01-03T08:31:05.000000000' '2018-01-08T08:34:01.000000000'
'2018-01-13T08:30:41.000000000' '2018-01-18T08:30:42.000000000'
'2018-01-23T08:33:58.000000000' '2018-01-28T08:30:20.000000000'
'2018-02-07T08:30:53.000000000' '2018-02-12T08:31:43.000000000'
'2018-02-17T08:23:09.000000000' '2018-02-17T08:35:40.000000000'
'2018-02-22T08:34:52.000000000' '2018-02-27T08:31:36.000000000']

[9]:

<xarray.DataArray 'green' (y: 601, x: 483)>
array([[1214, 1232, 1406, ..., 3436, 4252, 4300],
[1214, 1334, 1378, ..., 2006, 2602, 4184],
[1274, 1340, 1554, ..., 2436, 1858, 1890],
...,
[1142, 1086, 1202, ..., 1096, 1074, 1092],
[1188, 1258, 1190, ..., 1058, 1138, 1138],
[1152, 1134, 1074, ..., 1086, 1116, 1100]], dtype=uint16)
Coordinates:
time         datetime64[ns] 2018-01-03T08:31:05
* y            (y) float64 -2.519e+06 -2.519e+06 ... -2.507e+06 -2.507e+06
* x            (x) float64 2.378e+06 2.378e+06 ... 2.388e+06 2.388e+06
spatial_ref  int32 ...
Attributes:
units:         1
nodata:        0
crs:           EPSG:6933
grid_mapping:  spatial_ref

Or the sel() approach, used for selecting data based on its dimension of label value.

[10]:

ss = ds.green.sel(time='2018-01-08')
ss

[10]:

<xarray.DataArray 'green' (time: 1, y: 601, x: 483)>
array([[[1270, 1280, ..., 4228, 3950],
[1266, 1332, ..., 3880, 4372],
...,
[1172, 1180, ..., 1154, 1190],
[1242, 1204, ..., 1192, 1170]]], dtype=uint16)
Coordinates:
* time         (time) datetime64[ns] 2018-01-08T08:34:01
* y            (y) float64 -2.519e+06 -2.519e+06 ... -2.507e+06 -2.507e+06
* x            (x) float64 2.378e+06 2.378e+06 ... 2.388e+06 2.388e+06
spatial_ref  int32 ...
Attributes:
units:         1
nodata:        0
crs:           EPSG:6933
grid_mapping:  spatial_ref

Slicing data is also used to select a subset of data.

[11]:

ss.x.values[100]

[11]:

2380390.0

[12]:

ss = ds.green.sel(time='2018-01-08', x=slice(2378390,2380390))
ss

[12]:

<xarray.DataArray 'green' (time: 1, y: 601, x: 101)>
array([[[1270, 1280, ..., 1416, 1290],
[1266, 1332, ..., 1368, 1274],
...,
[1172, 1180, ..., 1086,  991],
[1242, 1204, ..., 1019,  986]]], dtype=uint16)
Coordinates:
* time         (time) datetime64[ns] 2018-01-08T08:34:01
* y            (y) float64 -2.519e+06 -2.519e+06 ... -2.507e+06 -2.507e+06
* x            (x) float64 2.378e+06 2.378e+06 2.378e+06 ... 2.38e+06 2.38e+06
spatial_ref  int32 ...
Attributes:
units:         1
nodata:        0
crs:           EPSG:6933
grid_mapping:  spatial_ref

Xarray exposes lots of functions to easily transform and analyse Datasets and DataArrays. For example, to calculate the spatial mean, standard deviation or sum of the green band:

[13]:

print("Mean of green band:", ds.green.mean().values)
print("Standard deviation of green band:", ds.green.std().values)
print("Sum of green band:", ds.green.sum().values)

Mean of green band: 4141.488778766468
Standard deviation of green band: 3775.5536474649584
Sum of green band: 14426445446


### Plotting data with Matplotlib¶

Plotting is also conveniently integrated in the library.

[14]:

ds["green"].isel(time=0).plot()

[14]:

<matplotlib.collections.QuadMesh at 0x7f88d342bd30>


but we still can do things manually using numpy and matplotlib if you choose:

[15]:

rgb = np.dstack((ds.red.isel(time=0).values, ds.green.isel(time=0).values, ds.blue.isel(time=0).values))
rgb = np.clip(rgb, 0, 2000) / 2000

plt.imshow(rgb);


But compare the above to elegantly chaining operations within xarray:

[16]:

ds[['red', 'green', 'blue']].isel(time=0).to_array().plot.imshow(robust=True, figsize=(6, 6));


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.

Last Tested:

[17]:

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

[17]:

'2023-08-11'