# Xarray for multidimensional gridded data¶

In last week's lecture, we saw how Pandas provided a way to keep track of additional "metadata" surrounding tabular datasets, including "indexes" for each row and labels for each column. These features, together with Pandas' many useful routines for all kinds of data munging and analysis, have made Pandas one of the most popular python packages in the world.

However, not all Earth science datasets easily fit into the "tabular" model (i.e. rows and columns) imposed by Pandas. In particular, we often deal with multidimensional data. By multidimensional data (also often called N-dimensional), I mean data with many independent dimensions or axes. For example, we might represent Earth's surface temperature $T$ as a three dimensional variable

$$T(x, y, t)$$

where $x$ is longitude, $y$ is latitude, and $t$ is time.

The point of xarray is to provide pandas-level convenience for working with this type of data.

## Learning Goals for Xarray¶

Because of the importance of xarray for data analysis in geoscience, we are going to spend a long time on it. The goals of the next three lessons include the following.

### Lesson 1: Xarray Fundamentals¶

#### Dataset Creation¶

1. Describe the core xarray data structures, the DataArray and the Dataset, and the components that make them up, including: Data Variables, Dimensions, Coordinates, Indexes, and Attributes
2. Create xarray DataArrays and DataSets out of raw numpy arrays
3. Create xarray objects with and without indexes
4. Load xarray datasets from netCDF files and openDAP servers
5. View and set attributes

#### Basic Indexing and Interpolation¶

1. Select data by position using .isel with values or slices
2. Select data by label using .sel with values or slices
3. Select timeseries data by date/time with values or slices
4. Use nearest-neighbor lookups with .sel
5. Mask data with .where
6. Interpolate data in one and several dimensions

#### Basic Computation¶

1. Do basic arithmetic with DataArrays and Datasets
2. Use numpy universal function on DataArrays and Datasets, or use corresponding built-in xarray methods
3. Combine multiple xarray objects in arithmetic operations and understand how they are broadcasted / aligned
4. Perform aggregation (reduction) along one or multiple dimensions of a DataArray or Dataset

#### Basic Plotting¶

1. Use built-in xarray plotting for 1D and 2D DataArrays
2. Customize plots with options

#### Xarray's groupby, resample, and rolling¶

1. Split xarray objects into groups using groupby
2. Apply reduction operations to groups (e.g. mean)
3. Apply non-reducing functions to groups (e.g. standardize)
4. Use groupby with time coordinates (e.g. to create climatologies)
5. Use artimetic between GroupBy objects and regular DataArrays / Datasets
6. Use groupby_bins to aggregate data in bins
7. Use resample on time dimensions
8. Use rolling to apply rolling aggregations

#### Merging Combining Datasets¶

1. Concatentate DataArrays and Datasets along a new or existing dimension
2. Merge multiple datasets with different variables
3. Add a new data variable to an existing Dataset

#### Reshaping Data¶

1. Transpose dimension order
2. Swap coordinates
3. Expand and squeeze dimensions
4. Convert between DataArray and Dataset
5. Use stack and unstack to transform data

1. Use differentiate to take derivatives of data
2. Use apply_ufunc to apply custom or specialized operations to data

#### Plotting¶

1. Show multiple line plots over a dimension using the hue keyword
2. Create multiple 2D plots using faceting

# Lesson 1: Xarray Fundamentals¶

## Xarray data structures¶

Like Pandas, xarray has two fundamental data structures:

• a DataArray, which holds a single multi-dimensional variable and its coordinates
• a Dataset, which holds multiple variables that potentially share the same coordinates

### DataArray¶

A DataArray has four essential attributes:

• values: a numpy.ndarray holding the array’s values
• dims: dimension names for each axis (e.g., ('x', 'y', 'z'))
• coords: a dict-like container of arrays (coordinates) that label each point (e.g., 1-dimensional arrays of numbers, datetime objects or strings)
• attrs: an OrderedDict to hold arbitrary metadata (attributes)

Let's start by constructing some DataArrays manually

In [1]:
import numpy as np
import xarray as xr
from matplotlib import pyplot as plt
%matplotlib inline
plt.rcParams['figure.figsize'] = (8,5)


A simple DataArray without dimensions or coordinates isn't much use.

In [2]:
da = xr.DataArray([9, 0, 2, 1, 0])
da

Out[2]:

array([9, 0, 2, 1, 0])
Dimensions without coordinates: dim_0

We can add a dimension name...

In [3]:
da = xr.DataArray([9, 0, 2, 1, 0], dims=['x'])
da

Out[3]:

array([9, 0, 2, 1, 0])
Dimensions without coordinates: x

But things get most interesting when we add a coordinate:

In [4]:
da = xr.DataArray([9, 0, 2, 1, 0],
dims=['x'],
coords={'x': [10, 20, 30, 40, 50]})
da

Out[4]:

array([9, 0, 2, 1, 0])
Coordinates:
* x        (x) int64 10 20 30 40 50

Xarray has built-in plotting, like pandas.

In [5]:
da.plot(marker='o')

Out[5]:
[]

### Multidimensional DataArray¶

If we are just dealing with 1D data, Pandas and Xarray have very similar capabilities. Xarray's real potential comes with multidimensional data.

Let's go back to the multidimensional ARGO data we loaded in the numpy lesson. If you haven't already downloaded it, you can do so at the command line with

curl -O http://www.ldeo.columbia.edu/~rpa/argo_float_4901412.npz


We reload this data and examine its keys.

In [6]:
argo_data = np.load('argo_float_4901412.npz')
argo_data.keys()

Out[6]:
['S', 'T', 'levels', 'lon', 'date', 'P', 'lat']

The values of the argo_data object are numpy arrays.

In [7]:
S = argo_data.f.S
T = argo_data.f.T
P = argo_data.f.P
levels = argo_data.f.levels
lon = argo_data.f.lon
lat = argo_data.f.lat
date = argo_data.f.date
print(S.shape, lon.shape, date.shape)

(78, 75) (75,) (75,)


Let's organize the data and coordinates of the salinity variable into a DataArray.

In [8]:
da_salinity = xr.DataArray(S, dims=['level', 'date'],
coords={'level': levels,
'date': date},)
da_salinity

Out[8]:

array([[ 35.638939,  35.514957,  35.572971, ...,  35.820938,  35.777939,
35.668911],
[ 35.633938,  35.521957,  35.573971, ...,  35.810932,  35.583897,
35.667912],
[ 35.681946,  35.525959,  35.572971, ...,  35.795929,  35.662907,
35.665913],
...,
[ 34.915859,  34.923904,  34.923904, ...,  34.934811,  34.940811,
34.946808],
[ 34.915859,  34.923904,  34.921906, ...,  34.932808,  34.93681 ,
34.94381 ],
[ 34.917858,  34.923904,  34.923904, ...,        nan,  34.93681 ,
nan]])
Coordinates:
* level    (level) int64 0 1 2 3 4 5 6 7 8 9 ... 68 69 70 71 72 73 74 75 76 77
* date     (date) datetime64[ns] 2012-07-13T22:33:06.019200 ... 2014-07-24T03:02:33.014400
In [9]:
da_salinity.plot(yincrease=False)

Out[9]:

Attributes can be used to store metadata. What metadata should you store? The CF Conventions are a great resource for thinking about climate metadata. Below we define two of the required CF-conventions attributes.

In [10]:
da_salinity.attrs['units'] = 'PSU'
da_salinity.attrs['standard_name'] = 'sea_water_salinity'
da_salinity

Out[10]:

array([[ 35.638939,  35.514957,  35.572971, ...,  35.820938,  35.777939,
35.668911],
[ 35.633938,  35.521957,  35.573971, ...,  35.810932,  35.583897,
35.667912],
[ 35.681946,  35.525959,  35.572971, ...,  35.795929,  35.662907,
35.665913],
...,
[ 34.915859,  34.923904,  34.923904, ...,  34.934811,  34.940811,
34.946808],
[ 34.915859,  34.923904,  34.921906, ...,  34.932808,  34.93681 ,
34.94381 ],
[ 34.917858,  34.923904,  34.923904, ...,        nan,  34.93681 ,
nan]])
Coordinates:
* level    (level) int64 0 1 2 3 4 5 6 7 8 9 ... 68 69 70 71 72 73 74 75 76 77
* date     (date) datetime64[ns] 2012-07-13T22:33:06.019200 ... 2014-07-24T03:02:33.014400
Attributes:
units:          PSU
standard_name:  sea_water_salinity

### Datasets¶

A Dataset holds many DataArrays which potentially can share coordinates. In analogy to pandas:

pandas.Series : pandas.Dataframe :: xarray.DataArray : xarray.Dataset



Constructing Datasets manually is a bit more involved in terms of syntax. The Dataset constructor takes three arguments:

• data_vars should be a dictionary with each key as the name of the variable and each value as one of:
• A DataArray or Variable
• A tuple of the form (dims, data[, attrs]), which is converted into arguments for Variable
• A pandas object, which is converted into a DataArray
• A 1D array or list, which is interpreted as values for a one dimensional coordinate variable along the same dimension as it’s name
• coords should be a dictionary of the same form as data_vars.
• attrs should be a dictionary.

Let's put together a Dataset with temperature, salinity and pressure all together

In [11]:
argo = xr.Dataset(
data_vars={'salinity':    (('level', 'date'), S),
'temperature': (('level', 'date'), T),
'pressure':    (('level', 'date'), P)},
coords={'level': levels,
'date': date})
argo

Out[11]:

Dimensions:      (date: 75, level: 78)
Coordinates:
* level        (level) int64 0 1 2 3 4 5 6 7 8 ... 69 70 71 72 73 74 75 76 77
* date         (date) datetime64[ns] 2012-07-13T22:33:06.019200 ... 2014-07-24T03:02:33.014400
Data variables:
salinity     (level, date) float64 35.64 35.51 35.57 35.4 ... nan 34.94 nan
temperature  (level, date) float64 18.97 18.44 19.1 19.79 ... nan 3.714 nan
pressure     (level, date) float64 6.8 6.1 6.5 5.0 ... 2e+03 nan 2e+03 nan

What about lon and lat? We forgot them in the creation process, but we can add them after the fact.

In [12]:
argo.coords['lon'] = lon
argo

Out[12]:

Dimensions:      (date: 75, level: 78, lon: 75)
Coordinates:
* level        (level) int64 0 1 2 3 4 5 6 7 8 ... 69 70 71 72 73 74 75 76 77
* date         (date) datetime64[ns] 2012-07-13T22:33:06.019200 ... 2014-07-24T03:02:33.014400
* lon          (lon) float64 -39.13 -37.28 -36.9 ... -33.83 -34.11 -34.38
Data variables:
salinity     (level, date) float64 35.64 35.51 35.57 35.4 ... nan 34.94 nan
temperature  (level, date) float64 18.97 18.44 19.1 19.79 ... nan 3.714 nan
pressure     (level, date) float64 6.8 6.1 6.5 5.0 ... 2e+03 nan 2e+03 nan

That was not quite right...we want lon to have dimension date:

In [13]:
del argo['lon']
argo.coords['lon'] = ('date', lon)
argo.coords['lat'] = ('date', lat)
argo

Out[13]:

Dimensions:      (date: 75, level: 78)
Coordinates:
* level        (level) int64 0 1 2 3 4 5 6 7 8 ... 69 70 71 72 73 74 75 76 77
* date         (date) datetime64[ns] 2012-07-13T22:33:06.019200 ... 2014-07-24T03:02:33.014400
lon          (date) float64 -39.13 -37.28 -36.9 ... -33.83 -34.11 -34.38
lat          (date) float64 47.19 46.72 46.45 46.23 ... 42.6 42.46 42.38
Data variables:
salinity     (level, date) float64 35.64 35.51 35.57 35.4 ... nan 34.94 nan
temperature  (level, date) float64 18.97 18.44 19.1 19.79 ... nan 3.714 nan
pressure     (level, date) float64 6.8 6.1 6.5 5.0 ... 2e+03 nan 2e+03 nan

### Coordinates vs. Data Variables¶

Data variables can be modified through arithmentic operations or other functions. Coordinates are always keept the same.

In [14]:
argo * 10000

Out[14]:

Dimensions:      (date: 75, level: 78)
Coordinates:
* level        (level) int64 0 1 2 3 4 5 6 7 8 ... 69 70 71 72 73 74 75 76 77
* date         (date) datetime64[ns] 2012-07-13T22:33:06.019200 ... 2014-07-24T03:02:33.014400
lat          (date) float64 47.19 46.72 46.45 46.23 ... 42.6 42.46 42.38
lon          (date) float64 -39.13 -37.28 -36.9 ... -33.83 -34.11 -34.38
Data variables:
salinity     (level, date) float64 3.564e+05 3.551e+05 ... 3.494e+05 nan
temperature  (level, date) float64 1.897e+05 1.844e+05 ... 3.714e+04 nan
pressure     (level, date) float64 6.8e+04 6.1e+04 6.5e+04 ... nan 2e+07 nan

Clearly lon and lat are coordinates rather than data variables. We can change their status as follows:

In [15]:
argo = argo.set_coords(['lon', 'lat'])
argo

Out[15]:

Dimensions:      (date: 75, level: 78)
Coordinates:
* level        (level) int64 0 1 2 3 4 5 6 7 8 ... 69 70 71 72 73 74 75 76 77
* date         (date) datetime64[ns] 2012-07-13T22:33:06.019200 ... 2014-07-24T03:02:33.014400
lon          (date) float64 -39.13 -37.28 -36.9 ... -33.83 -34.11 -34.38
lat          (date) float64 47.19 46.72 46.45 46.23 ... 42.6 42.46 42.38
Data variables:
salinity     (level, date) float64 35.64 35.51 35.57 35.4 ... nan 34.94 nan
temperature  (level, date) float64 18.97 18.44 19.1 19.79 ... nan 3.714 nan
pressure     (level, date) float64 6.8 6.1 6.5 5.0 ... 2e+03 nan 2e+03 nan

The * symbol in the representation above indicates that level and date are "dimension coordinates" (they describe the coordinates associated with data variable axes) while lon and lat are "non-dimension coordinates". We can make any variable a non-dimension coordiante.

Alternatively, we could have assigned directly to coords as follows:

In [16]:
argo.coords['lon'] = ('date', lon)
argo.coords['lat'] = ('date', lat)


## Working with Labeled Data¶

Xarray's labels make working with multidimensional data much easier.

### Selecting Data (Indexing)¶

We can always use regular numpy indexing and slicing on DataArrays

In [17]:
argo.salinity[2].plot()

Out[17]:
[]
In [18]:
argo.salinity[:, 10].plot()

Out[18]:
[]

However, it is often much more powerful to use xarray's .sel() method to use label-based indexing.

In [19]:
argo.salinity.sel(level=2).plot()

Out[19]:
[]
In [20]:
argo.salinity.sel(date='2012-10-22').plot(y='level', yincrease=False)

Out[20]:
[]

.sel() also supports slicing. Unfortunately we have to use a somewhat awkward syntax, but it still works.

In [21]:
argo.salinity.sel(date=slice('2012-10-01', '2012-12-01')).plot()

Out[21]:

.sel() also works on the whole Dataset

In [22]:
argo.sel(date='2012-10-22')

Out[22]:

Dimensions:      (date: 1, level: 78)
Coordinates:
* level        (level) int64 0 1 2 3 4 5 6 7 8 ... 69 70 71 72 73 74 75 76 77
* date         (date) datetime64[ns] 2012-10-22T02:50:32.006400
lon          (date) float64 -32.97
lat          (date) float64 44.13
Data variables:
salinity     (level, date) float64 35.47 35.47 35.47 ... 34.93 34.92 nan
temperature  (level, date) float64 17.13 17.13 17.13 ... 3.736 3.639 nan
pressure     (level, date) float64 6.4 10.3 15.4 ... 1.9e+03 1.951e+03 nan

### Computation¶

Xarray dataarrays and datasets work seamlessly with arithmetic operators and numpy array functions.

In [23]:
temp_kelvin = argo.temperature + 273.15
temp_kelvin.plot(yincrease=False)

Out[23]:

We can also combine multiple xarray datasets in arithemtic operations

In [24]:
g = 9.8
buoyancy = g * (2e-4 * argo.temperature - 7e-4 * argo.salinity)
buoyancy.plot(yincrease=False)

Out[24]:

Broadcasting arrays in numpy is a nightmare. It is much easier when the data axes are labeled!

This is a useless calculation, but it illustrates how perfoming an operation on arrays with differenty coordinates will result in automatic broadcasting

In [25]:
level_times_lat = argo.level * argo.lat
level_times_lat

Out[25]:

array([[    0.   ,     0.   ,     0.   , ...,     0.   ,     0.   ,     0.   ],
[   47.187,    46.716,    46.45 , ...,    42.601,    42.457,    42.379],
[   94.374,    93.432,    92.9  , ...,    85.202,    84.914,    84.758],
...,
[ 3539.025,  3503.7  ,  3483.75 , ...,  3195.075,  3184.275,  3178.425],
[ 3586.212,  3550.416,  3530.2  , ...,  3237.676,  3226.732,  3220.804],
[ 3633.399,  3597.132,  3576.65 , ...,  3280.277,  3269.189,  3263.183]])
Coordinates:
* level    (level) int64 0 1 2 3 4 5 6 7 8 9 ... 68 69 70 71 72 73 74 75 76 77
* date     (date) datetime64[ns] 2012-07-13T22:33:06.019200 ... 2014-07-24T03:02:33.014400
lon      (date) float64 -39.13 -37.28 -36.9 -36.89 ... -33.83 -34.11 -34.38
lat      (date) float64 47.19 46.72 46.45 46.23 ... 42.72 42.6 42.46 42.38
In [26]:
level_times_lat.plot()

Out[26]:

### Reductions¶

Just like in numpy, we can reduce xarray DataArrays along any number of axes:

In [27]:
argo.temperature.mean(axis=0).dims

Out[27]:
('date',)
In [28]:
argo.temperature.mean(axis=1).dims

Out[28]:
('level',)

However, rather than performing reductions on axes (as in numpy), we can perform them on dimensions. This turns out to be a huge convenience

In [29]:
argo_mean = argo.mean(dim='date')
argo_mean

Out[29]:

Dimensions:      (level: 78)
Coordinates:
* level        (level) int64 0 1 2 3 4 5 6 7 8 ... 69 70 71 72 73 74 75 76 77
Data variables:
salinity     (level) float64 35.91 35.9 35.9 35.9 ... 34.94 34.94 34.93
temperature  (level) float64 17.6 17.57 17.51 17.42 ... 3.789 3.73 3.662
pressure     (level) float64 6.435 10.57 15.54 ... 1.95e+03 1.999e+03
In [30]:
argo_mean.salinity.plot(y='level', yincrease=False)

Out[30]:
[]
In [31]:
argo_std = argo.std(dim='date')
argo_std.salinity.plot(y='level', yincrease=False)

Out[31]:
[]

NetCDF (Network Common Data Format) is the most widely used format for distributing geoscience data. NetCDF is maintained by the Unidata organization.

Below we quote from the NetCDF website:

NetCDF (network Common Data Form) is a set of interfaces for array-oriented data access and a freely distributed collection of data access libraries for C, Fortran, C++, Java, and other languages. The netCDF libraries support a machine-independent format for representing scientific data. Together, the interfaces, libraries, and format support the creation, access, and sharing of scientific data.

NetCDF data is:

• Self-Describing. A netCDF file includes information about the data it contains.
• Portable. A netCDF file can be accessed by computers with different ways of storing integers, characters, and floating-point numbers.
• Scalable. A small subset of a large dataset may be accessed efficiently.
• Appendable. Data may be appended to a properly structured netCDF file without copying the dataset or redefining its structure.
• Sharable. One writer and multiple readers may simultaneously access the same netCDF file.
• Archivable. Access to all earlier forms of netCDF data will be supported by current and future versions of the software.

Xarray was designed to make reading netCDF files in python as easy, powerful, and flexible as possible. (See xarray netCDF docs for more details.)

In [32]:
! wget https://data.giss.nasa.gov/pub/gistemp/gistemp1200_ERSSTv5.nc.gz

In [33]:
ds = xr.open_dataset('gistemp1200_ERSSTv5.nc')
ds

Out[33]:

Dimensions:      (lat: 90, lon: 180, nv: 2, time: 1665)
Coordinates:
* lat          (lat) float32 -89.0 -87.0 -85.0 -83.0 ... 83.0 85.0 87.0 89.0
* lon          (lon) float32 -179.0 -177.0 -175.0 -173.0 ... 175.0 177.0 179.0
* time         (time) datetime64[ns] 1880-01-15 1880-02-15 ... 2018-09-15
Dimensions without coordinates: nv
Data variables:
time_bnds    (time, nv) int32 ...
tempanomaly  (time, lat, lon) float32 ...
Attributes:
title:        GISTEMP Surface Temperature Analysis
institution:  NASA Goddard Institute for Space Studies
source:       http://data.giss.nasa.gov/gistemp/
Conventions:  CF-1.6
history:      Created 2018-10-12 16:36:37 by SBBX_to_nc 2.0 - ILAND=1200,...
In [34]:
ds.tempanomaly.isel(time=-1).plot()

Out[34]:
In [35]:
ds.tempanomaly.mean(dim=('lon', 'lat')).plot()

Out[35]:
[]
In [ ]: