ObservationRecipe

When handling weather and climate data, it can be tedious and error-prone when setting up the observation for calibration with EnsembleKalmanProcesses (or EKP for short). As such, ClimaCalibrate provides recipes for estimating the noise covariance matrix and packaging it, together with the observation and its metadata, into an EKP.Observation.

We start with a SampleCollection, which the SampleBuilder module produces by flattening one or more ClimaAnalysis.OutputVars into a matrix of samples. ObservationRecipe then takes that SampleCollection, estimates the noise covariance matrix from the samples, selects the sample at the specified index as the observation, and builds the EKP.Observation used in the calibration.

How do I use this to set up observation for calibration with EKP?

All functions assume that any data preprocessing is done with ClimaAnalysis.

Covariance Estimators

There are currently three covariance estimators, ScalarCovariance, SeasonalDiagonalCovariance, and SVDplusDCovariance, which are subtypes of AbstractCovarianceEstimator. Each estimates the noise covariance matrix from the matrix of samples in a SampleCollection. ScalarCovariance approximates the observation noise covariance as a scalar diagonal matrix. SeasonalDiagonalCovariance approximates the observation noise covariance as a diagonal of variances across all the seasons for each observation, neglecting correlations between points. SVDplusDCovariance additionally approximates the correlations between points from, often limited, time series observations. Because SeasonalDiagonalCovariance and SVDplusDCovariance estimate variances and correlations across samples, they require at least two samples (columns), whereas ScalarCovariance works with a single sample.

Necessary data preprocessing

In most cases, the OutputVars represent time series data of summary statistics. For example, to compute seasonal averages of a OutputVar, one can use ClimaAnalysis.average_season_across_time, which will produce a OutputVar that can be used with either SeasonalDiagonalCovariance or SVDplusDCovariance.

import ClimaAnalysis

obs_var = ClimaAnalysis.OutputVar(
    "precip.mon.mean.nc",
    "precip",
    new_start_date = start_date,
    shift_by = Dates.firstdayofmonth,
)

# -- preprocessing for units, times, grid, etc. --

seasonal_averages = ClimaAnalysis.average_season_across_time(obs_var)

Observation

After preprocessing the OutputVars so that they represent time series data of summary statistics, you build the samples with functions provided by the SampleBuilder module and pass a covariance estimator, the resulting SampleCollection, and the index of the sample to use as the observation to ObservationRecipe.observation, as shown below. See the Sample Builder page for the details of build_samples_by_times.

import ClimaAnalysis
import EnsembleKalmanProcesses as EKP
import ClimaCalibrate
import ClimaCalibrate.ObservationRecipe
import ClimaCalibrate.SampleBuilder

# Vars are OutputVars preprocessed to ensure consistent units, times,
# and grid as the diagnostics produced from the model.
# In this example, we want to calibrate with seasonal averages, so we use
# ClimaAnalysis.average_season_across_time
vars = ClimaAnalysis.average_season_across_time.(vars)

# We need the start and end dates of each sample. To find these, we can use the
# function below. In this example, the dates in `vars` are all the same. For
# debugging, it is helpful to use `ClimaAnalysis.dates(var)`.
sample_date_ranges =
    ObservationRecipe.seasonally_aligned_yearly_sample_date_ranges(first(vars))

# Build the samples (one per date range)
sample_collection =
    SampleBuilder.build_samples_by_times(
        vars,
        sample_date_ranges;
        FT = Float32
    )

# We choose SVDplusDCovariance.
covar_estimator = ObservationRecipe.SVDplusDCovariance(
    model_error_scale = Float32(0.05),
    regularization = Float32(1e-6),
)

# Finally, we form the observation, using the first sample as the observation
# in the calibration
obs = ObservationRecipe.observation(covar_estimator, sample_collection, 1)

Metadata

When creating an observation with observation, metadata is extracted from the OutputVars and attached to the observation. The metadata for each observation can be accessed with EKP.get_metadata(obs::EKP.Observation) and the metadata for each iteration can be accessed with ClimaCalibrate.get_metadata_for_nth_iteration. The metadata can be used with ClimaAnalysis.unflatten to reconstruct the original OutputVar before flattening. See the ClimaAnalysis documentation about ClimaAnalysis.FlatVar for more information.

Debugging observational and simulation data

When setting up a calibration, it may be helpful to visualize the EKP.Observations or inspect the observational data and metadata together. To help with this, ObservationRecipe provides several functions that reconstruct the underlying data back into OutputVars:

• reconstruct_vars reconstructs the samples of an EKP.Observation as a vector of OutputVars.
• reconstruct_diag_cov reconstructs the diagonal of the covariance matrix of an EKP.Observation as a vector of OutputVars (only supported for diagonal covariance matrices).
• reconstruct_g reconstructs the G ensemble matrix of the itth iteration as a matrix of OutputVars.
• reconstruct_g_mean reconstructs the mean forward model evaluation of the itth iteration as a vector of OutputVars.
• reconstruct_g_mean_final reconstructs the mean forward model evaluation of the last iteration as a vector of OutputVars.

# obs is an EKP.Observation
# ekp is the EKP.EnsembleKalmanProcess
# it is the iteration index

ObservationRecipe.reconstruct_vars(obs)
# Reconstructing the diagonal of a covariance matrix as an OutputVar is only
# supported for diagonal covariance matrices
ObservationRecipe.reconstruct_diag_cov(obs)
ObservationRecipe.reconstruct_g(ekp, it)
ObservationRecipe.reconstruct_g_mean(ekp, it)
ObservationRecipe.reconstruct_g_mean_final(ekp)

Creating custom covariance estimators

In the cases where the provided covariance estimators are not sufficient, it is possible to create your own covariance estimator using the functionality provided by ClimaCalibrate and ClimaAnalysis.

The steps are:

1. Define a struct that subtypes ObservationRecipe.AbstractCovarianceEstimator. Any fields of this struct will be available when implementing the covariance method.
2. Implement a method of covariance that dispatches on your struct, with the signature ObservationRecipe.covariance(estimator::YourType, sample_collection). It must return a noise covariance matrix.

What you can use

The sample_collection argument is a SampleCollection storing a matrix of samples and a matrix of metadata. Here's a collection of helpful functions from the SampleBuilder module when creating the covariance matrix.

• get_samples: return the full sample matrix.
• get_metadata: return the full metadata matrix. Each column holds one ClimaAnalysis.Var.Metadata per variable.
• num_samples: return the number of samples.

The covariance you return must be a square matrix whose side length equals the number of rows of get_samples(sample_collection) and must not contain NaN or Inf.

Example: per-variable constant variance

The following estimator gives each variable its own constant variance, filling that variable's block of a diagonal covariance matrix.

import ClimaAnalysis
import ClimaCalibrate.ObservationRecipe
import ClimaCalibrate.SampleBuilder
import LinearAlgebra: Diagonal

# One variance per variable, in the same order as the OutputVars used to build
# the samples.
struct PerVariableScalar <: ObservationRecipe.AbstractCovarianceEstimator
    variances::Vector{Float64}
end

function ObservationRecipe.covariance(
    estimator::PerVariableScalar,
    sample_collection,
)
    FT = eltype(SampleBuilder.get_samples(sample_collection))
    # Only the flattened lengths are needed here, and those are the same across
    # the samples, so the first column is enough
    metadata = view(SampleBuilder.get_metadata(sample_collection), :, 1)
    # Repeat each variance over its variable's flattened block, then stack.
    diag_cov = reduce(
        vcat,
        fill(FT(v), ClimaAnalysis.flattened_length(m)) for
        (v, m) in zip(estimator.variances, metadata)
    )
    return Diagonal(diag_cov)
end

Using it is the same as for the built-in estimators:

julia> sample_collectionSampleCollection (9×1 matrix of Float32)
1 sample(s), each 9 value(s) from 2 variable(s)
Short name  Units   Indices  Dimensions
----------------------------------------------
pr          mm/day  1:6      lat (2), time (3)
rsut        W m-2   7:9      lat (3), time (1)
julia> estimator = PerVariableScalar([0.1, 0.5])Main.PerVariableScalar([0.1, 0.5])
julia> obs = ObservationRecipe.observation(estimator, sample_collection, 1);
julia> EKP.get_covs(obs)1-element Vector{LinearAlgebra.Diagonal{Float32, Vector{Float32}}}:
 [0.1 0.0 … 0.0 0.0; 0.0 0.1 … 0.0 0.0; … ; 0.0 0.0 … 0.5 0.0; 0.0 0.0 … 0.0 0.5]

Frequently asked questions

Q: I need to compute g_ensemble and I do not know how the data of the OutputVars is flattened.

A: When forming the sample, the data in a OutputVar is flattened using ClimaAnalysis.flatten. See ClimaAnalysis.flatten in the ClimaAnalysis documentation for more information. The order of the variables in the observation is the same as the order of the OutputVars when creating the EKP.Observation using ObservationRecipe.observation. If you are using ObservationRecipe, it is recommended that you also use GEnsembleBuilder which simplifies building the G ensemble matrix.

Q: How is the name of the observation determined?

A: By default, the name of the observation is determined by the short name in the attributes of the OutputVar. If there are multiple OutputVars, then the name is all the short names separated by semicolons. If a OutputVar has no short name, the empty string will be used instead and a warning is emitted. You can override this by passing the name keyword argument to ObservationRecipe.observation.

Q: What is regularization and model_error_scale when making a covariance matrix?

A: The model error scale and regularization terms are used to inflate the diagonal of the observation covariance matrix to reflect estimates of measurement error. You can add a fixed percentage inflation of the noise due to the model error to the covariance matrix with the model_error_scale keyword argument. Additionally, to prevent very small variance along the diagonal of the covariance matrix, you can add a regularization with the regularization keyword argument. For SVDplusDCovariance, the regularization keyword argument can also be a QuantileRegularization, which sets the regularization from a quantile of the model error scale instead of a fixed value.

Q: How do I apply latitude weighting to the covariance matrix?

A: All three covariance estimators accept a use_latitude_weights keyword argument. This accounts for the varying area of grid cells with latitude. The min_cosd_lat keyword argument (default 0.1) sets the minimum value of cosd(lat) used in the weight, which prevents very small values along the diagonal that can cause issues when inverting the covariance matrix. This requires the OutputVars to have a latitude dimension.