Evaluate against sparse observations¶

This guide will walk through the basics of how to evaluate a gridded forecast against a sparse ground truth dataset, in this case METAR weather station data.

Sparse data is tabular in nature. The METAR data on the WeatherBench cloud bucket uses Parquet, a tabular cloud-optimized data format.

import numpy as np
from weatherbenchX import interpolations
from weatherbenchX import binning
from weatherbenchX import aggregation
from weatherbenchX.metrics import base as metrics_base
from weatherbenchX.metrics import deterministic
from weatherbenchX.data_loaders import sparse_parquet
from weatherbenchX.data_loaders import xarray_loaders

variables = ['2m_temperature', '10m_wind_speed']
init_times = np.array(['2020-01-01T00', '2020-01-01T12'], dtype='datetime64[ns]')
lead_times = np.array([6, 12], dtype='timedelta64[h]').astype('timedelta64[ns]')

target_data_loader = sparse_parquet.METARFromParquet(
    path='gs://weatherbench2/datasets/metar/metar-timeNominal-by-month/',
    variables=variables,
    partitioned_by='month',
    time_dim='timeNominal',
    add_nan_mask=True
)

target_chunk = target_data_loader.load_chunk(init_times, lead_times)
target_chunk

{'2m_temperature': <xarray.DataArray '2m_temperature' (index: 33026)> Size: 132kB
 array([273.15, 282.15, 291.15, ..., 274.15, 268.15, 272.85], dtype=float32)
 Coordinates:
     latitude     (index) float32 132kB -77.87 -53.8 -33.38 ... 46.55 49.82 49.83
     longitude    (index) float32 132kB 167.0 292.2 289.2 ... 299.0 285.0 295.7
     elevation    (index) float32 132kB 8.0 22.0 476.0 141.0 ... 13.0 381.0 53.0
     stationName  (index) object 264kB 'NZCM' 'SAWE' 'SCEL' ... 'CWUK' 'CWBY'
     valid_time   (index) datetime64[ns] 264kB 2020-01-01T06:00:00 ... 2020-01-02
     init_time    (index) datetime64[ns] 264kB 2020-01-01 ... 2020-01-01T12:00:00
     lead_time    (index) timedelta64[ns] 264kB 06:00:00 06:00:00 ... 12:00:00
   * index        (index) int64 264kB 0 1 2 3 4 ... 33021 33022 33023 33024 33025
     mask         (index) bool 33kB True True True True ... True True True True,
 '10m_wind_speed': <xarray.DataArray '10m_wind_speed' (index: 33026)> Size: 132kB
 array([4.1, 5.1, 2.1, ..., 9.3, 2.1, 2.1], dtype=float32)
 Coordinates:
     latitude     (index) float32 132kB -77.87 -53.8 -33.38 ... 46.55 49.82 49.83
     longitude    (index) float32 132kB 167.0 292.2 289.2 ... 299.0 285.0 295.7
     elevation    (index) float32 132kB 8.0 22.0 476.0 141.0 ... 13.0 381.0 53.0
     stationName  (index) object 264kB 'NZCM' 'SAWE' 'SCEL' ... 'CWUK' 'CWBY'
     valid_time   (index) datetime64[ns] 264kB 2020-01-01T06:00:00 ... 2020-01-02
     init_time    (index) datetime64[ns] 264kB 2020-01-01 ... 2020-01-01T12:00:00
     lead_time    (index) timedelta64[ns] 264kB 06:00:00 06:00:00 ... 12:00:00
   * index        (index) int64 264kB 0 1 2 3 4 ... 33021 33022 33023 33024 33025
     mask         (index) bool 33kB True True True True ... True True True True}

The data now only has an index dimension with init and lead time being non-dimension coordinates. This is because for each init and lead time, there will be a different number of observations.

When loaded raw, the data has missing values. One way to deal with this is using the NaN mask as done in this example. Another way to deal with this for sparse data is to set the dropna argument in the data loader. With the split_variables argument set to True, this will then return a dictionary of DataArrays for each variable with NaNs dropped. With split_variables=False, this will drop all values where any variable is NaN.

There also are several measurement for the same station for the same nominal time (which always is a full hour; the raw observations will have different observation times). We cloud also drop those using remove_duplicates.

Next, we need to interpolate the gridded forecast to the station locations. For this, we can use the InterpolateToReferenceCoords interpolation, which is passed to the gridded prediction loader.

interpolation = interpolations.InterpolateToReferenceCoords(
    method='nearest',
    wrap_longitude=True
)

prediction_path = 'gs://weatherbench2/datasets/hres/2016-2022-0012-64x32_equiangular_conservative.zarr'
prediction_data_loader = xarray_loaders.PredictionsFromXarray(
    path=prediction_path,
    variables=variables,
    interpolation=interpolation,
)

In addition to the init and lead times, we now also pass the target chunk to the prediction loader. The interpolator will then use the coordinates on the target chunks to interpolate the gridded predictions.

prediction_chunk = prediction_data_loader.load_chunk(init_times, lead_times, reference=target_chunk)
prediction_chunk

{'2m_temperature': <xarray.DataArray '2m_temperature' (index: 33026)> Size: 132kB
 array([271.56082, 282.60428, 289.23752, ..., 273.2698 , 269.18704,
        273.2698 ], dtype=float32)
 Coordinates:
     init_time    (index) datetime64[ns] 264kB 2020-01-01 ... 2020-01-01T12:00:00
     lead_time    (index) timedelta64[ns] 264kB 06:00:00 06:00:00 ... 12:00:00
     longitude    (index) float32 132kB 167.0 292.2 289.2 ... 299.0 285.0 295.7
     latitude     (index) float32 132kB -77.87 -53.8 -33.38 ... 46.55 49.82 49.83
     elevation    (index) float32 132kB 8.0 22.0 476.0 141.0 ... 13.0 381.0 53.0
     stationName  (index) object 264kB 'NZCM' 'SAWE' 'SCEL' ... 'CWUK' 'CWBY'
     valid_time   (index) datetime64[ns] 264kB 2020-01-01T06:00:00 ... 2020-01-02
   * index        (index) int64 264kB 0 1 2 3 4 ... 33021 33022 33023 33024 33025
     mask         (index) bool 33kB True True True True ... True True True True
 Attributes:
     long_name:      2 metre temperature
     short_name:     t2m
     standard_name:  unknown
     units:          K,
 '10m_wind_speed': <xarray.DataArray '10m_wind_speed' (index: 33026)> Size: 132kB
 array([2.2358263, 5.0393696, 5.332086 , ..., 7.8019896, 3.5218635,
        7.8019896], dtype=float32)
 Coordinates:
     init_time    (index) datetime64[ns] 264kB 2020-01-01 ... 2020-01-01T12:00:00
     lead_time    (index) timedelta64[ns] 264kB 06:00:00 06:00:00 ... 12:00:00
     longitude    (index) float32 132kB 167.0 292.2 289.2 ... 299.0 285.0 295.7
     latitude     (index) float32 132kB -77.87 -53.8 -33.38 ... 46.55 49.82 49.83
     elevation    (index) float32 132kB 8.0 22.0 476.0 141.0 ... 13.0 381.0 53.0
     stationName  (index) object 264kB 'NZCM' 'SAWE' 'SCEL' ... 'CWUK' 'CWBY'
     valid_time   (index) datetime64[ns] 264kB 2020-01-01T06:00:00 ... 2020-01-02
   * index        (index) int64 264kB 0 1 2 3 4 ... 33021 33022 33023 33024 33025
     mask         (index) bool 33kB True True True True ... True True True True}

Now the target and prediction data is aligned and we can proceed as usual, apart from the aggregation. Because lead_time is now not a dimension any more, the aggregator would simply average out all lead_times. Typically, however, we want results as a function of lead time. To achieve this, we need to add a binning instance that adds a lead_time bin.

metrics = {
  'rmse': deterministic.RMSE(),
  'mae': deterministic.MAE(),
}
statistics = metrics_base.compute_unique_statistics_for_all_metrics(
  metrics, prediction_chunk, target_chunk
)
bin_by = [binning.ByExactCoord('lead_time')]
aggregator = aggregation.Aggregator(
  reduce_dims=['index'],
  bin_by=bin_by,
  masked=True
)
aggregation_state = aggregator.aggregate_statistics(statistics)