Lecture 8: 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 and environmental 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.

Xarray Data Model¶

Here is an example of how we might structure a dataset for a weather forecast:

No description has been provided for this image

You'll notice multiple data variables (temperature, precipitation), coordinate variables (latitude, longitude), and dimensions (x, y, t). We'll cover how these fit into Xarray's data structures below.

Xarray functionality¶

Xarray doesn’t just keep track of labels on arrays – it uses them to provide a powerful and concise interface. For example:

Apply operations over dimensions by name: x.sum('time').
Select values by label (or logical location) instead of integer location: x.loc['2014-01-01'] or x.sel(time='2014-01-01').
Mathematical operations (e.g., x - y) vectorize across multiple dimensions (array broadcasting) based on dimension names, not shape.
Easily use the split-apply-combine paradigm with groupby: x.groupby('time.dayofyear').mean().
Database-like alignment based on coordinate labels that smoothly handles missing values: x, y = xr.align(x, y, join='outer').
Keep track of arbitrary metadata in the form of a Python dictionary: x.attrs.

The N-dimensional nature of xarray’s data structures makes it suitable for dealing with multi-dimensional scientific data, and its use of dimension names instead of axis labels (dim='time' instead of axis=0) makes such arrays much more manageable than the raw numpy ndarray: with xarray, you don’t need to keep track of the order of an array’s dimensions or insert dummy dimensions of size 1 to align arrays (e.g., using np.newaxis).

The immediate payoff of using xarray is that you’ll write less code. The long-term payoff is that you’ll understand what you were thinking when you come back to look at it weeks or months later.

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

Install Xarray:¶

$ conda install -c conda-forge xarray dask netCDF4 bottleneck python-graphviz

In [1]:

Copied!





# First import xarray 
import matplotlib.pyplot as plt
import numpy as np
import xarray as xr
# First import xarray 
import matplotlib.pyplot as plt
import numpy as np
import xarray as xr

Xarray Data Structure¶

Xarray has a few small real-world tutorial datasets hosted in this GitHub repository https://github.com/pydata/xarray-data. We'll use the xarray.tutorial.load_dataset convenience function to download and open the air_temperature (National Centers for Environmental Prediction) Dataset by name.

We can access "layers" of the Dataset (individual DataArrays) with dictionary syntax

We can save some typing by using the "attribute" or "dot" notation. This won't work for variable names that clash with built-in method names (for example, mean).

What is included in xarray `Dataset`?¶

The output consists of:

a summary of all dimensions of the Dataset (lat: 25, time: 2920, lon: 53): this tells us that the first dimension is named lat and has a size of 25, the second dimension is named time and has a size of 2920, and the third dimension is named lon and has a size of 53. Because we will access the dimensions by name, the order doesn't matter.
an unordered list of coordinates or dimensions with coordinates with one item per line. Each item has a name, one or more dimensions in parentheses, a dtype and a preview of the values.
an alphabetically sorted list of dimensions without coordinates (if there are any)
an unordered list of attributes, or metadata

DataArray¶

The DataArray class consists of an array (data) and its associated dimension names, labels, and attributes (metadata).

We can also access the data array directly:

In [6]:

Copied!

ds.air.data
ds.air.data

Out[6]:

array([[[241.2    , 242.5    , 243.5    , ..., 232.79999, 235.5    ,
         238.59999],
        [243.79999, 244.5    , 244.7    , ..., 232.79999, 235.29999,
         239.29999],
        [250.     , 249.79999, 248.89   , ..., 233.2    , 236.39   ,
         241.7    ],
        ...,
        [296.6    , 296.19998, 296.4    , ..., 295.4    , 295.1    ,
         294.69998],
        [295.9    , 296.19998, 296.79   , ..., 295.9    , 295.9    ,
         295.19998],
        [296.29   , 296.79   , 297.1    , ..., 296.9    , 296.79   ,
         296.6    ]],

       [[242.09999, 242.7    , 243.09999, ..., 232.     , 233.59999,
         235.79999],
        [243.59999, 244.09999, 244.2    , ..., 231.     , 232.5    ,
         235.7    ],
        [253.2    , 252.89   , 252.09999, ..., 230.79999, 233.39   ,
         238.5    ],
        ...,
        [296.4    , 295.9    , 296.19998, ..., 295.4    , 295.1    ,
         294.79   ],
        [296.19998, 296.69998, 296.79   , ..., 295.6    , 295.5    ,
         295.1    ],
        [296.29   , 297.19998, 297.4    , ..., 296.4    , 296.4    ,
         296.6    ]],

       [[242.29999, 242.2    , 242.29999, ..., 234.29999, 236.09999,
         238.7    ],
        [244.59999, 244.39   , 244.     , ..., 230.29999, 232.     ,
         235.7    ],
        [256.19998, 255.5    , 254.2    , ..., 231.2    , 233.2    ,
         238.2    ],
        ...,
        [295.6    , 295.4    , 295.4    , ..., 296.29   , 295.29   ,
         295.     ],
        [296.19998, 296.5    , 296.29   , ..., 296.4    , 296.     ,
         295.6    ],
        [296.4    , 296.29   , 296.4    , ..., 297.     , 297.     ,
         296.79   ]],

       ...,

       [[243.48999, 242.98999, 242.09   , ..., 244.18999, 244.48999,
         244.89   ],
        [249.09   , 248.98999, 248.59   , ..., 240.59   , 241.29   ,
         242.68999],
        [262.69   , 262.19   , 261.69   , ..., 239.39   , 241.68999,
         245.18999],
        ...,
        [294.79   , 295.29   , 297.49   , ..., 295.49   , 295.38998,
         294.69   ],
        [296.79   , 297.88998, 298.29   , ..., 295.49   , 295.49   ,
         294.79   ],
        [298.19   , 299.19   , 298.79   , ..., 296.09   , 295.79   ,
         295.79   ]],

       [[245.79   , 244.79   , 243.48999, ..., 243.29   , 243.98999,
         244.79   ],
        [249.89   , 249.29   , 248.48999, ..., 241.29   , 242.48999,
         244.29   ],
        [262.38998, 261.79   , 261.29   , ..., 240.48999, 243.09   ,
         246.89   ],
        ...,
        [293.69   , 293.88998, 295.38998, ..., 295.09   , 294.69   ,
         294.29   ],
        [296.29   , 297.19   , 297.59   , ..., 295.29   , 295.09   ,
         294.38998],
        [297.79   , 298.38998, 298.49   , ..., 295.69   , 295.49   ,
         295.19   ]],

       [[245.09   , 244.29   , 243.29   , ..., 241.68999, 241.48999,
         241.79   ],
        [249.89   , 249.29   , 248.39   , ..., 239.59   , 240.29   ,
         241.68999],
        [262.99   , 262.19   , 261.38998, ..., 239.89   , 242.59   ,
         246.29   ],
        ...,
        [293.79   , 293.69   , 295.09   , ..., 295.29   , 295.09   ,
         294.69   ],
        [296.09   , 296.88998, 297.19   , ..., 295.69   , 295.69   ,
         295.19   ],
        [297.69   , 298.09   , 298.09   , ..., 296.49   , 296.19   ,
         295.69   ]]], dtype=float32)

Named dimensions¶

.dims are the named axes of your data. They may (dimension coordinates) or may not (dimensions without coordinates) have associated values.

In this case we have 2 spatial dimensions (latitude and longitude are stored with shorthand names lat and lon) and one temporal dimension (time).

In [7]:

Copied!

ds.air.dims
ds.air.dims

Out[7]:

('time', 'lat', 'lon')

Coordinates¶

.coords is a simple dict-like data container for mapping coordinate names to values.

Here we see the actual timestamps and spatial positions of our air temperature data:

In [8]:

Copied!

ds.air.coords
ds.air.coords

Out[8]:

Coordinates:
  * lat      (lat) float32 75.0 72.5 70.0 67.5 65.0 ... 25.0 22.5 20.0 17.5 15.0
  * lon      (lon) float32 200.0 202.5 205.0 207.5 ... 322.5 325.0 327.5 330.0
  * time     (time) datetime64[ns] 2013-01-01 ... 2014-12-31T18:00:00

Attributes¶

.attrs is a dictionary that can contain arbitrary Python objects (strings, lists, integers, dictionaries, etc.) containing information about your data. Your only limitation is that some attributes may not be writeable to certain file formats.

In [9]:

Copied!

ds.air.attrs
ds.air.attrs

Out[9]:

{'long_name': '4xDaily Air temperature at sigma level 995',
 'units': 'degK',
 'precision': 2,
 'GRIB_id': 11,
 'GRIB_name': 'TMP',
 'var_desc': 'Air temperature',
 'dataset': 'NMC Reanalysis',
 'level_desc': 'Surface',
 'statistic': 'Individual Obs',
 'parent_stat': 'Other',
 'actual_range': array([185.16, 322.1 ], dtype='>f4')}

In [10]:

Copied!

# assign your own attributes!
ds.air.attrs["new_attr"] = "xarray"
ds.air.attrs
# assign your own attributes!
ds.air.attrs["new_attr"] = "xarray"
ds.air.attrs

Out[10]:

{'long_name': '4xDaily Air temperature at sigma level 995',
 'units': 'degK',
 'precision': 2,
 'GRIB_id': 11,
 'GRIB_name': 'TMP',
 'var_desc': 'Air temperature',
 'dataset': 'NMC Reanalysis',
 'level_desc': 'Surface',
 'statistic': 'Individual Obs',
 'parent_stat': 'Other',
 'actual_range': array([185.16, 322.1 ], dtype='>f4'),
 'new_attr': 'xarray'}

Working with Labelled Data¶

Xarray's labels make working with multidimensional data much easier. Metadata provides context and provides code that is more legible. This reduces the likelihood of errors from typos and makes analysis more intuitive and fun!

In [11]:

Copied!

## Without Xarray, let's try to use matplotlib to plot the data

temp = ds.air.data
plt.pcolormesh(temp[0,:,:])
## Without Xarray, let's try to use matplotlib to plot the data

temp = ds.air.data
plt.pcolormesh(temp[0,:,:])

Out[11]:

<matplotlib.collections.QuadMesh at 0x1945afee0>

In [12]:

Copied!





### Add Lat and Lon

temp = ds.air.data
lat = ds.air.lat.data
lon = ds.air.lon.data
plt.pcolormesh(lon, lat, temp[0,:,:])
### Add Lat and Lon

temp = ds.air.data
lat = ds.air.lat.data
lon = ds.air.lon.data
plt.pcolormesh(lon, lat, temp[0,:,:])

Out[12]:

<matplotlib.collections.QuadMesh at 0x1947868b0>

In [13]:

Copied!

temp.mean(axis=1)  ## what did I just do? I can't tell by looking at this line.
temp.mean(axis=1)  ## what did I just do? I can't tell by looking at this line.

Out[13]:

array([[279.39798, 279.6664 , 279.66122, ..., 279.9508 , 280.31522,
        280.6624 ],
       [279.05722, 279.538  , 279.7296 , ..., 279.77563, 280.27002,
        280.79764],
       [279.0104 , 279.2808 , 279.5508 , ..., 279.682  , 280.19763,
        280.81403],
       ...,
       [279.63   , 279.934  , 280.534  , ..., 279.802  , 280.346  ,
        280.77798],
       [279.398  , 279.66602, 280.31796, ..., 279.766  , 280.34198,
        280.834  ],
       [279.27   , 279.354  , 279.88202, ..., 279.42596, 279.96997,
        280.48196]], dtype=float32)

In [14]:

Copied!

ds.air.isel(time=0).plot()
ds.air.isel(time=0).plot()

Out[14]:

<matplotlib.collections.QuadMesh at 0x194849940>

In [15]:

Copied!

ds.air.mean(dim="time").plot()
ds.air.mean(dim="time").plot()

Out[15]:

<matplotlib.collections.QuadMesh at 0x194942490>

Selecting Data (Indexing)¶

We can always use regular numpy indexing and slicing on DataArrays

Position-based Indexing¶

Indexing a DataArray directly works (mostly) just like it does for numpy ndarrays, except that the returned object is always another DataArray:

This approach however does not take advantage of the dimension names and coordinate location information that is present in a Xarray object.

Positional Indexing Using Dimension Names¶

Remembering the axis order can be challenging even with 2D arrays:

is np_array[0,3] the first row and third column or first column and third row?
or did I store these samples by row or by column when I saved the data?!.

The difficulty is compounded with added dimensions.

Xarray objects eliminate much of the mental overhead by allowing indexing using dimension names instead of axes numbers:

In [18]:

Copied!

da.isel(lat=0, lon=0).plot();
da.isel(lat=0, lon=0).plot();

Slicing is also possible similarly:

In [19]:

Copied!

da.isel(time=slice(0, 20), lat=0, lon=0).plot();
da.isel(time=slice(0, 20), lat=0, lon=0).plot();

{note}
Using the `isel` method, the user can choose/slice the specific elements from a Dataset or DataArray.

So far, we have explored positional indexing, which relies on knowing the exact indices. But, what if you wanted to select data specifically for a particular latitude? It becomes challenging to determine the corresponding indices in such cases. Xarray reduce this complexity by introducing label-based indexing.

Label-based Indexing¶

To select data by coordinate labels instead of integer indices we can use sel instead of isel:

For example, let's select all data for Lat 25 °N and Lon 210 °E using sel :

In [20]:

hide-output

Copied!

da.sel(lat=25, lon=210).plot();
da.sel(lat=25, lon=210).plot();

Similarly we can do slicing or filter a range using the .slice function:

Dropping using `drop_sel`¶

If instead of selecting data we want to drop it, we can use drop_sel method with syntax similar to sel:

So far, all the above will require us to specify exact coordinate values, but what if we don't have the exact values? We can use nearest neighbor lookups to address this issue:

Nearest Neighbor Lookups¶

The label based selection methods sel() support method and tolerance keyword argument. The method parameter allows for enabling nearest neighbor (inexact) lookups by use of the methods ffill (propagate last valid index forward), backfill or nearest:

tolerance argument limits the maximum distance for valid matches with an inexact lookup:

{tip}
All of these indexing methods work on the dataset too!

We can also use these methods to index all variables in a dataset simultaneously, returning a new dataset:

Datetime Indexing¶

Datetime indexing is a critical feature when working with time series data, which is a common occurrence in atmospheric and environmental sciences. Essentially, datetime indexing allows you to select data points or a series of data points that correspond to certain date or time criteria. This becomes essential for time-series analysis where the date or time information associated with each data point can be as critical as the data point itself.

Let's see some of the techniques to perform datetime indexing in Xarray:

Selecting data based on single datetime¶

Let's say we have a Dataset ds and we want to select data at a particular date and time, for instance, '2013-01-01' at 6AM. We can do this by using the sel (select) method, like so:

By default, datetime selection will return a range of values that match the provided string. For e.g. time="2013-01-01" will return all timestamps for that day (4 of them here):

We can use this feature to select all points in a year:

or a month:

Selecting data for a range of dates¶

Now, let's say we want to select data between a certain range of dates. We can still use the sel method, but this time we will combine it with slice:

The slice function takes two arguments, start and stop, to make a slice that includes these endpoints. When we use slice with the sel method, it provides an efficient way to select a range of dates. The above example shows the usage of slice for datetime indexing.

Indexing with a DatetimeIndex or date string list¶

Another technique is to use a list of datetime objects or date strings for indexing. For example, you could select data for specific, non-contiguous dates like this:

Fancy indexing based on year, month, day, or other datetime components¶

In addition to the basic datetime indexing techniques, Xarray also supports "fancy" indexing options, which can provide more flexibility and efficiency in your data analysis tasks. You can directly access datetime components such as year, month, day, hour, etc. using the .dt accessor. Here is an example of selecting all data points from July across all years:

Or, if you wanted to select data from a specific day of each month, you could use:

Xarray Computation¶

Similar to Pandas, xarray supports different kinds of arithmetic operations

In [39]:

Copied!

ds['air_C'] = ds['air'] - 273.15
ds['air_C'] = ds['air'] - 273.15

In [43]:

Copied!

(ds['air_C'].mean(dim='time') >  0).plot()
(ds['air_C'].mean(dim='time') >  0).plot()

Out[43]:

<matplotlib.collections.QuadMesh at 0x194bbaa30>

In [44]:

Copied!

# Apply the mask:
ds.where(ds['air_C'].mean(dim='time') >  0).mean(dim = 'time').air.plot()
# Apply the mask:
ds.where(ds['air_C'].mean(dim='time') >  0).mean(dim = 'time').air.plot()

Out[44]:

<matplotlib.collections.QuadMesh at 0x194c9b100>

Broadcasting: expanding data¶

In [45]:

Copied!

ds['air_C'].mean(dim='time').plot()
ds['air_C'].mean(dim='time').plot()

Out[45]:

<matplotlib.collections.QuadMesh at 0x194dd1130>

In [46]:

Copied!

(ds['air_C'].mean(dim='time')*np.cos(np.deg2rad(ds.lat))).plot()
(ds['air_C'].mean(dim='time')*np.cos(np.deg2rad(ds.lat))).plot()

Out[46]:

<matplotlib.collections.QuadMesh at 0x194e8ccd0>

In [47]:

Copied!

# Make cos(lat) two dimensional
(xr.ones_like(ds['air_C'].mean(dim='time'))*np.cos(np.deg2rad(ds.lat))).plot()
# Make cos(lat) two dimensional
(xr.ones_like(ds['air_C'].mean(dim='time'))*np.cos(np.deg2rad(ds.lat))).plot()

Out[47]:

<matplotlib.collections.QuadMesh at 0x194f7a3d0>

In [48]:

Copied!

# Make cos(lat) two dimensional
(np.cos(np.deg2rad(ds.lat))*xr.ones_like(ds.lon)).plot()
# Make cos(lat) two dimensional
(np.cos(np.deg2rad(ds.lat))*xr.ones_like(ds.lon)).plot()

Out[48]:

<matplotlib.collections.QuadMesh at 0x194aa2f10>

Xarray Groupby and Resample¶

groupby¶

In [50]:

Copied!

# seasonal groups
ds.groupby("time.season")
# seasonal groups
ds.groupby("time.season")

Out[50]:

DatasetGroupBy, grouped over 'season'
4 groups with labels 'DJF', 'JJA', 'MAM', 'SON'.

The seasons are out of order (they are alphabetically sorted). This is a common annoyance. The solution is to use .sel to change the order of labels

In [53]:

Copied!

seasonal_mean.air.plot(col="season")
seasonal_mean.air.plot(col="season")

Out[53]:

<xarray.plot.facetgrid.FacetGrid at 0x1950926a0>

In [54]:

Copied!

# Make the figure to 2x2 plots:

# facet the seasonal_mean
seasonal_mean.air.plot(col="season", col_wrap=2);
# Make the figure to 2x2 plots:

# facet the seasonal_mean
seasonal_mean.air.plot(col="season", col_wrap=2);

In [55]:

Copied!

# Calculate zonal average
seasonal_mean.air.mean("lon").plot.line(hue="season", y="lat");
# Calculate zonal average
seasonal_mean.air.mean("lon").plot.line(hue="season", y="lat");

In [ ]:

resample¶

In [57]:

Copied!

ds.resample(time="M").mean().mean(dim =['lat','lon']).air.plot()
ds.resample(time="M").mean().mean(dim =['lat','lon']).air.plot()

Out[57]:

[<matplotlib.lines.Line2D at 0x1954c0730>]

In [ ]:

Lecture 8: Xarray for multidimensional gridded data¶

Xarray Data Model¶

Xarray functionality¶

Xarray data structures¶

DataArray¶

Install Xarray:¶

Xarray Data Structure¶

What is included in xarray Dataset?¶

DataArray¶

Named dimensions¶

Coordinates¶

Attributes¶

Working with Labelled Data¶

Selecting Data (Indexing)¶

Position-based Indexing¶

Positional Indexing Using Dimension Names¶

Label-based Indexing¶

Dropping using drop_sel¶

Nearest Neighbor Lookups¶

Datetime Indexing¶

Selecting data based on single datetime¶

Selecting data for a range of dates¶

Indexing with a DatetimeIndex or date string list¶

Fancy indexing based on year, month, day, or other datetime components¶

Xarray Computation¶

Broadcasting: expanding data¶

Xarray Groupby and Resample¶

groupby¶

resample¶

What is included in xarray `Dataset`?¶

Dropping using `drop_sel`¶