Benchmark DUACS Lagrangian Cumulative Distance

2023-04-27 DUACS_SSH_BENCHMARK_DEMO

---

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

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 trajectories of independent drifters distributed by CMEMS </h5>

1. Import packages

import os
import numpy as np
import xarray as xr
import matplotlib
import matplotlib.pylab as plt
import matplotlib.ticker as mticker
import cmocean
from scipy.interpolate import RegularGridInterpolator
import pickle
import gc
from matplotlib import cm
import warnings
warnings.filterwarnings('ignore')

%matplotlib inline

import sys
sys.path.append('..')
from src.mod_plot import *
from src.mod_traj import *
from src.mod_compare import regional_zoom

2. Setup parameters

method_name = 'DUACS_global'

time_min = '2019-01-01'                                        # time min for analysis
time_max = '2019-12-31'                                        # time max for analysis

dir_out = f'../results/'                               # output directory path
prefix_out = 'dict_drifter_adv'
results_out = 'deviat_uv_'+method_name+'.nc'

if not os.path.exists(dir_out):
    os.system('mkdir '+dir_out)

# Global
lon_min = 0                                          # domain min longitude
lon_max = 360                                          # domain max longitude
lat_min = -90.                                          # domain min latitude
lat_max = 90.                                          # domain max latitude

3. Maps to evaluate

3.1 Download maps

path_maps =  f'../data/maps/DUACS_global_allsat-alg/dt*'

maps = xr.open_mfdataset(path_maps)

maps = maps.sel({'longitude':slice(lon_min-0.5,lon_max+0.5)})
maps = maps.sel({'latitude':slice(lat_min-0.5,lat_max+0.5)})
maps

<xarray.Dataset>
Dimensions:    (latitude: 720, longitude: 1440, time: 365)
Coordinates:
  * latitude   (latitude) float32 -89.88 -89.62 -89.38 ... 89.38 89.62 89.88
  * longitude  (longitude) float64 0.125 0.375 0.625 0.875 ... 359.4 359.6 359.9
  * time       (time) datetime64[ns] 2019-01-01 2019-01-02 ... 2019-12-31
Data variables:
    sla        (time, latitude, longitude) float64 dask.array<chunksize=(1, 720, 1440), meta=np.ndarray>
    ugosa      (time, latitude, longitude) float64 dask.array<chunksize=(1, 720, 1440), meta=np.ndarray>
    vgosa      (time, latitude, longitude) float64 dask.array<chunksize=(1, 720, 1440), meta=np.ndarray>
    adt        (time, latitude, longitude) float64 dask.array<chunksize=(1, 720, 1440), meta=np.ndarray>
    ugos       (time, latitude, longitude) float64 dask.array<chunksize=(1, 720, 1440), meta=np.ndarray>
    vgos       (time, latitude, longitude) float64 dask.array<chunksize=(1, 720, 1440), meta=np.ndarray>

xarray.Dataset

• Dimensions:
  • latitude: 720
  • longitude: 1440
  • time: 365
• Coordinates: (3)
  • latitude
    (latitude)
    float32
    -89.88 -89.62 ... 89.62 89.88

    axis :

    Y

    bounds :

    lat_bnds

    long_name :

    Latitude

    standard_name :

    latitude

    units :

    degrees_north

    valid_max :

    89.875

    valid_min :

    -89.875

    array([-89.875, -89.625, -89.375, ...,  89.375,  89.625,  89.875],
          dtype=float32)

  • longitude
    (longitude)
    float64
    0.125 0.375 0.625 ... 359.6 359.9

    valid_max :

    359.875

    valid_min :

    0.125

    bounds :

    lon_bnds

    axis :

    X

    long_name :

    Longitude

    standard_name :

    longitude

    units :

    degrees_east

    array([1.25000e-01, 3.75000e-01, 6.25000e-01, ..., 3.59375e+02, 3.59625e+02,
           3.59875e+02])

  • 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=(1, 720, 1440), 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 2.82 GiB 7.91 MiB Shape (365, 720, 1440) (1, 720, 1440) Dask graph 365 chunks in 731 graph layers Data type float64 numpy.ndarray

  • ugosa
    (time, latitude, longitude)
    float64
    dask.array<chunksize=(1, 720, 1440), 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 2.82 GiB 7.91 MiB Shape (365, 720, 1440) (1, 720, 1440) Dask graph 365 chunks in 731 graph layers Data type float64 numpy.ndarray

  • vgosa
    (time, latitude, longitude)
    float64
    dask.array<chunksize=(1, 720, 1440), 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 2.82 GiB 7.91 MiB Shape (365, 720, 1440) (1, 720, 1440) Dask graph 365 chunks in 731 graph layers Data type float64 numpy.ndarray

  • adt
    (time, latitude, longitude)
    float64
    dask.array<chunksize=(1, 720, 1440), 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 2.82 GiB 7.91 MiB Shape (365, 720, 1440) (1, 720, 1440) Dask graph 365 chunks in 731 graph layers Data type float64 numpy.ndarray

  • ugos
    (time, latitude, longitude)
    float64
    dask.array<chunksize=(1, 720, 1440), 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 2.82 GiB 7.91 MiB Shape (365, 720, 1440) (1, 720, 1440) Dask graph 365 chunks in 731 graph layers Data type float64 numpy.ndarray

  • vgos
    (time, latitude, longitude)
    float64
    dask.array<chunksize=(1, 720, 1440), 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 2.82 GiB 7.91 MiB Shape (365, 720, 1440) (1, 720, 1440) Dask graph 365 chunks in 731 graph layers Data type float64 numpy.ndarray

• Indexes: (3)
  • latitude
    PandasIndex

    PandasIndex(Float64Index([-89.875, -89.625, -89.375, -89.125, -88.875, -88.625, -88.375,
                  -88.125, -87.875, -87.625,
                  ...
                   87.625,  87.875,  88.125,  88.375,  88.625,  88.875,  89.125,
                   89.375,  89.625,  89.875],
                 dtype='float64', name='latitude', length=720))

  • longitude
    PandasIndex

    PandasIndex(Float64Index([  0.125,   0.375,   0.625,   0.875,   1.125,   1.375,   1.625,
                    1.875,   2.125,   2.375,
                  ...
                  357.625, 357.875, 358.125, 358.375, 358.625, 358.875, 359.125,
                  359.375, 359.625, 359.875],
                 dtype='float64', name='longitude', length=1440))

  • 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)

3.2 Retrieve maps info

time = maps.time.values
lon = maps.longitude.values
lat = maps.latitude.values

time = (time.astype('datetime64[h]') - np.datetime64('2019-01-01')).astype(int)
Nt = time[-1]

u_maps = maps.ugos
v_maps = maps.vgos

4. Independant drifters

4.1 Download drifters

path_drifters = '../data/independent_drifters/uv*'
ds_drifters = xr.open_mfdataset(path_drifters,concat_dim='time',combine='nested')
ds_drifters

<xarray.Dataset>
Dimensions:    (time: 2156405)
Coordinates:
  * time       (time) datetime64[ns] 2019-01-01 ... 2019-12-31
    latitude   (time) float32 dask.array<chunksize=(5559,), meta=np.ndarray>
    longitude  (time) float32 dask.array<chunksize=(5559,), meta=np.ndarray>
Data variables:
    EWCT       (time) float32 dask.array<chunksize=(5559,), meta=np.ndarray>
    NSCT       (time) float32 dask.array<chunksize=(5559,), meta=np.ndarray>
    sensor_id  (time) float64 dask.array<chunksize=(5559,), meta=np.ndarray>
Attributes: (12/46)
    data_type:                   OceanSITES trajectory data
    format_version:              2.0
    platform_code:               116275
    date_update:                 2020-10-13T12:17:40Z
    institution:                 AOML
    institution_edmo_code:       1799
    ...                          ...
    deployment_lat:              -58.44
    last_longitude_observation:  82.75
    last_latitude_observation:   -18.49
    date_drog_lost:              2017-01-21T03:37:00Z
    death_type:                  stop transmitting
    last_date_observation:       2019-01-16T01:51:00Z

xarray.Dataset

• Dimensions:
  • time: 2156405
• Coordinates: (3)
  • time
    (time)
    datetime64[ns]
    2019-01-01 ... 2019-12-31

    long_name :

    Time

    standard_name :

    time

    valid_min :

    0.0

    valid_max :

    90000.0

    QC_indicator :

    1.0

    QC_procedure :

    1.0

    uncertainty :

    comment :

    axis :

    T

    array(['2019-01-01T00:00:00.000000000', '2019-01-01T06:00:00.000000000',
           '2019-01-01T12:00:00.000000000', ..., '2019-12-31T00:00:00.000000000',
           '2019-12-31T00:00:00.000000000', '2019-12-31T00:00:00.000000000'],
          dtype='datetime64[ns]')

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

    Array Chunk Bytes 8.23 MiB 26.33 kiB Shape (2156405,) (6741,) Dask graph 365 chunks in 731 graph layers Data type float32 numpy.ndarray

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

    Array Chunk Bytes 8.23 MiB 26.33 kiB Shape (2156405,) (6741,) Dask graph 365 chunks in 731 graph layers Data type float32 numpy.ndarray

• Data variables: (3)
  • EWCT
    (time)
    float32
    dask.array<chunksize=(5559,), meta=np.ndarray>

    long_name :

    West-east current component at the drog depht

    standard_name :

    eastward_sea_water_velocity

    units :

    m s-1

    Array Chunk Bytes 8.23 MiB 26.33 kiB Shape (2156405,) (6741,) Dask graph 365 chunks in 731 graph layers Data type float32 numpy.ndarray

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

    long_name :

    South-north current component at the drog depht

    standard_name :

    northward_sea_water_velocity

    units :

    m s-1

    Array Chunk Bytes 8.23 MiB 26.33 kiB Shape (2156405,) (6741,) Dask graph 365 chunks in 731 graph layers Data type float32 numpy.ndarray

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

    Array Chunk Bytes 16.45 MiB 52.66 kiB Shape (2156405,) (6741,) Dask graph 365 chunks in 731 graph layers Data type float64 numpy.ndarray

• Indexes: (1)
  • time
    PandasIndex

    PandasIndex(DatetimeIndex(['2019-01-01 00:00:00', '2019-01-01 06:00:00',
                   '2019-01-01 12:00:00', '2019-01-01 18:00:00',
                   '2019-01-01 00:00:00', '2019-01-01 06:00:00',
                   '2019-01-01 12:00:00', '2019-01-01 18:00:00',
                   '2019-01-01 00:00:00', '2019-01-01 06:00:00',
                   ...
                   '2019-12-31 00:00:00', '2019-12-31 00:00:00',
                   '2019-12-31 00:00:00', '2019-12-31 00:00:00',
                   '2019-12-31 00:00:00', '2019-12-31 00:00:00',
                   '2019-12-31 00:00:00', '2019-12-31 00:00:00',
                   '2019-12-31 00:00:00', '2019-12-31 00:00:00'],
                  dtype='datetime64[ns]', name='time', length=2156405, freq=None))

• Attributes: (46)

  data_type :

  OceanSITES trajectory data

  format_version :

  2.0

  platform_code :

  116275

  date_update :

  2020-10-13T12:17:40Z

  institution :

  AOML

  institution_edmo_code :

  1799

  wmo_platform_code :

  3300883

  platform_name :

  DRIFTING BUOY

  history :

  data_mode :

  D

  references :

  Product User Manual http://dx.doi.org/10.13155/41257; Quality Information Document http://dx.doi.org/10.13155/41256; http://marine.copernicus.eu; http://www.coriolis.eu.org

  comment :

  5254 SVPB Name: BUOY_TYPE, dtype: object: Standard Surface Velocity Program drifter with Barometer

  Conventions :

  CF-1.6 OceanSITES-Manual-1.2 Copernicus-InSituTAC-SRD-1.3 Copernicus-InSituTAC-ParametersList-3.0.0

  netcdf_version :

  netCDF-4 classic model

  title :

  Global Ocean- Delayed Mode in-situ observations of ocean surface currents

  summary :

  In-situ observation yearly delivery in delayed mode of Ocean surface currents. The In Situ delayed mode product designed for reanalysis purposes integrates the best available version of in situ data for Ocean surface currents. The data are collected from the Surface Drifter Data Assembly Centre (SD-DAC at NOAA AOML). All surface drifters data have been processed to check for drogue loss. Drogued and undrogued drifting buoy surface ocean currents are provided with a drogue presence flag as well as a wind slippage correction for undrogued buoy. The product is designed to be assimilated into or for validation purposes of operational models operated by ocean forecasting centers for reanalysis purposes or for research community.

  naming_authority :

  OceanSITES, Copernicus Marine In Situ TAC

  citation :

  These data were collected and made freely available by the Copernicus project and the programs that contribute to it

  doi :

  10.17882/41334

  id :

  cdm_data_type :

  Time-series

  area :

  Global Ocean

  geospatial_lat_min :

  -58.379

  geospatial_lat_max :

  -18.479

  geospatial_lon_min :

  -62.86400000000003

  geospatial_lon_max :

  104.132

  geospatial_vertical_min :

  0.5

  geospatial_vertical_max :

  15

  time_coverage_start :

  2016-01-06T12:00:00Z

  time_coverage_end :

  2019-01-16T00:00:00Z

  contact :

  cmems-service@ifremer.fr

  author :

  Copernicus In Situ Thematic Assembly Center

  data_assembly_center :

  AOML

  pi_name :

  Helene Etienne : hetienne@groupcls.com

  distribution_statement :

  These data follow Copernicus standards; they are public and free of charge. User assumes all risk for use of data. User must display citation in any publication or product using data. User must contact PI prior to any commercial use of data.

  update_interval :

  biannual

  qc_manual :

  OceanSITES Users Manual v1.2

  AOML_experiment_number :

  6129

  deployment_date :

  2016-01-06T07:11:00Z

  deployment_lon :

  -63.02999999999997

  deployment_lat :

  -58.44

  last_longitude_observation :

  82.75

  last_latitude_observation :

  -18.49

  date_drog_lost :

  2017-01-21T03:37:00Z

  death_type :

  stop transmitting

  last_date_observation :

  2019-01-16T01:51:00Z

4.2 Prepare drifter data

%%time
ind, time_drifter, lon_drifter, lat_drifter, id_drifter= prepare_drifter_data(ds_drifters,maps)

CPU times: user 1min 54s, sys: 15.9 s, total: 2min 10s
Wall time: 2min

5. Interpolate maps at drifters locations

%%time

# Create interpolation function
fu = RegularGridInterpolator((time, lat, lon), u_maps.values)
fv = RegularGridInterpolator((time, lat, lon), v_maps.values)

CPU times: user 30.4 s, sys: 12.4 s, total: 42.8 s
Wall time: 49.7 s

6. Artificial drifter trajectories

6.1 Trajectory parameters and initialization

nday = 5
dt_h = 24
mode = 'euler'
write = False
Np = time_drifter.size

horizons = np.arange(0,nday*24/dt_h+1)

dict_drifter_adv_maps = {0:{'time':time_drifter,
                           'lon':lon_drifter,
                           'lat':lat_drifter,
                            'id':id_drifter}
                       }
horizons

array([0., 1., 2., 3., 4., 5.])

6.2 Compute trajectories

%%time
dict_drifter_adv_maps = compute_traj(dict_drifter_adv_maps, horizons, fu, fv ,Np, Nt, dt_h, method_name, dir_out, prefix_out, mode, write)

h 1.0
h 2.0.0
h 3.0.0
h 4.0.0
h 5.0.0
CPU times: user 143 ms, sys: 644 ms, total: 787 ms
Wall time: 1.23 s

7. Compute deviation

# Deviation = Cumulative distance between artificial and real trajectories

# Average deviation by bins
lon_out=np.arange(lon_min, lon_max, 1)
lat_out=np.arange(lat_min, lat_max, 1)

horizon_days = [1,2,3,4,5]

dev_maps, var_dev_maps, dmean = compute_deviation(dict_drifter_adv_maps, Nt, dt_h, horizon_days,lon_out,lat_out, dir_out, method_name, results_out)

8. Plots

8.1 Plot deviation maps at different horizons

horizon_plots = [1,2,3,4,5]
horizon_plot_max = [40,80,125,150,200]

plot_traj_deviation_maps( lon_out, lat_out, dev_maps, dt_h,horizon_plots,horizon_plot_max,method_name)

8.2 Plot spatially averaged deviation as a function of horizons

plot_meantraj_deviation(dir_out)

plot_meantraj_deviation(dir_out, just_basin='global')