AlongTrack Evaluation

AlongTrack Evaluation#

In this notebook, we look at how to evaluate alongtrack data. In the OSSE setting in previous notebooks, we saw that we can evaluate our method on the alongtracks and then compare our learned field with the field from the simulation. In the OSE setting, we don’t have the original field. So we need to evaluate the field directly on the leave-one-out alongtrack. We demonstrate some simple statistics in the real and spectral space that we can perform.

import autoroot
import typing as tp
from dataclasses import dataclass
import functools as ft
import numpy as np
import pandas as pd
import xarray as xr
import einops
from metpy.units import units
import pint_xarray
import xarray_dataclasses as xrdataclass
from oceanbench._src.geoprocessing.validation import validate_latlon, validate_time, decode_cf_time, validate_ssh
from oceanbench._src.preprocessing.alongtrack import alongtrack_ssh
from oceanbench._src.geoprocessing.subset import where_slice
from oceanbench._src.preprocessing.alongtrack import remove_swath_dimension
import matplotlib.pyplot as plt
import matplotlib.colors as colors
import matplotlib.ticker as ticker
import seaborn as sns

sns.reset_defaults()
sns.set_context(context="talk", font_scale=0.7)

%load_ext autoreload
%autoreload 2

Field#

Now, we will load a field from a method that interpolated the data.

# !ls /gpfswork/rech/yrf/commun/data_challenges/dc20a_osse/staging/natl60/
!ls /gpfswork/rech/yrf/commun/data_challenges/dc21a_ose/test/results

leaderboard.csv			  OSE_ssh_mapping_DUACS.nc
OSE_ssh_mapping_4dvarNet_2022.nc  OSE_ssh_mapping_DYMOST.nc
OSE_ssh_mapping_4dvarNet.nc	  OSE_ssh_mapping_MIOST.nc
OSE_ssh_mapping_BASELINE.nc	  results.csv
OSE_ssh_mapping_BFN.nc

def preprocess_field(da):
    
    # validate coordinates
    da = validate_latlon(da)
    da = validate_time(da)
    da = validate_ssh(da)
    
    # select region + time period
    da = da.sel(
        time=slice("2017-01-01", "2017-12-31"),
        lon=slice(-65., -55.),
        lat=slice(33., 43.),
    )
    
    return da

Real Statistics#

Statistics#

# psd_calc.psd_welch_score(ds_segments, "ssh", "ssh_interp", delta_x=delta_x, nperseg=nperseg).score

Spectral Statistics#

Power Spectrum#

# SSH AlongTrack
ds_psd = psd_calc.psd_welch(ds_segments, variable="ssh", delta_x=delta_x, nperseg=nperseg).ssh

# SSH Interpolated Map
ds_psd_interp = psd_calc.psd_welch(ds_segments, variable="ssh_interp", delta_x=delta_x, nperseg=nperseg).ssh_interp

# PSD Error
ds_psd_error = psd_calc.psd_welch_error(ds_segments, variable_ref="ssh", variable="ssh_interp", delta_x=delta_x, nperseg=nperseg).error

# PSD Score
ds_psd_score = psd_calc.psd_welch_score(ds_segments, variable_ref="ssh", variable="ssh_interp", delta_x=delta_x, nperseg=nperseg).score

ds_psd

<xarray.DataArray 'ssh' (wavenumber: 79)>
array([1.02373037e+01, 7.14277725e+01, 4.22874489e+01, 2.21006603e+01,
       8.25816345e+00, 2.71950841e+00, 1.40674925e+00, 6.88535631e-01,
       3.68134439e-01, 2.17998356e-01, 1.30187139e-01, 8.60557631e-02,
       6.66696206e-02, 5.22705801e-02, 4.14400324e-02, 3.09867822e-02,
       2.78125107e-02, 2.32093651e-02, 2.11372841e-02, 1.81758646e-02,
       1.80478077e-02, 1.65228639e-02, 1.55075006e-02, 1.27541861e-02,
       1.12466970e-02, 1.24145634e-02, 1.20474407e-02, 1.09734135e-02,
       1.00068841e-02, 1.06369033e-02, 1.12501457e-02, 1.16684642e-02,
       9.71117616e-03, 9.58298892e-03, 9.61237121e-03, 9.55455657e-03,
       9.76762362e-03, 1.02562634e-02, 9.34584532e-03, 8.92794598e-03,
       9.51405521e-03, 9.44814924e-03, 9.20096785e-03, 9.09277238e-03,
       8.38163868e-03, 8.41708854e-03, 8.41432996e-03, 7.14715524e-03,
       7.98082165e-03, 8.83140601e-03, 7.99081381e-03, 8.02045316e-03,
       7.78246671e-03, 8.94837547e-03, 8.09417851e-03, 8.00068583e-03,
       7.34653696e-03, 8.02976545e-03, 7.93727767e-03, 7.48019805e-03,
       7.49919470e-03, 7.87938386e-03, 7.77373137e-03, 7.88705889e-03,
       8.31951480e-03, 7.64988083e-03, 7.48078013e-03, 7.67526543e-03,
       8.13236553e-03, 7.90031720e-03, 6.98836194e-03, 7.61722075e-03,
       7.75263691e-03, 7.75971264e-03, 7.80881057e-03, 7.37115275e-03,
       7.29903765e-03, 7.00414460e-03, 3.10720433e-03], dtype=float32)
Coordinates:
  * wavenumber  (wavenumber) float64 0.0 0.001004 0.002007 ... 0.07728 0.07829

xarray.DataArray

'ssh'

wavenumber: 79

10.24 71.43 42.29 22.1 8.258 ... 0.007371 0.007299 0.007004 0.003107

array([1.02373037e+01, 7.14277725e+01, 4.22874489e+01, 2.21006603e+01,
       8.25816345e+00, 2.71950841e+00, 1.40674925e+00, 6.88535631e-01,
       3.68134439e-01, 2.17998356e-01, 1.30187139e-01, 8.60557631e-02,
       6.66696206e-02, 5.22705801e-02, 4.14400324e-02, 3.09867822e-02,
       2.78125107e-02, 2.32093651e-02, 2.11372841e-02, 1.81758646e-02,
       1.80478077e-02, 1.65228639e-02, 1.55075006e-02, 1.27541861e-02,
       1.12466970e-02, 1.24145634e-02, 1.20474407e-02, 1.09734135e-02,
       1.00068841e-02, 1.06369033e-02, 1.12501457e-02, 1.16684642e-02,
       9.71117616e-03, 9.58298892e-03, 9.61237121e-03, 9.55455657e-03,
       9.76762362e-03, 1.02562634e-02, 9.34584532e-03, 8.92794598e-03,
       9.51405521e-03, 9.44814924e-03, 9.20096785e-03, 9.09277238e-03,
       8.38163868e-03, 8.41708854e-03, 8.41432996e-03, 7.14715524e-03,
       7.98082165e-03, 8.83140601e-03, 7.99081381e-03, 8.02045316e-03,
       7.78246671e-03, 8.94837547e-03, 8.09417851e-03, 8.00068583e-03,
       7.34653696e-03, 8.02976545e-03, 7.93727767e-03, 7.48019805e-03,
       7.49919470e-03, 7.87938386e-03, 7.77373137e-03, 7.88705889e-03,
       8.31951480e-03, 7.64988083e-03, 7.48078013e-03, 7.67526543e-03,
       8.13236553e-03, 7.90031720e-03, 6.98836194e-03, 7.61722075e-03,
       7.75263691e-03, 7.75971264e-03, 7.80881057e-03, 7.37115275e-03,
       7.29903765e-03, 7.00414460e-03, 3.10720433e-03], dtype=float32)

Coordinates: (1)

wavenumber

(wavenumber)

float64

0.0 0.001004 ... 0.07728 0.07829

array([0.      , 0.001004, 0.002007, 0.003011, 0.004015, 0.005018, 0.006022,
       0.007026, 0.008029, 0.009033, 0.010037, 0.01104 , 0.012044, 0.013048,
       0.014051, 0.015055, 0.016059, 0.017062, 0.018066, 0.01907 , 0.020073,
       0.021077, 0.022081, 0.023084, 0.024088, 0.025092, 0.026095, 0.027099,
       0.028103, 0.029106, 0.03011 , 0.031114, 0.032117, 0.033121, 0.034125,
       0.035128, 0.036132, 0.037136, 0.038139, 0.039143, 0.040147, 0.04115 ,
       0.042154, 0.043158, 0.044161, 0.045165, 0.046169, 0.047172, 0.048176,
       0.04918 , 0.050183, 0.051187, 0.052191, 0.053194, 0.054198, 0.055202,
       0.056206, 0.057209, 0.058213, 0.059217, 0.06022 , 0.061224, 0.062228,
       0.063231, 0.064235, 0.065239, 0.066242, 0.067246, 0.06825 , 0.069253,
       0.070257, 0.071261, 0.072264, 0.073268, 0.074272, 0.075275, 0.076279,
       0.077283, 0.078286])

Indexes: (1)

wavenumber

PandasIndex

PandasIndex(Float64Index([                  0.0,  0.001003669810264894,
               0.002007339620529788, 0.0030110094307946823,
               0.004014679241059576,   0.00501834905132447,
               0.006022018861589365,  0.007025688671854258,
               0.008029358482119152,  0.009033028292384046,
                0.01003669810264894,  0.011040367912913834,
                0.01204403772317873,  0.013047707533443623,
               0.014051377343708517,   0.01505504715397341,
               0.016058716964238304,    0.0170623867745032,
               0.018066056584768092,  0.019069726395032988,
                0.02007339620529788,  0.021077066015562775,
               0.022080735825827667,  0.023084405636092563,
                0.02408807544635746,   0.02509174525662235,
               0.026095415066887246,  0.027099084877152138,
               0.028102754687417034,  0.029106424497681926,
                0.03011009430794682,  0.031113764118211713,
                0.03211743392847661,  0.033121103738741504,
                 0.0341247735490064,   0.03512844335927129,
               0.036132113169536184,   0.03713578297980108,
               0.038139452790065975,  0.039143122600330864,
                0.04014679241059576,  0.041150462220860655,
                0.04215413203112555,  0.043157801841390446,
               0.044161471651655335,   0.04516514146192023,
               0.046168811272185126,   0.04717248108245002,
                0.04817615089271492,  0.049179820702979805,
                 0.0501834905132447,  0.051187160323509596,
                0.05219083013377449,   0.05319449994403938,
               0.054198169754304276,   0.05520183956456917,
                0.05620550937483407,   0.05720917918509896,
                0.05821284899536385,   0.05921651880562875,
                0.06022018861589364,   0.06122385842615854,
                0.06222752823642343,   0.06323119804668832,
                0.06423486785695322,   0.06523853766721811,
                0.06624220747748301,    0.0672458772877479,
                 0.0682495470980128,   0.06925321690827768,
                0.07025688671854258,   0.07126055652880747,
                0.07226422633907237,   0.07326789614933726,
                0.07427156595960216,   0.07527523576986705,
                0.07627890558013195,   0.07728257539039685,
                0.07828624520066173],
             dtype='float64', name='wavenumber'))

Attributes: (0)

Figure#

class PlotPSDIsotropic:
        
    def init_fig(self, ax=None, figsize=None):
        if ax is None:
            figsize = (5,4) if figsize is None else figsize
            self.fig, self.ax = plt.subplots(figsize=figsize)
        else:
            self.ax = ax
            self.fig = plt.gcf()
        
    def plot_wavenumber(self, da,freq_scale=1.0, units=None, **kwargs):
        
        if units is not None:
            xlabel = f"Wavenumber [cycles/{units}]"
        else:
            xlabel = f"Wavenumber"
            
        dim = list(da.dims)[0]
        
        self.ax.plot(da[dim] * freq_scale, da, **kwargs)

        self.ax.set(
            yscale="log", xscale="log",
            xlabel=xlabel,
            ylabel=f"PSD [{da.name}]",
            xlim=[10**(-3) - 0.00025, 10**(-1) +0.025]
        )

        self.ax.legend()
        self.ax.grid(which="both", alpha=0.5)
        
    def plot_wavelength(self, da, freq_scale=1.0, units=None, **kwargs):
        
        if units is not None:
            xlabel = f"Wavelength [{units}]"
        else:
            xlabel = f"Wavelength"
        
        dim = list(da.dims)[0]
        self.ax.plot(1/(da[dim] * freq_scale), da, **kwargs)
        
        self.ax.set(
            yscale="log", xscale="log",
            xlabel=xlabel,
            ylabel=f"PSD [{da.name}]"
        )

        self.ax.xaxis.set_major_formatter("{x:.0f}")
        self.ax.invert_xaxis()
        
        self.ax.legend()
        self.ax.grid(which="both", alpha=0.5)
                
    def plot_both(self, da, freq_scale=1.0, units=None, **kwargs):
        
        if units is not None:
            xlabel = f"Wavelength [{units}]"
        else:
            xlabel = f"Wavelength"
        
        self.plot_wavenumber(da=da, units=units, freq_scale=freq_scale, **kwargs)
        
        self.secax = self.ax.secondary_xaxis(
            "top", functions=(lambda x: 1 / (x + 1e-20), lambda x: 1 / (x + 1e-20))
        )
        self.secax.xaxis.set_major_formatter("{x:.0f}")
        self.secax.set(xlabel=xlabel)

psd_iso_plot = PlotPSDIsotropic()
psd_iso_plot.init_fig()
psd_iso_plot.plot_both(
    ds_psd,
    freq_scale=1, 
    units="km",
    label="NADIR",
    color="black",
)
psd_iso_plot.plot_both(
    ds_psd_interp,
    freq_scale=1, 
    units="km",
    label="DUACS",
    color="tab:green",
)

# set custom bounds
psd_iso_plot.ax.set_xlim((10**(-3) - 0.00025, 10**(-1) +0.025))
psd_iso_plot.ax.set_ylabel("PSD [SSH]")

plt.show()

../../_images/1f1e516d4717bb019ddc305044653b8ecc3d12de0ae8c6c75c8e19099ddcfcb8.png

Power Spectrum Score#

\[ \text{PSD}(u,\hat{u}) = 1 - \frac{\mathcal{F}(u) - \mathcal{F}(\hat{u})}{\mathcal{F}(u)} \]

where \(\mathcal{F}\) is the power spectrum transformation.

from oceanbench._src.metrics.utils import find_intercept_1D, find_intercept_2D

Resolved Scale#

We can assign an arbitrary threshhold which we decide as the “resolved scale” for the quantity of interest. In this case, we decided 0.5 which is the minimum signal-to-noise ratio that we deem trustworthy.

find_intercept_1D(
            x=ds_psd_score.values, 
    y=1./(ds_psd_score.wavenumber), 
    level=0.5,
    kind="slinear",
    fill_value="extrapolate"
)        

array(152.67922533)

Figure#

class PlotPSDScoreIsotropic(PlotPSDIsotropic):
    
    def _add_score(
        self,
        da,
        freq_scale=1.0, 
        units=None, 
        threshhold: float=0.5, 
        threshhold_color="k",
        name=""
):
        
        dim = da.dims[0]
        self.ax.set(ylabel="PSD Score", yscale="linear")
        self.ax.set_ylim((0,1.0))
        self.ax.set_xlim((
            10**(-3) - 0.00025,
            10**(-1) +0.025,
        ))
        
        resolved_scale = freq_scale / find_intercept_1D(
            x=da.values, y=1./(da[dim].values+1e-15), level=threshhold
        )        
        self.ax.vlines(
            x=resolved_scale, 
            ymin=0, ymax=threshhold, 
            color=threshhold_color,
            linewidth=2, linestyle="--",
        )
        self.ax.hlines(
            y=threshhold,
            xmin=np.ma.min(np.ma.masked_invalid(da[dim].values * freq_scale)),
            xmax=resolved_scale, color=threshhold_color,
            linewidth=2, linestyle="--"
        )        
        label = f"{name}: {1/resolved_scale:.0f} {units} "
        self.ax.scatter(
            resolved_scale, threshhold,
            color=threshhold_color, marker=".",
            linewidth=5, label=label,
            zorder=3
        )
        
        
    def plot_score(
        self, 
        da, 
        freq_scale=1.0, 
        units=None, 
        threshhold: float=0.5, 
        threshhold_color="k",
        name="",
        **kwargs
    ):
        
        self.plot_both(da=da, freq_scale=freq_scale, units=units, **kwargs)
        self._add_score(
            da=da, 
            freq_scale=freq_scale,
            units=units,
            threshhold=threshhold, 
            threshhold_color=threshhold_color,
            name=name
        )

psd_iso_plot = PlotPSDScoreIsotropic()
psd_iso_plot.init_fig()
psd_iso_plot.plot_score(
    ds_psd_score,
    freq_scale=1, 
    units="km",
    name="DUACS",
    color="black",
    threshhold=0.5,
    threshhold_color="tab:blue"
    
)
plt.legend()
plt.show()

No artists with labels found to put in legend.  Note that artists whose label start with an underscore are ignored when legend() is called with no argument.

../../_images/b650ce049ed5b4ea0a66a0ff96134ac068458cdc7c66e28c8fda1f0347e08908.png

AlongTrack Evaluation

Contents

AlongTrack Evaluation#

Data#

Field#

Gridding#

Regridding#

AlongTrack Segments#

Real Statistics#

Statistics#

Spectral Statistics#

Power Spectrum#

Figure#

Power Spectrum Score#

Resolved Scale#

Figure#