Gulf Stream: Benchmark DUACS sea surface height maps

2023-04-27 DUACS_SSH_BENCHMARK_DEMO

---

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

Gulf Stream: Benchmark of DUACS sea surface height maps

The notebook aims to evaluate the sea surface height maps produced by the DUACS 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”. ***

Learning outcomes

At the end of this notebook you will know:

• How you can evaluated sea surface height maps with independent alongtrack data: statistical and spectral analysis

from glob import glob
import numpy as np
import os

import sys
sys.path.append('..')
from src.mod_plot import *
from src.mod_stat import *
from src.mod_spectral import *

import logging
logger = logging.getLogger()
logger.setLevel(logging.INFO)

0. Parameters

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}')

# Gulf Stream region
region = 'GS'
lon_min = 295                                          # domain min longitude
lon_max = 305                                          # domain max longitude
lat_min = 33.                                          # domain min latitude
lat_max = 43.                                          # domain max latitude
box_lonlat_GS = {'lon_min':lon_min,'lon_max':lon_max,'lat_min':lat_min,'lat_max':lat_max}

method_name = 'DUACS'

stat_output_filename = f'{output_dir}/stat_sla_duacs_geos_GS.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_duacs_geos_GS.nc'    # output spectral analysis filename
segment_lenght = 1000.                                         # spectral parameer: along-track segment lenght in kilometer to consider in the spectral analysis

1. Input files

Sea Surface Height from Saral/AltiKa

%%time
list_of_file = sorted(glob('../data/independent_alongtrack/indep_nadir_GS.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

CPU times: user 68.3 ms, sys: 12.9 ms, total: 81.2 ms
Wall time: 85.5 ms

<xarray.Dataset>
Dimensions:         (time: 37177)
Coordinates:
  * time            (time) datetime64[ns] 2019-01-01T22:10:01.073445888 ... 2...
    longitude       (time) float64 dask.array<chunksize=(37177,), meta=np.ndarray>
    latitude        (time) float64 dask.array<chunksize=(37177,), meta=np.ndarray>
Data variables:
    cycle           (time) float64 dask.array<chunksize=(37177,), meta=np.ndarray>
    track           (time) float64 dask.array<chunksize=(37177,), meta=np.ndarray>
    sla_unfiltered  (time) float32 dask.array<chunksize=(37177,), meta=np.ndarray>
    sla_filtered    (time) float32 dask.array<chunksize=(37177,), meta=np.ndarray>
    dac             (time) float32 dask.array<chunksize=(37177,), meta=np.ndarray>
    ocean_tide      (time) float32 dask.array<chunksize=(37177,), meta=np.ndarray>
    internal_tide   (time) float32 dask.array<chunksize=(37177,), meta=np.ndarray>
    lwe             (time) float32 dask.array<chunksize=(37177,), meta=np.ndarray>
    mdt             (time) float32 dask.array<chunksize=(37177,), meta=np.ndarray>
    tpa_correction  (time) float32 dask.array<chunksize=(37177,), meta=np.ndarray>
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...

xarray.Dataset

• Dimensions:
  • time: 37177
• Coordinates: (3)
  • time
    (time)
    datetime64[ns]
    2019-01-01T22:10:01.073445888 .....

    axis :

    T

    long_name :

    Time of measurement

    standard_name :

    time

    array(['2019-01-01T22:10:01.073445888', '2019-01-01T22:10:02.129446144',
           '2019-01-01T22:10:03.185446144', ..., '2019-12-30T09:33:55.091854080',
           '2019-12-30T09:33:56.147854080', '2019-12-30T09:33:57.203853824'],
          dtype='datetime64[ns]')

  • longitude
    (time)
    float64
    dask.array<chunksize=(37177,), meta=np.ndarray>

    long_name :

    Longitude of measurement

    standard_name :

    longitude

    units :

    degrees_east

    Array Chunk Bytes 290.45 kiB 290.45 kiB Shape (37177,) (37177,) Dask graph 1 chunks in 2 graph layers Data type float64 numpy.ndarray

  • latitude
    (time)
    float64
    dask.array<chunksize=(37177,), meta=np.ndarray>

    long_name :

    Latitude of measurement

    standard_name :

    latitude

    units :

    degrees_north

    Array Chunk Bytes 290.45 kiB 290.45 kiB Shape (37177,) (37177,) Dask graph 1 chunks in 2 graph layers Data type float64 numpy.ndarray

• Data variables: (10)
  • cycle
    (time)
    float64
    dask.array<chunksize=(37177,), meta=np.ndarray>

    long_name :

    Cycle the measurement belongs to

    units :

    1

    Array Chunk Bytes 290.45 kiB 290.45 kiB Shape (37177,) (37177,) Dask graph 1 chunks in 4 graph layers Data type float64 numpy.ndarray

  • track
    (time)
    float64
    dask.array<chunksize=(37177,), meta=np.ndarray>

    long_name :

    Track in cycle the measurement belongs to

    units :

    1

    Array Chunk Bytes 290.45 kiB 290.45 kiB Shape (37177,) (37177,) Dask graph 1 chunks in 4 graph layers Data type float64 numpy.ndarray

  • sla_unfiltered
    (time)
    float32
    dask.array<chunksize=(37177,), meta=np.ndarray>

    comment :

    The sea level anomaly is the sea surface height above mean sea surface height; the uncorrected sla can be computed as follows: [uncorrected sla]=[sla from product]+[dac]+[ocean_tide]+[internal_tide]-[lwe]; see the product user manual for details

    long_name :

    Sea level anomaly not-filtered not-subsampled with dac, ocean_tide and lwe correction applied

    standard_name :

    sea_surface_height_above_sea_level

    units :

    m

    Array Chunk Bytes 145.22 kiB 145.22 kiB Shape (37177,) (37177,) Dask graph 1 chunks in 4 graph layers Data type float32 numpy.ndarray

  • sla_filtered
    (time)
    float32
    dask.array<chunksize=(37177,), meta=np.ndarray>

    comment :

    The sea level anomaly is the sea surface height above mean sea surface height; the uncorrected sla can be computed as follows: [uncorrected sla]=[sla from product]+[dac]+[ocean_tide]+[internal_tide]-[lwe]; see the product user manual for details

    long_name :

    Sea level anomaly filtered not-subsampled with dac, ocean_tide and lwe correction applied

    standard_name :

    sea_surface_height_above_sea_level

    units :

    m

    Array Chunk Bytes 145.22 kiB 145.22 kiB Shape (37177,) (37177,) Dask graph 1 chunks in 4 graph layers Data type float32 numpy.ndarray

  • dac
    (time)
    float32
    dask.array<chunksize=(37177,), meta=np.ndarray>

    comment :

    The sla in this file is already corrected for the dac; the uncorrected sla can be computed as follows: [uncorrected sla]=[sla from product]+[dac]; see the product user manual for details

    long_name :

    Dynamic Atmospheric Correction

    units :

    m

    Array Chunk Bytes 145.22 kiB 145.22 kiB Shape (37177,) (37177,) Dask graph 1 chunks in 4 graph layers Data type float32 numpy.ndarray

  • ocean_tide
    (time)
    float32
    dask.array<chunksize=(37177,), meta=np.ndarray>

    comment :

    The sla in this file is already corrected for the ocean_tide; the uncorrected sla can be computed as follows: [uncorrected sla]=[sla from product]+[ocean_tide]; see the product user manual for details

    long_name :

    Ocean tide model

    units :

    m

    Array Chunk Bytes 145.22 kiB 145.22 kiB Shape (37177,) (37177,) Dask graph 1 chunks in 4 graph layers Data type float32 numpy.ndarray

  • internal_tide
    (time)
    float32
    dask.array<chunksize=(37177,), meta=np.ndarray>

    comment :

    The sla in this file is already corrected for the internal_tide; the uncorrected sla can be computed as follows: [uncorrected sla]=[sla from product]+[internal_tide]; see the product user manual for details

    long_name :

    Internal tide correction

    units :

    m

    Array Chunk Bytes 145.22 kiB 145.22 kiB Shape (37177,) (37177,) Dask graph 1 chunks in 4 graph layers Data type float32 numpy.ndarray

  • lwe
    (time)
    float32
    dask.array<chunksize=(37177,), meta=np.ndarray>

    comment :

    The sla in this file is already corrected for the lwe; the uncorrected sla can be computed as follows: [uncorrected sla]=[sla from product]-[lwe]; see the product user manual for details

    long_name :

    Long wavelength error

    units :

    m

    Array Chunk Bytes 145.22 kiB 145.22 kiB Shape (37177,) (37177,) Dask graph 1 chunks in 4 graph layers Data type float32 numpy.ndarray

  • mdt
    (time)
    float32
    dask.array<chunksize=(37177,), meta=np.ndarray>

    comment :

    The mean dynamic topography is the sea surface height above geoid; it is used to compute the absolute dynamic tyopography adt=sla+mdt

    long_name :

    Mean dynamic topography

    standard_name :

    sea_surface_height_above_geoid

    units :

    m

    Array Chunk Bytes 145.22 kiB 145.22 kiB Shape (37177,) (37177,) Dask graph 1 chunks in 4 graph layers Data type float32 numpy.ndarray

  • tpa_correction
    (time)
    float32
    dask.array<chunksize=(37177,), meta=np.ndarray>

    comment :

    This variable can be added to the SLA to correct for the observed instrumental drift during the lifetime of the TOPEX-A mission (the correction is null after this period). This is a global correction to be added a posteriori (and not before) on the global mean sea level estimate derived from the along-track product. It can be applied at regional or local scale as a best estimate (better than no correction, since the regional variation of the instrumental drift is unknown). See product manual for more details.

    long_name :

    TOPEX-A instrumental drift correction derived from altimetry and tide gauges global comparisons (WCRP Sea Level Budget Group, 2018)

    standard_name :

    sea_surface_height_above_sea_level

    units :

    m

    Array Chunk Bytes 145.22 kiB 145.22 kiB Shape (37177,) (37177,) Dask graph 1 chunks in 4 graph layers Data type float32 numpy.ndarray

• Indexes: (1)
  • time
    PandasIndex

    PandasIndex(DatetimeIndex(['2019-01-01 22:10:01.073445888',
                   '2019-01-01 22:10:02.129446144',
                   '2019-01-01 22:10:03.185446144',
                   '2019-01-01 22:10:04.241445888',
                   '2019-01-01 22:10:05.297445888',
                   '2019-01-01 22:10:06.353445888',
                   '2019-01-01 22:10:07.409445888',
                   '2019-01-01 22:10:08.465445888',
                   '2019-01-01 22:10:09.521446144',
                   '2019-01-01 22:10:10.577446144',
                   ...
                   '2019-12-30 09:33:47.699854080',
                   '2019-12-30 09:33:48.755854080',
                   '2019-12-30 09:33:49.811853824',
                   '2019-12-30 09:33:50.867853824',
                   '2019-12-30 09:33:51.923853824',
                   '2019-12-30 09:33:52.979853824',
                   '2019-12-30 09:33:54.035853824',
                   '2019-12-30 09:33:55.091854080',
                   '2019-12-30 09:33:56.147854080',
                   '2019-12-30 09:33:57.203853824'],
                  dtype='datetime64[ns]', name='time', length=37177, freq=None))

• Attributes: (44)

  Conventions :

  CF-1.6

  Metadata_Conventions :

  Unidata Dataset Discovery v1.0

  cdm_data_type :

  Swath

  comment :

  Sea surface height measured by altimeters referenced to the [1993, 2012] period; with additional corrections; the proposed sla is already corrected for dac, ocean_tide and lwe; [uncorrected sla]=[sla from product]+[dac]+[ocean_tide]-[lwe]

  contact :

  servicedesk.cmems@mercator-ocean.eu

  creator_email :

  servicedesk.cmems@mercator-ocean.eu

  creator_name :

  CMEMS - Sea Level Thematic Assembly Center

  creator_url :

  http://marine.copernicus.eu

  date_created :

  2021-09-08T13:54:51Z

  date_issued :

  2021-09-08T13:54:51Z

  date_modified :

  2021-09-08T13:54:51Z

  geospatial_lat_max :

  81.443179

  geospatial_lat_min :

  -78.158928

  geospatial_lat_resolution :

  0.0556490000000025

  geospatial_lat_units :

  degrees_north

  geospatial_lon_max :

  359.998793

  geospatial_lon_min :

  0.010756

  geospatial_lon_resolution :

  0.017816999999979544

  geospatial_lon_units :

  degrees_east

  geospatial_vertical_max :

  0.0

  geospatial_vertical_min :

  0.0

  geospatial_vertical_positive :

  down

  geospatial_vertical_resolution :

  point

  geospatial_vertical_units :

  m

  history :

  2021-09-08T13:54:51Z: Creation

  institution :

  CLS, CNES

  keywords :

  Oceans > Ocean Topography > Sea Surface Height

  keywords_vocabulary :

  NetCDF COARDS Climate and Forecast Standard Names

  license :

  http://marine.copernicus.eu/web/27-service-commitments-and-licence.php

  platform :

  Altika Drifting Phase

  processing_level :

  L3

  product_version :

  vDec2021

  project :

  COPERNICUS MARINE ENVIRONMENT MONITORING SERVICE (CMEMS)

  references :

  http://marine.copernicus.eu

  software_version :

  7.0_DUACS_DT2021_baseline

  source :

  Altika Drifting Phase measurements

  ssalto_duacs_comment :

  The reference mission used for the altimeter inter-calibration processing is Topex/Poseidon between 1993-01-01 and 2002-04-23, Jason-1 between 2002-04-24 and 2008-10-18, OSTM/Jason-2 between 2008-10-19 and 2016-06-25, Jason-3 since 2016-06-25.

  standard_name_vocabulary :

  NetCDF Climate and Forecast (CF) Metadata Convention Standard Name Table v37

  summary :

  SSALTO/DUACS Delayed-Time Level-3 sea surface height measured by Altika Drifting Phase altimetry observations over Global Ocean.

  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 Along track SSALTO/DUACS Sea Surface Height L3 product

Sea Surface Height maps to evaluate

%%time
list_of_maps = sorted(glob('../data/maps/DUACS_GS.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

CPU times: user 39.1 ms, sys: 6.14 ms, total: 45.3 ms
Wall time: 69.9 ms

<xarray.Dataset>
Dimensions:    (latitude: 44, longitude: 44, time: 365)
Coordinates:
  * latitude   (latitude) float32 32.62 32.88 33.12 33.38 ... 42.88 43.12 43.38
  * longitude  (longitude) float64 294.6 294.9 295.1 295.4 ... 304.9 305.1 305.4
  * time       (time) datetime64[ns] 2019-01-01 2019-01-02 ... 2019-12-31
Data variables:
    sla        (time, latitude, longitude) float64 dask.array<chunksize=(365, 44, 44), meta=np.ndarray>
    ugosa      (time, latitude, longitude) float64 dask.array<chunksize=(365, 44, 44), meta=np.ndarray>
    vgosa      (time, latitude, longitude) float64 dask.array<chunksize=(365, 44, 44), meta=np.ndarray>
    adt        (time, latitude, longitude) float64 dask.array<chunksize=(365, 44, 44), meta=np.ndarray>
    ugos       (time, latitude, longitude) float64 dask.array<chunksize=(365, 44, 44), meta=np.ndarray>
    vgos       (time, latitude, longitude) float64 dask.array<chunksize=(365, 44, 44), meta=np.ndarray>

xarray.Dataset

• Dimensions:
  • latitude: 44
  • longitude: 44
  • time: 365
• Coordinates: (3)
  • latitude
    (latitude)
    float32
    32.62 32.88 33.12 ... 43.12 43.38

    axis :

    Y

    bounds :

    lat_bnds

    long_name :

    Latitude

    standard_name :

    latitude

    units :

    degrees_north

    valid_max :

    89.875

    valid_min :

    -89.875

    array([32.625, 32.875, 33.125, 33.375, 33.625, 33.875, 34.125, 34.375, 34.625,
           34.875, 35.125, 35.375, 35.625, 35.875, 36.125, 36.375, 36.625, 36.875,
           37.125, 37.375, 37.625, 37.875, 38.125, 38.375, 38.625, 38.875, 39.125,
           39.375, 39.625, 39.875, 40.125, 40.375, 40.625, 40.875, 41.125, 41.375,
           41.625, 41.875, 42.125, 42.375, 42.625, 42.875, 43.125, 43.375],
          dtype=float32)

  • longitude
    (longitude)
    float64
    294.6 294.9 295.1 ... 305.1 305.4

    valid_max :

    359.875

    valid_min :

    0.125

    bounds :

    lon_bnds

    axis :

    X

    long_name :

    Longitude

    standard_name :

    longitude

    units :

    degrees_east

    array([294.625, 294.875, 295.125, 295.375, 295.625, 295.875, 296.125, 296.375,
           296.625, 296.875, 297.125, 297.375, 297.625, 297.875, 298.125, 298.375,
           298.625, 298.875, 299.125, 299.375, 299.625, 299.875, 300.125, 300.375,
           300.625, 300.875, 301.125, 301.375, 301.625, 301.875, 302.125, 302.375,
           302.625, 302.875, 303.125, 303.375, 303.625, 303.875, 304.125, 304.375,
           304.625, 304.875, 305.125, 305.375])

  • time
    (time)
    datetime64[ns]
    2019-01-01 ... 2019-12-31

    axis :

    T

    long_name :

    Time

    standard_name :

    time

    array(['2019-01-01T00:00:00.000000000', '2019-01-02T00:00:00.000000000',
           '2019-01-03T00:00:00.000000000', ..., '2019-12-29T00:00:00.000000000',
           '2019-12-30T00:00:00.000000000', '2019-12-31T00:00:00.000000000'],
          dtype='datetime64[ns]')

• Data variables: (6)
  • sla
    (time, latitude, longitude)
    float64
    dask.array<chunksize=(365, 44, 44), meta=np.ndarray>

    ancillary_variables :

    err_sla

    comment :

    The sea level anomaly is the sea surface height above mean sea surface; it is referenced to the [1993, 2012] period; see the product user manual for details

    grid_mapping :

    crs

    long_name :

    Sea level anomaly

    standard_name :

    sea_surface_height_above_sea_level

    units :

    m

    Array Chunk Bytes 5.39 MiB 5.39 MiB Shape (365, 44, 44) (365, 44, 44) Dask graph 1 chunks in 2 graph layers Data type float64 numpy.ndarray

  • ugosa
    (time, latitude, longitude)
    float64
    dask.array<chunksize=(365, 44, 44), meta=np.ndarray>

    ancillary_variables :

    err_ugosa

    comment :

    The geostrophic velocity anomalies are referenced to the [1993, 2012] period

    grid_mapping :

    crs

    long_name :

    Geostrophic velocity anomalies: zonal component

    standard_name :

    surface_geostrophic_eastward_sea_water_velocity_assuming_sea_level_for_geoid

    units :

    m/s

    Array Chunk Bytes 5.39 MiB 5.39 MiB Shape (365, 44, 44) (365, 44, 44) Dask graph 1 chunks in 2 graph layers Data type float64 numpy.ndarray

  • vgosa
    (time, latitude, longitude)
    float64
    dask.array<chunksize=(365, 44, 44), meta=np.ndarray>

    ancillary_variables :

    err_vgosa

    comment :

    The geostrophic velocity anomalies are referenced to the [1993, 2012] period

    grid_mapping :

    crs

    long_name :

    Geostrophic velocity anomalies: meridian component

    standard_name :

    surface_geostrophic_northward_sea_water_velocity_assuming_sea_level_for_geoid

    units :

    m/s

    Array Chunk Bytes 5.39 MiB 5.39 MiB Shape (365, 44, 44) (365, 44, 44) Dask graph 1 chunks in 2 graph layers Data type float64 numpy.ndarray

  • adt
    (time, latitude, longitude)
    float64
    dask.array<chunksize=(365, 44, 44), meta=np.ndarray>

    comment :

    The absolute dynamic topography is the sea surface height above geoid; the adt is obtained as follows: adt=sla+mdt where mdt is the mean dynamic topography; see the product user manual for details

    grid_mapping :

    crs

    long_name :

    Absolute dynamic topography

    standard_name :

    sea_surface_height_above_geoid

    units :

    m

    Array Chunk Bytes 5.39 MiB 5.39 MiB Shape (365, 44, 44) (365, 44, 44) Dask graph 1 chunks in 2 graph layers Data type float64 numpy.ndarray

  • ugos
    (time, latitude, longitude)
    float64
    dask.array<chunksize=(365, 44, 44), meta=np.ndarray>

    grid_mapping :

    crs

    long_name :

    Absolute geostrophic velocity: zonal component

    standard_name :

    surface_geostrophic_eastward_sea_water_velocity

    units :

    m/s

    Array Chunk Bytes 5.39 MiB 5.39 MiB Shape (365, 44, 44) (365, 44, 44) Dask graph 1 chunks in 2 graph layers Data type float64 numpy.ndarray

  • vgos
    (time, latitude, longitude)
    float64
    dask.array<chunksize=(365, 44, 44), meta=np.ndarray>

    grid_mapping :

    crs

    long_name :

    Absolute geostrophic velocity: meridian component

    standard_name :

    surface_geostrophic_northward_sea_water_velocity

    units :

    m/s

    Array Chunk Bytes 5.39 MiB 5.39 MiB Shape (365, 44, 44) (365, 44, 44) Dask graph 1 chunks in 2 graph layers Data type float64 numpy.ndarray

• Indexes: (3)
  • latitude
    PandasIndex

    PandasIndex(Float64Index([32.625, 32.875, 33.125, 33.375, 33.625, 33.875, 34.125, 34.375,
                  34.625, 34.875, 35.125, 35.375, 35.625, 35.875, 36.125, 36.375,
                  36.625, 36.875, 37.125, 37.375, 37.625, 37.875, 38.125, 38.375,
                  38.625, 38.875, 39.125, 39.375, 39.625, 39.875, 40.125, 40.375,
                  40.625, 40.875, 41.125, 41.375, 41.625, 41.875, 42.125, 42.375,
                  42.625, 42.875, 43.125, 43.375],
                 dtype='float64', name='latitude'))

  • longitude
    PandasIndex

    PandasIndex(Float64Index([294.625, 294.875, 295.125, 295.375, 295.625, 295.875, 296.125,
                  296.375, 296.625, 296.875, 297.125, 297.375, 297.625, 297.875,
                  298.125, 298.375, 298.625, 298.875, 299.125, 299.375, 299.625,
                  299.875, 300.125, 300.375, 300.625, 300.875, 301.125, 301.375,
                  301.625, 301.875, 302.125, 302.375, 302.625, 302.875, 303.125,
                  303.375, 303.625, 303.875, 304.125, 304.375, 304.625, 304.875,
                  305.125, 305.375],
                 dtype='float64', name='longitude'))

  • time
    PandasIndex

    PandasIndex(DatetimeIndex(['2019-01-01', '2019-01-02', '2019-01-03', '2019-01-04',
                   '2019-01-05', '2019-01-06', '2019-01-07', '2019-01-08',
                   '2019-01-09', '2019-01-10',
                   ...
                   '2019-12-22', '2019-12-23', '2019-12-24', '2019-12-25',
                   '2019-12-26', '2019-12-27', '2019-12-28', '2019-12-29',
                   '2019-12-30', '2019-12-31'],
                  dtype='datetime64[ns]', name='time', length=365, freq=None))

• Attributes: (0)

2. Statistical & Spectral Analysis

2.1 Interpolate sea surface height maps onto along-track positions

ds_interp = run_interpolation(ds_maps, ds_alg)
ds_interp = ds_interp.dropna('time')
ds_interp

2023-07-04 14:34:05 INFO     fetch data from 2019-01-01 00:00:00 to 2019-02-01 00:00:00
2023-07-04 14:34:05 INFO     fetch data from 2019-01-31 00:00:00 to 2019-03-01 00:00:00
2023-07-04 14:34:05 INFO     fetch data from 2019-02-28 00:00:00 to 2019-04-01 00:00:00
2023-07-04 14:34:05 INFO     fetch data from 2019-03-31 00:00:00 to 2019-05-01 00:00:00
2023-07-04 14:34:05 INFO     fetch data from 2019-04-30 00:00:00 to 2019-06-01 00:00:00
2023-07-04 14:34:05 INFO     fetch data from 2019-05-31 00:00:00 to 2019-07-01 00:00:00
2023-07-04 14:34:05 INFO     fetch data from 2019-06-30 00:00:00 to 2019-08-01 00:00:00
2023-07-04 14:34:05 INFO     fetch data from 2019-07-31 00:00:00 to 2019-09-01 00:00:00
2023-07-04 14:34:05 INFO     fetch data from 2019-08-31 00:00:00 to 2019-10-01 00:00:00
2023-07-04 14:34:05 INFO     fetch data from 2019-09-30 00:00:00 to 2019-11-01 00:00:00
2023-07-04 14:34:05 INFO     fetch data from 2019-10-31 00:00:00 to 2019-12-01 00:00:00
2023-07-04 14:34:05 INFO     fetch data from 2019-11-30 00:00:00 to 2019-12-31 00:00:00

<xarray.Dataset>
Dimensions:            (time: 37073)
Coordinates:
  * time               (time) datetime64[ns] 2019-01-01T22:10:01.073445888 .....
Data variables: (12/13)
    cycle              (time) float64 126.0 126.0 126.0 ... 136.0 136.0 136.0
    track              (time) float64 36.0 36.0 36.0 36.0 ... 391.0 391.0 391.0
    sla_unfiltered     (time) float32 0.124 0.177 0.163 ... 0.152 0.181 0.176
    sla_filtered       (time) float32 0.151 0.161 0.166 ... 0.151 0.169 0.178
    dac                (time) float32 0.058 0.059 0.059 ... -0.059 -0.059 -0.06
    ocean_tide         (time) float32 0.172 0.172 0.173 ... -0.339 -0.342 -0.349
    ...                 ...
    lwe                (time) float32 -0.044 -0.043 -0.043 ... -0.047 -0.047
    mdt                (time) float32 -0.174 -0.172 -0.17 ... -0.173 -0.17
    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 305.0 304.9 304.9 ... 297.6 297.6 297.6
    latitude           (time) float64 42.99 42.93 42.87 ... 42.86 42.92 42.98
    msla_interpolated  (time) float64 0.1328 0.1322 0.1284 ... 0.106 0.1141

xarray.Dataset

• Dimensions:
  • time: 37073
• Coordinates: (1)
  • time
    (time)
    datetime64[ns]
    2019-01-01T22:10:01.073445888 .....

    array(['2019-01-01T22:10:01.073445888', '2019-01-01T22:10:02.129446144',
           '2019-01-01T22:10:03.185446144', ..., '2019-12-30T09:33:55.091854080',
           '2019-12-30T09:33:56.147854080', '2019-12-30T09:33:57.203853824'],
          dtype='datetime64[ns]')

• Data variables: (13)
  • cycle
    (time)
    float64
    126.0 126.0 126.0 ... 136.0 136.0

    array([126., 126., 126., ..., 136., 136., 136.])

  • track
    (time)
    float64
    36.0 36.0 36.0 ... 391.0 391.0

    array([ 36.,  36.,  36., ..., 391., 391., 391.])

  • sla_unfiltered
    (time)
    float32
    0.124 0.177 0.163 ... 0.181 0.176

    array([0.12400001, 0.177     , 0.163     , ..., 0.15200001, 0.18100001,
           0.17600001], dtype=float32)

  • sla_filtered
    (time)
    float32
    0.151 0.161 0.166 ... 0.169 0.178

    array([0.15100001, 0.16100001, 0.16600001, ..., 0.15100001, 0.16900001,
           0.178     ], dtype=float32)

  • dac
    (time)
    float32
    0.058 0.059 0.059 ... -0.059 -0.06

    array([ 0.058,  0.059,  0.059, ..., -0.059, -0.059, -0.06 ], dtype=float32)

  • ocean_tide
    (time)
    float32
    0.172 0.172 0.173 ... -0.342 -0.349

    array([ 0.172     ,  0.172     ,  0.17300001, ..., -0.33900002,
           -0.342     , -0.349     ], dtype=float32)

  • internal_tide
    (time)
    float32
    0.001 -0.001 -0.002 ... 0.004 0.003

    array([ 0.001, -0.001, -0.002, ...,  0.005,  0.004,  0.003], dtype=float32)

  • lwe
    (time)
    float32
    -0.044 -0.043 ... -0.047 -0.047

    array([-0.044, -0.043, -0.043, ..., -0.047, -0.047, -0.047], dtype=float32)

  • mdt
    (time)
    float32
    -0.174 -0.172 ... -0.173 -0.17

    array([-0.17400001, -0.172     , -0.17      , ..., -0.17500001,
           -0.17300001, -0.17      ], dtype=float32)

  • tpa_correction
    (time)
    float32
    0.0 0.0 0.0 0.0 ... 0.0 0.0 0.0 0.0

    array([0., 0., 0., ..., 0., 0., 0.], dtype=float32)

  • longitude
    (time)
    float64
    305.0 304.9 304.9 ... 297.6 297.6

    array([304.959191, 304.937747, 304.916336, ..., 297.643509, 297.622114,
           297.600686])

  • latitude
    (time)
    float64
    42.99 42.93 42.87 ... 42.92 42.98

    array([42.989853, 42.92905 , 42.868244, ..., 42.862757, 42.923512,
           42.984265])

  • msla_interpolated
    (time)
    float64
    0.1328 0.1322 ... 0.106 0.1141

    array([0.13280733, 0.13222027, 0.12840952, ..., 0.10332474, 0.10596411,
           0.11408924])

• Indexes: (1)
  • time
    PandasIndex

    PandasIndex(DatetimeIndex(['2019-01-01 22:10:01.073445888',
                   '2019-01-01 22:10:02.129446144',
                   '2019-01-01 22:10:03.185446144',
                   '2019-01-01 22:10:04.241445888',
                   '2019-01-01 22:10:05.297445888',
                   '2019-01-01 22:10:06.353445888',
                   '2019-01-01 22:10:07.409445888',
                   '2019-01-01 22:10:08.465445888',
                   '2019-01-01 22:10:09.521446144',
                   '2019-01-01 22:10:10.577446144',
                   ...
                   '2019-12-30 09:33:47.699854080',
                   '2019-12-30 09:33:48.755854080',
                   '2019-12-30 09:33:49.811853824',
                   '2019-12-30 09:33:50.867853824',
                   '2019-12-30 09:33:51.923853824',
                   '2019-12-30 09:33:52.979853824',
                   '2019-12-30 09:33:54.035853824',
                   '2019-12-30 09:33:55.091854080',
                   '2019-12-30 09:33:56.147854080',
                   '2019-12-30 09:33:57.203853824'],
                  dtype='datetime64[ns]', name='time', length=37073, freq=None))

• Attributes: (0)

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

Once the maps have been interpolated to the position of the along-track, it is possible to calculate different statistics on the time series SLA_alongtrack and SLA_maps.

We propose below the following statistics: error variance maps (static by 1°x1° box), explained variance maps, temporal evolution of global error variance and explained variance.

These statistics are also applied to the filtered signals focusing on the 65-500km scale range. A bandpass filter is applied before calculating the scores.

compute_stat_scores(ds_interp, lambda_min, lambda_max, stat_output_filename,method_name)

2023-07-04 14:34:05 INFO     Compute mapping error all scales
2023-07-04 14:34:05 INFO     Compute mapping error for scales between 65.0 and 500.0 km
2023-07-04 14:34:06 INFO     Compute binning statistics
2023-07-04 14:34:07 INFO     Compute statistics by oceanic regime
2023-07-04 14:34:31 INFO     Stat file saved as: ../results/stat_sla_duacs_geos_GS.nc

# Plot gridded stats
# Hvplot
# plot_stat_score_map(stat_output_filename)
# Matplotlib
plot_stat_score_map_png(stat_output_filename,region=region,box_lonlat=box_lonlat_GS)

plot_stat_score_timeseries(stat_output_filename)

plot_stat_by_regimes(stat_output_filename)

	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.001673	0.012362	0.000474	0.008133	0.864663	0.941740
offshore_highvar	0.004331	0.075310	0.001669	0.034653	0.942486	0.951841
offshore_lowvar	0.001736	0.017152	0.000561	0.010251	0.898770	0.945300
equatorial_band	NaN	NaN	NaN	NaN	NaN	NaN
arctic	NaN	NaN	NaN	NaN	NaN	NaN
antarctic	NaN	NaN	NaN	NaN	NaN	NaN

2.3 Compute Spectral scores

compute_psd_scores_v2(ds_interp, psd_output_filename, lenght_scale=segment_lenght,method_name=method_name)

2023-07-04 14:34:38 INFO     Segment computation...
2023-07-04 14:34:38 INFO     Spectral analysis...
2023-07-04 14:34:41 INFO     Saving ouput...
2023-07-04 14:34:45 INFO     PSD file saved as: ../results/psd_sla_duacs_geos_GS.nc

# Plot effective resolution
# Hvplot
# plot_effective_resolution(psd_output_filename)
# Matplotlib
plot_effective_resolution_png(psd_output_filename,region=region,box_lonlat=box_lonlat_GS)

plot_psd_scores(psd_output_filename)