Core dimensions#
Previously we learned to use apply_ufunc
on simple functions that acted element by element.
Here we move on to slightly more complex functions like np.mean
that can act along a subset of an input array’s dimensions.
Such operations involve the concept of “core dimensions”.
Our learning goals are:
Learn how to identify “core dimensions” for the function you’re applying.
Learn that “core dimensions” are automatically moved or transposed to the end of the array.
Introduction#
For using more complex operations that consider some array values collectively,
it’s important to understand the idea of core dimensions.
Usually, they correspond to the fundamental dimensions over
which an operation is defined, e.g., the summed axis in np.sum
. One way to think about core dimensions
is to consider the smallest dimensionality of data that the function acts on.
Important
A good clue that core dimensions are needed is the presence of an axis
argument on the
corresponding NumPy function.
Setup#
%xmode minimal
import numpy as np
import xarray as xr
# limit the amount of information printed to screen
xr.set_options(display_expand_data=False)
np.set_printoptions(threshold=10, edgeitems=2)
Exception reporting mode: Minimal
Let’s load a dataset
ds = xr.tutorial.load_dataset("air_temperature")
ds
<xarray.Dataset> Size: 31MB Dimensions: (lat: 25, time: 2920, lon: 53) Coordinates: * lat (lat) float32 100B 75.0 72.5 70.0 67.5 65.0 ... 22.5 20.0 17.5 15.0 * lon (lon) float32 212B 200.0 202.5 205.0 207.5 ... 325.0 327.5 330.0 * time (time) datetime64[ns] 23kB 2013-01-01 ... 2014-12-31T18:00:00 Data variables: air (time, lat, lon) float64 31MB 241.2 242.5 243.5 ... 296.2 295.7 Attributes: Conventions: COARDS title: 4x daily NMC reanalysis (1948) description: Data is from NMC initialized reanalysis\n(4x/day). These a... platform: Model references: http://www.esrl.noaa.gov/psd/data/gridded/data.ncep.reanaly...
Reducing with np.mean
#
Let’s write a function that computes the mean along time
for a provided xarray object.
This function requires one core dimension time
. For ds.air
note that time
is the 0th axis.
ds.air.dims
('time', 'lat', 'lon')
get_axis_num
is a useful method.
ds.air.get_axis_num("time")
0
np.mean(ds.air, axis=ds.air.get_axis_num("time"))
<xarray.DataArray 'air' (lat: 25, lon: 53)> Size: 11kB 260.4 260.2 259.9 259.5 259.0 258.6 ... 298.0 297.9 297.8 297.3 297.3 297.3 Coordinates: * lat (lat) float32 100B 75.0 72.5 70.0 67.5 65.0 ... 22.5 20.0 17.5 15.0 * lon (lon) float32 212B 200.0 202.5 205.0 207.5 ... 325.0 327.5 330.0
np.mean(ds.air.data, axis=0)
array([[260.37644178, 260.18305137, ..., 251.93811644, 253.43804795],
[262.73439384, 262.79397603, ..., 251.58575685, 254.35926027],
...,
[298.12920205, 297.93700685, ..., 296.7770274 , 296.44383562],
[298.36615068, 298.38573973, ..., 297.28144521, 297.30510274]])
Let’s try to use apply_ufunc
to replicate np.mean(ds.air.data, axis=0)
xr.apply_ufunc(
# function to apply
np.mean,
# object with data to pass to function
ds,
# keyword arguments to pass to np.mean
kwargs={"axis": 0},
)
ValueError: applied function returned data with an unexpected number of dimensions. Received 2 dimension(s) but expected 3 dimensions with names ('time', 'lat', 'lon'), from:
array([[260.376442, 260.183051, 259.886627, ..., 250.815901, 251.938116,
253.438048],
[262.734394, 262.793976, 262.749339, ..., 249.755904, 251.585757,
254.35926 ],
[264.768764, 264.327308, 264.061695, ..., 250.60789 , 253.58351 ,
257.715599],
...,
[297.649863, 296.953332, 296.629315, ..., 296.810925, 296.287962,
295.816455],
[298.129202, 297.937007, 297.470394, ..., 296.859548, 296.777027,
296.443836],
[298.366151, 298.38574 , 298.114144, ..., 297.338205, 297.281445,
297.305103]])
The error here
applied function returned data with unexpected number of dimensions.
Received 2 dimension(s) but expected 3 dimensions with names: ('time', 'lat', 'lon')
means that while np.mean
did indeed reduce one dimension, we did not tell apply_ufunc
that this would happen. That is, we need to specify the core dimensions on the input.
Do that by passing a list of dimension names for each input object. For this function we have one input : ds
and with a single core dimension "time"
so we have input_core_dims=[["time"]]
xr.apply_ufunc(
np.mean,
ds,
# specify core dimensions as a list of lists
# here 'time' is the core dimension on `ds`
input_core_dims=[
["time"], # core dimension for ds
],
kwargs={"axis": 0},
)
ValueError: size of dimension 'lat' on inputs was unexpectedly changed by applied function from 25 to 53. Only dimensions specified in ``exclude_dims`` with xarray.apply_ufunc are allowed to change size. The data returned was:
array([[279.398 , 279.0572, 279.0104, ..., 279.63 , 279.398 , 279.27 ],
[279.6664, 279.538 , 279.2808, ..., 279.934 , 279.666 , 279.354 ],
[279.6612, 279.7296, 279.5508, ..., 280.534 , 280.318 , 279.882 ],
...,
[279.9508, 279.7756, 279.682 , ..., 279.802 , 279.766 , 279.426 ],
[280.3152, 280.27 , 280.1976, ..., 280.346 , 280.342 , 279.97 ],
[280.6624, 280.7976, 280.814 , ..., 280.778 , 280.834 , 280.482 ]])
This next error is a little confusing.
size of dimension 'lat' on inputs was unexpectedly changed by applied function from 25 to 53.
Only dimensions specified in ``exclude_dims`` with xarray.apply_ufunc are allowed to change size.
A good trick here is to pass a little wrapper function to apply_ufunc
instead and inspect the shapes of data received by the wrapper.
def wrapper(array, **kwargs):
print(f"received {type(array)} shape: {array.shape}, kwargs: {kwargs}")
result = np.mean(array, **kwargs)
print(f"result.shape: {result.shape}")
return result
xr.apply_ufunc(
wrapper,
ds,
# specify core dimensions as a list of lists
# here 'time' is the core dimension on `ds`
input_core_dims=[["time"]],
kwargs={"axis": 0},
)
received <class 'numpy.ndarray'> shape: (25, 53, 2920), kwargs: {'axis': 0}
result.shape: (53, 2920)
ValueError: size of dimension 'lat' on inputs was unexpectedly changed by applied function from 25 to 53. Only dimensions specified in ``exclude_dims`` with xarray.apply_ufunc are allowed to change size. The data returned was:
array([[279.398 , 279.0572, 279.0104, ..., 279.63 , 279.398 , 279.27 ],
[279.6664, 279.538 , 279.2808, ..., 279.934 , 279.666 , 279.354 ],
[279.6612, 279.7296, 279.5508, ..., 280.534 , 280.318 , 279.882 ],
...,
[279.9508, 279.7756, 279.682 , ..., 279.802 , 279.766 , 279.426 ],
[280.3152, 280.27 , 280.1976, ..., 280.346 , 280.342 , 279.97 ],
[280.6624, 280.7976, 280.814 , ..., 280.778 , 280.834 , 280.482 ]])
Now we see the issue:
received <class 'numpy.ndarray'> shape: (25, 53, 2920), kwargs: {'axis': 0}
result.shape: (53, 2920)
The time
dimension is of size 2920
and is now the last axis of the array but was initially the first axis
ds.air.get_axis_num("time")
0
Important
This illustrates an important concept. Arrays are transposed so that core dimensions are at the end.
With apply_ufunc
, core dimensions are recognized by name, and then moved to
the last dimension of any input arguments before applying the given function.
This means that for functions that accept an axis
argument, you usually need
to set axis=-1
Such behaviour means that our functions (like wrapper
or np.mean
) do not need to know the exact order of dimensions. They can rely on the core dimensions being at the end allowing us to write very general code!
We can fix our apply_ufunc
call by specifying axis=-1
instead.
def wrapper(array, **kwargs):
print(f"received {type(array)} shape: {array.shape}, kwargs: {kwargs}")
result = np.mean(array, **kwargs)
print(f"result.shape: {result.shape}")
return result
xr.apply_ufunc(
wrapper,
ds,
input_core_dims=[["time"]],
kwargs={"axis": -1},
)
received <class 'numpy.ndarray'> shape: (25, 53, 2920), kwargs: {'axis': -1}
result.shape: (25, 53)
<xarray.Dataset> Size: 11kB Dimensions: (lat: 25, lon: 53) Coordinates: * lat (lat) float32 100B 75.0 72.5 70.0 67.5 65.0 ... 22.5 20.0 17.5 15.0 * lon (lon) float32 212B 200.0 202.5 205.0 207.5 ... 325.0 327.5 330.0 Data variables: air (lat, lon) float64 11kB 260.4 260.2 259.9 ... 297.3 297.3 297.3
Use apply_ufunc
to apply scipy.integrate.trapezoid
along the time
axis.
Solution to Exercise 24
import scipy as sp
import scipy.integrate
xr.apply_ufunc(scipy.integrate.trapezoid, ds, input_core_dims=[["time"]], kwargs={"axis": -1})