Handling dimensions with scores

Contents

Interactive online version: Binder badge. Download notebook.

Handling dimensions with scores#

Scores is primarily designed to work with xarray data. xarray allows us to work with multidimensional data with just a single implementation of a metric. In this notebook we illustrate how scores can work with data with different dimensions.

We recommend that you visit the xarray documentation to learn about xarray if you aren’t familiar with xarray.

First, let’s illustrate how scores can be used with data with different dimensions using synthetic data.

[2]:

import xarray as xr
import numpy as np
import pandas as pd

from scores.continuous import mse

np.random.seed(100)

[3]:

# First, let's set up some coordinates
valid_times = pd.date_range("2024-01-01", "2024-01-31")
lead_time = np.arange(1, 11)
station_numbers = np.arange(1, 11)
lat = np.arange(-90, 90, 15)
lon = np.arange(-180, 180, 15)

Data could be 1-dimensional (e.g., a timeseries at an airport)

[4]:

obs_timeseries = xr.DataArray(
    data=np.random.uniform(0, 20, size=len(valid_times)), dims=["valid_time"], coords={"valid_time": valid_times}
)
fcst_timeseries = obs_timeseries + np.random.normal(size=len(valid_times))
fcst_timeseries

[4]:

<xarray.DataArray (valid_time: 31)>
array([ 9.99728569,  5.79631977,  8.08231187, 16.79159751,  1.66155591,
        2.92840713, 14.57085403, 18.35566678,  4.26985841, 11.75686432,
       16.98228183,  3.20109898,  3.40046656,  3.00603842,  2.70310169,
       20.72364935, 15.20432613,  2.92782807, 13.96422444,  5.58507007,
       10.3729014 , 20.04247224, 16.48540033,  7.30003296,  1.93230337,
        6.16384669, -0.54614964,  4.17452229, 15.22369956, -0.23038042,
       13.50482056])
Coordinates:
  * valid_time  (valid_time) datetime64[ns] 2024-01-01 2024-01-02 ... 2024-01-31
[16]:

# Calculate MSE
mse(fcst_timeseries, obs_timeseries)

[16]:

<xarray.DataArray ()>
array(1.21810508)

Data could be 2-dimensional. For example, forecasts and observations at weather stations might have the dimensions valid_time and station_number.

[6]:

obs_stations = xr.DataArray(
    data=np.random.uniform(0, 20, size=[len(valid_times), len(station_numbers)]),
    dims=["valid_time", "station_number"],
    coords={"valid_time": valid_times, "station_number": station_numbers},
)
fcst_stations = obs_stations + np.random.normal(size=[len(valid_times), len(station_numbers)])
fcst_stations

[6]:

<xarray.DataArray (valid_time: 31, station_number: 10)>
array([[ 7.45491025e+00,  1.13596085e+00,  9.55356309e+00,
         9.25567774e+00,  1.95808607e+00,  1.16998470e+01,
         2.26261486e+01,  7.50361829e+00,  4.88419520e+00,
         1.62100534e+01],
       [ 1.46161421e+01,  6.26155468e+00,  1.26185279e+01,
         1.00887950e+01,  5.93870112e+00,  2.34637784e+00,
         6.79355816e+00,  8.70726799e+00,  1.51762850e+01,
         4.19755632e+00],
       [ 1.33905197e+01,  4.89084256e+00,  1.50774715e+01,
         1.43755762e+01,  1.78662887e+01,  1.29870635e+01,
         6.77606683e+00,  1.49704571e+01,  1.72497642e+01,
         1.45631540e+01],
       [ 1.75553287e+01,  1.96506539e+01,  4.20757053e+00,
         1.03412492e+00,  2.22919631e+00,  1.87921640e+01,
         2.03776121e+01,  3.77641474e+00,  6.30355657e+00,
         1.03712249e+01],
       [ 5.97467154e+00,  9.83056417e+00,  1.27625714e+01,
         9.66740698e+00,  1.33053421e+01,  1.22353186e+01,
        -9.89714911e-01,  4.94891473e+00,  1.18427111e+01,
         9.44907435e+00],
...
       [ 1.38757704e+01,  1.83863686e+01,  1.03139981e+01,
         1.00442693e+01,  1.28952417e+01,  1.22352772e+01,
         9.71756996e+00,  6.57696791e+00,  1.10432330e+01,
         2.03478908e+01],
       [ 1.94901078e+01,  7.72110555e+00,  1.70172641e+01,
         1.16445538e+01,  4.56336210e+00,  1.42529800e+00,
        -6.77478357e-01,  1.55661625e+00,  3.83132564e+00,
         1.72249589e+01],
       [ 5.71416646e+00,  1.63894575e+01,  9.23409025e+00,
         3.54024807e+00,  1.63192885e+01,  1.82041529e+01,
         1.93933734e+01,  1.31416482e+01,  1.43947255e+01,
         1.08362201e+01],
       [ 9.24551408e+00,  2.16468824e+01,  5.20689186e-01,
         1.94635556e+01,  5.16272310e+00,  1.26412064e+01,
         1.04677155e+01,  1.23261874e+01,  1.18393369e+01,
         4.79039382e+00],
       [ 4.16452314e+00,  2.43650695e+00,  4.29515334e-01,
         2.05556446e+01,  6.45836848e+00,  2.87991004e+00,
         6.15523884e+00,  7.99932744e+00,  1.24329978e+01,
         2.98415626e+00]])
Coordinates:
  * valid_time      (valid_time) datetime64[ns] 2024-01-01 ... 2024-01-31
  * station_number  (station_number) int64 1 2 3 4 5 6 7 8 9 10

If we want to calculate the MSE aggregated across time and space, there are a couple of different ways that we can call the MSE function that will give us the same result.

  1. Not using any dim handling arg will aggregate data and take the mean across all dimensions

  2. Use reduce_dims="all" to take the mean across all dimensions

[7]:

mse(fcst_stations, obs_stations)

[7]:

<xarray.DataArray ()>
array(1.07490488)
[8]:

mse(fcst_stations, obs_stations, reduce_dims="all")

[8]:

<xarray.DataArray ()>
array(1.07490488)

If we want to calculate the MSE aggregated across space (station_number dimension), but not time (valid_time dimension), then there are also two options.

  1. Use preserve_dims=["valid_time"], or

  2. Use reduce_dims=["station_number']

when the score is called.

You cannot supply both the reduce_dims and preserve_dims args at the same time.

[9]:

mse(fcst_stations, obs_stations, preserve_dims=["valid_time"])

[9]:

<xarray.DataArray (valid_time: 31)>
array([1.44266719, 0.60170873, 1.76365055, 0.96647668, 2.17553097,
       1.04172737, 0.59422138, 0.3350336 , 1.28693609, 1.19349423,
       0.73360716, 1.36191764, 0.84777202, 0.34682081, 0.9437102 ,
       1.38537663, 0.73261705, 0.3140486 , 0.75206193, 0.72817601,
       1.01798424, 1.62898803, 1.53654807, 1.39244355, 2.54658843,
       0.97392879, 0.93276668, 0.55192585, 0.40048297, 1.5661211 ,
       1.22671862])
Coordinates:
  * valid_time  (valid_time) datetime64[ns] 2024-01-01 2024-01-02 ... 2024-01-31
[10]:

mse(fcst_stations, obs_stations, reduce_dims=["station_number"])

[10]:

<xarray.DataArray (valid_time: 31)>
array([1.44266719, 0.60170873, 1.76365055, 0.96647668, 2.17553097,
       1.04172737, 0.59422138, 0.3350336 , 1.28693609, 1.19349423,
       0.73360716, 1.36191764, 0.84777202, 0.34682081, 0.9437102 ,
       1.38537663, 0.73261705, 0.3140486 , 0.75206193, 0.72817601,
       1.01798424, 1.62898803, 1.53654807, 1.39244355, 2.54658843,
       0.97392879, 0.93276668, 0.55192585, 0.40048297, 1.5661211 ,
       1.22671862])
Coordinates:
  * valid_time  (valid_time) datetime64[ns] 2024-01-01 2024-01-02 ... 2024-01-31

Data can have many more dimensions. Additionally, forecast and observations can have different dimensions (e.g. forecast with a lead_time dimension). When this occurs, forecast and observations will be broadcast against each other following numpy broadcasting rules.

We will show an example with gridded data.

The forecast has dimensions [valid_time, lat, lon, lead_time], while the observations only have dimensions [valid_time, lat, lon]

[11]:

obs_grid = xr.DataArray(
    data=np.random.uniform(0, 20, size=[len(valid_times), len(lat), len(lon)]),
    dims=["valid_time", "lat", "lon"],
    coords={"valid_time": valid_times, "lat": lat, "lon": lon},
)
fcst_grid = obs_grid + np.random.normal(size=[len(valid_times), len(lat), len(lon)])
fcst_grid = fcst_grid * xr.DataArray(data=np.arange(1, 1.5, 0.05), dims="lead_time", coords={"lead_time": lead_time})

View the dimensions of the synthetic gridded data that we created

[12]:

fcst_grid.dims

[12]:

('valid_time', 'lat', 'lon', 'lead_time')

[13]:

obs_grid.dims

[13]:

('valid_time', 'lat', 'lon')

Perhaps we want to calculate MSE while preserving the lead_time dimension to see how performance changes by lead time.

[14]:

mse(fcst_grid, obs_grid, preserve_dims=["lead_time"])

[14]:

<xarray.DataArray (lead_time: 10)>
array([ 1.01199121,  1.47394501,  2.60395301,  4.40201522,  6.86813163,
       10.00230224, 13.80452706, 18.27480608, 23.41313931, 29.21952673])
Coordinates:
  * lead_time  (lead_time) int64 1 2 3 4 5 6 7 8 9 10

We may want to preserve all dimensions. This is particularly important for calculating confidence intervals. We can do so using two approaches:

  1. Set preserve_dims="all, or

  2. Set preserve_dims=["valid_time", "lat", "lon", "lead_time]

when the score is called.

[15]:

mse(fcst_grid, obs_grid, preserve_dims="all")

[15]:

<xarray.DataArray (valid_time: 31, lat: 12, lon: 24, lead_time: 10)>
array([[[[1.28647777e-02, 2.16718450e-01, 6.68531203e-01, ...,
          6.64698118e+00, 8.58654842e+00, 1.07740747e+01],
         [2.04897998e-04, 4.02193738e-01, 1.64529154e+00, ...,
          2.04774151e+01, 2.67671667e+01, 3.38980272e+01],
         [1.07744316e-01, 1.67325809e+00, 5.10238192e+00, ...,
          5.02021518e+01, 6.48129360e+01, 8.12873302e+01],
         ...,
         [2.45338991e-01, 5.76386971e-01, 1.04670498e+00, ...,
          5.48734552e+00, 6.79328372e+00, 8.23849196e+00],
         [1.04428114e+00, 1.63368497e+00, 2.35442318e+00, ...,
          7.92813013e+00, 9.43687470e+00, 1.10769537e+01],
         [2.87449866e+00, 1.20173728e+00, 2.47049673e-01, ...,
          6.24471824e+00, 9.59847327e+00, 1.36703021e+01]],

        [[8.70007369e-01, 1.94902328e+00, 3.45739298e+00, ...,
          1.74395483e+01, 2.15240408e+01, 2.60378870e+01],
         [1.52170687e-01, 3.65791013e-01, 6.71617242e-01, ...,
          3.58383693e+00, 4.44289857e+00, 5.39416612e+00],
         [1.77571143e+00, 1.01779157e+00, 4.69437238e-01, ...,
          8.71148603e-01, 1.58018748e+00, 2.49879189e+00],
...
         [8.04121243e-02, 4.52654992e-03, 2.22042891e-02, ...,
          1.51404269e+00, 2.09310031e+00, 2.76572124e+00],
         [6.94671713e+00, 1.30781221e+01, 2.11330960e+01, ...,
          9.02614978e+01, 1.09857885e+02, 1.31377840e+02],
         [7.69873368e+00, 1.12775249e+01, 1.55373621e+01, ...,
          4.70522362e+01, 5.53983486e+01, 6.44255070e+01]],

        [[9.21378635e-01, 3.37361470e-01, 4.07118698e-02, ...,
          2.86797733e+00, 4.29553311e+00, 6.01045646e+00],
         [4.41616170e-01, 5.71178527e-01, 7.17383596e-01, ...,
          1.69804960e+00, 1.94411094e+00, 2.20681498e+00],
         [2.58370376e+00, 2.73965251e+00, 2.90017080e+00, ...,
          3.77130518e+00, 3.95924064e+00, 4.15174564e+00],
         ...,
         [1.30112375e+00, 1.81307488e-01, 8.35578666e-02, ...,
          1.49258094e+01, 2.09604596e+01, 2.80171765e+01],
         [1.14134585e+00, 2.59280871e+00, 4.63154419e+00, ...,
          2.36343108e+01, 2.91966820e+01, 3.53463258e+01],
         [4.76874467e-05, 6.65423020e-01, 2.68427235e+00, ...,
          3.30806289e+01, 4.32203222e+01, 5.47134895e+01]]]])
Coordinates:
  * valid_time  (valid_time) datetime64[ns] 2024-01-01 2024-01-02 ... 2024-01-31
  * lat         (lat) int64 -90 -75 -60 -45 -30 -15 0 15 30 45 60 75
  * lon         (lon) int64 -180 -165 -150 -135 -120 ... 105 120 135 150 165
  * lead_time   (lead_time) int64 1 2 3 4 5 6 7 8 9 10

Final notes#

Some scores will handle dimensions in specific ways. For example when using scores.probability.crps_for_ensemble, you need to specify the name of the ensemble member dimension with the crps_for_ensemble arg. Another example is scores.continuous.flip_flop_index_proportion_exceeding where you specify what dimension to sample across with the sampling_dim arg. In both cases the reduce_dims and preserve_dims are still available.

Contents