header

Benchmark MIOST sea surface height maps (Geostrophic mode only)

2023-04-27 MIOST_SSH_GEOS_BENCHMARK_DEMO


Authors: CLS & Datlas Copyright: 2023 CLS & Datlas License: MIT

Benchmark of MIOST sea surface height maps

The notebook aims to evaluate the sea surface height maps produced by the MIOST system.

<h5> These maps are equivalent to the SEALEVEL_GLO_PHY_L4_MY_008_047 product distributed by the Copernicus Marine Service, except that a nadir altimeter (SARAL/Altika, SEALEVEL_GLO_PHY_L3_MY_008_062 product) has been excluded from the mapping. </h5>
    <h5> We provide below a demonstration of the validation of these maps against the independent SSH data from the Saral/AltiKa altimeter distributed by CMEMS </h5>

General Note 1: Execute each cell through the button from the top MENU (or keyboard shortcut Shift + Enter). General Note 2: If, for any reason, the kernel is not working anymore, in the top MENU, click on the button. Then, in the top MENU, click on “Cell” and select “Run All Above Selected Cell”. ***

[ ]:

[1]:
from glob import glob
import numpy as np
import os
[2]:
import sys
sys.path.append('..')
from src.mod_plot import *
from src.mod_stat import *
from src.mod_spectral import *
[3]:
import logging
logger = logging.getLogger()
logger.setLevel(logging.INFO)
[4]:
time_min = '2019-01-01'                                        # time min for analysis
time_max = '2019-12-31'                                        # time max for analysis
output_dir = '../results'                                      # output directory path
os.system(f'mkdir -p {output_dir}')
stat_output_filename = f'{output_dir}/stat_sla_miost_geos.nc'  # output statistical analysis filename
lambda_min = 65.                                               # minimun spatial scale in kilometer to consider on the filtered signal
lambda_max = 500.                                              # maximum spatial scale in kilometer to consider on the filtered signal
psd_output_filename = f'{output_dir}/psd_sla_miost_geos.nc'    # output spectral analysis filename
segment_lenght = 1000.                                         # spectral parameer: along-track segment lenght in kilometer to consider in the spectral analysis

Sea Surface Height from Saral/AltiKa

[5]:
list_of_file = sorted(glob('/work/ALT/swot/aval/wisa/data_challenge_ose/data/independent_alongtrack/alg/2019/*.nc'))
ds_alg = xr.open_mfdataset(list_of_file, combine='nested', concat_dim='time')
ds_alg = ds_alg.where((ds_alg.time >= np.datetime64(time_min)) & (ds_alg.time <=  np.datetime64(time_max)), drop=True)
ds_alg = ds_alg.sortby('time')
ds_alg
[5]:
<xarray.Dataset>
Dimensions:         (time: 14437696)
Coordinates:
  * time            (time) datetime64[ns] 2019-01-01T00:04:07.003014144 ... 2...
    longitude       (time) float64 dask.array<chunksize=(44621,), meta=np.ndarray>
    latitude        (time) float64 dask.array<chunksize=(44621,), meta=np.ndarray>
Data variables:
    cycle           (time) float64 dask.array<chunksize=(44621,), meta=np.ndarray>
    track           (time) float64 dask.array<chunksize=(44621,), meta=np.ndarray>
    sla_unfiltered  (time) float32 dask.array<chunksize=(44621,), meta=np.ndarray>
    sla_filtered    (time) float32 dask.array<chunksize=(44621,), meta=np.ndarray>
    dac             (time) float32 dask.array<chunksize=(44621,), meta=np.ndarray>
    ocean_tide      (time) float32 dask.array<chunksize=(44621,), meta=np.ndarray>
    internal_tide   (time) float32 dask.array<chunksize=(44621,), meta=np.ndarray>
    lwe             (time) float32 dask.array<chunksize=(44621,), meta=np.ndarray>
    mdt             (time) float32 dask.array<chunksize=(44621,), meta=np.ndarray>
    tpa_correction  (time) float32 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0 0.0
Attributes: (12/44)
    Conventions:                     CF-1.6
    Metadata_Conventions:            Unidata Dataset Discovery v1.0
    cdm_data_type:                   Swath
    comment:                         Sea surface height measured by altimeter...
    contact:                         servicedesk.cmems@mercator-ocean.eu
    creator_email:                   servicedesk.cmems@mercator-ocean.eu
    ...                              ...
    summary:                         SSALTO/DUACS Delayed-Time Level-3 sea su...
    time_coverage_duration:          P24H11M29.716548S
    time_coverage_end:               2019-01-01T23:36:52Z
    time_coverage_resolution:        P1S
    time_coverage_start:             2018-12-31T23:25:22Z
    title:                           DT Altika Drifting Phase Global Ocean Al...

Sea Surface Height maps to evaluate

[6]:
list_of_maps = sorted(glob('/work/ALT/swot/aval/wisa/data_challenge_ose/data/maps/MIOST_geos_global_allsat-alg/*.nc'))
ds_maps = xr.open_mfdataset(list_of_maps, combine='nested', concat_dim='time')
ds_maps = ds_maps.sel(time=slice(time_min, time_max))
ds_maps
[6]:
<xarray.Dataset>
Dimensions:    (longitude: 3600, latitude: 1701, time: 365)
Coordinates:
  * longitude  (longitude) float64 0.0 0.1 0.2 0.3 ... 359.6 359.7 359.8 359.9
  * latitude   (latitude) float32 -80.05 -79.95 -79.85 ... 89.75 89.85 89.95
  * time       (time) datetime64[ns] 2019-01-01 2019-01-02 ... 2019-12-31
Data variables:
    sla        (time, latitude, longitude) float64 dask.array<chunksize=(1, 1701, 3600), meta=np.ndarray>
    ugosa      (time, latitude, longitude) float64 dask.array<chunksize=(1, 1701, 3600), meta=np.ndarray>
    vgosa      (time, latitude, longitude) float64 dask.array<chunksize=(1, 1701, 3600), meta=np.ndarray>
    adt        (time, latitude, longitude) float64 dask.array<chunksize=(1, 1701, 3600), meta=np.ndarray>
    ugos       (time, latitude, longitude) float64 dask.array<chunksize=(1, 1701, 3600), meta=np.ndarray>
    vgos       (time, latitude, longitude) float64 dask.array<chunksize=(1, 1701, 3600), meta=np.ndarray>

2.1 Interpolate sea surface height maps onto along-track positions

[7]:
ds_interp = run_interpolation(ds_maps, ds_alg)
ds_interp = ds_interp.dropna('time')
ds_interp
2023-04-28 13:52:30 INFO     fetch data from 2019-01-01 00:00:00 to 2019-02-01 00:00:00
2023-04-28 13:52:41 INFO     fetch data from 2019-01-31 00:00:00 to 2019-03-01 00:00:00
2023-04-28 13:52:51 INFO     fetch data from 2019-02-28 00:00:00 to 2019-04-01 00:00:00
2023-04-28 13:53:02 INFO     fetch data from 2019-03-31 00:00:00 to 2019-05-01 00:00:00
2023-04-28 13:53:13 INFO     fetch data from 2019-04-30 00:00:00 to 2019-06-01 00:00:00
2023-04-28 13:53:23 INFO     fetch data from 2019-05-31 00:00:00 to 2019-07-01 00:00:00
2023-04-28 13:53:34 INFO     fetch data from 2019-06-30 00:00:00 to 2019-08-01 00:00:00
2023-04-28 13:53:45 INFO     fetch data from 2019-07-31 00:00:00 to 2019-09-01 00:00:00
2023-04-28 13:53:56 INFO     fetch data from 2019-08-31 00:00:00 to 2019-10-01 00:00:00
2023-04-28 13:54:06 INFO     fetch data from 2019-09-30 00:00:00 to 2019-11-01 00:00:00
2023-04-28 13:54:16 INFO     fetch data from 2019-10-31 00:00:00 to 2019-12-01 00:00:00
2023-04-28 13:54:27 INFO     fetch data from 2019-11-30 00:00:00 to 2019-12-31 00:00:00
[7]:
<xarray.Dataset>
Dimensions:            (time: 14370249)
Coordinates:
  * time               (time) datetime64[ns] 2019-01-01T00:04:07.003014144 .....
Data variables: (12/13)
    cycle              (time) float64 126.0 126.0 126.0 ... 136.0 136.0 136.0
    track              (time) float64 9.0 9.0 9.0 9.0 ... 408.0 408.0 408.0
    sla_unfiltered     (time) float32 0.127 0.101 0.07 ... 0.007 0.056 0.01
    sla_filtered       (time) float32 0.117 0.128 0.142 ... 0.03 0.031 0.027
    dac                (time) float32 0.303 0.299 0.296 ... 0.172 0.171 0.172
    ocean_tide         (time) float32 -0.091 -0.096 -0.099 ... -0.304 -0.314
    ...                 ...
    lwe                (time) float32 -0.019 -0.019 -0.019 ... 0.014 0.014 0.016
    mdt                (time) float32 -0.104 -0.106 -0.108 ... -1.138 -1.154
    tpa_correction     (time) float32 0.0 0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0
    longitude          (time) float64 64.42 64.33 64.25 ... 246.7 246.6 246.3
    latitude           (time) float64 69.9 69.95 70.01 ... -69.27 -69.32 -69.55
    msla_interpolated  (time) float64 0.1579 0.1586 0.1596 ... 0.02084 0.03111

2.2 Compute grid boxes statistics & statistics by regime (coastal, offshore low variability, offshore high variability)

[8]:
compute_stat_scores(ds_interp, lambda_min, lambda_max, stat_output_filename)
2023-04-28 13:54:42 INFO     Compute mapping error all scales
2023-04-28 13:54:42 INFO     Compute mapping error for scales between 65.0 and 500.0 km
2023-04-28 13:56:37 INFO     Compute binning statistics
2023-04-28 13:56:54 INFO     Compute statistics by oceanic regime
2023-04-28 13:57:35 INFO     Stat file saved as: ../results/stat_sla_miost_geos.nc
[9]:
# Plot gridded stats
# Hvplot
# plot_stat_score_map(stat_output_filename)
# Matplotlib
plot_stat_score_map_png(stat_output_filename)
../_images/gallery_ssh_scores_MIOST_geos_19_0.png
[10]:
plot_stat_score_timeseries(stat_output_filename)
[10]:
[11]:
plot_stat_by_regimes(stat_output_filename)
[11]:
mapping_err_var [m²] sla_unfiltered_var [m²] mapping_err_filtered_var [m²] sla_filtered_var [m²] var_score_allscale var_score_filtered
coastal 0.002576 0.008303 0.000549 0.001813 0.689822 0.697257
offshore_highvar 0.002921 0.053263 0.000957 0.021750 0.945153 0.956015
offshore_lowvar 0.001221 0.006299 0.000225 0.001775 0.806072 0.873085
equatorial_band 0.001449 0.005975 0.000309 0.000536 0.757493 0.423876
arctic 0.002477 0.007084 0.000413 0.001019 0.650382 0.594350
antarctic 0.003350 0.002344 0.000675 0.000753 -0.429089 0.103416

2.3 Compute Spectral scores

[ ]:
compute_psd_scores(ds_interp, psd_output_filename, lenght_scale=segment_lenght)
[12]:
compute_psd_scores_v2(ds_interp, psd_output_filename, lenght_scale=segment_lenght)
2023-04-28 13:57:47 INFO     Segment computation...
2023-04-28 13:58:02 INFO     Spectral analysis...
2023-04-28 14:12:31 INFO     Saving ouput...
2023-04-28 14:12:39 INFO     PSD file saved as: ../results/psd_sla_miost_geos.nc
[13]:
# Plot effective resolution
# Hvplot
# plot_effective_resolution(psd_output_filename)
# Matplotlib
plot_effective_resolution_png(psd_output_filename)
../_images/gallery_ssh_scores_MIOST_geos_25_0.png
[ ]:
plot_psd_scores(psd_output_filename)
[ ]:

[ ]: