Gulf Stream: Benchmark DUACS geostrophic currents maps

2023-04-27 DUACS_SSH_BENCHMARK_DEMO

---

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

Gulf Stream: Benchmark of DUACS geostrophic currents maps

The notebook aims to evaluate the surface current 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 current data from the drifters database distributed by CMEMS (INSITU_GLO_PHY_UV_DISCRETE_MY_013_044 product) </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 currents maps with drifters database: statistical and spectral analysis

from glob import glob
import numpy as np
import os
import warnings
warnings.filterwarnings("ignore")

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_uv_duacs_geos_GS.nc'   # output statistical analysis filename
psd_output_filename = f'{output_dir}/psd_uv_duacs_geos_GS.nc'     # output spectral analysis filename
segment_lenght = np.timedelta64(40, 'D')                      # spectral parameter: drifters segment lenght in days to consider in the spectral analysis

1. Input files

Sea Surface currents from Drifters database

filenames_drifters = sorted(glob('../data/independent_drifters/indep_drifters_GS.nc'))

ds_drifter = xr.open_mfdataset(filenames_drifters, combine='nested', concat_dim='time')
ds_drifter = ds_drifter.where((ds_drifter.time >= np.datetime64(time_min)) & (ds_drifter.time <=  np.datetime64(time_max)), drop=True)
ds_drifter

<xarray.Dataset>
Dimensions:    (time: 11235)
Coordinates:
  * time       (time) datetime64[ns] 2019-01-01 2019-01-01 ... 2019-12-31
    latitude   (time) float32 dask.array<chunksize=(11235,), meta=np.ndarray>
    longitude  (time) float32 dask.array<chunksize=(11235,), meta=np.ndarray>
Data variables:
    EWCT       (time) float32 dask.array<chunksize=(11235,), meta=np.ndarray>
    NSCT       (time) float32 dask.array<chunksize=(11235,), meta=np.ndarray>
    sensor_id  (time) float64 dask.array<chunksize=(11235,), 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: 11235
• 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-01T00:00:00.000000000',
           '2019-01-01T00: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=(11235,), meta=np.ndarray>

    Array Chunk Bytes 43.89 kiB 43.89 kiB Shape (11235,) (11235,) Dask graph 1 chunks in 2 graph layers Data type float32 numpy.ndarray

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

    Array Chunk Bytes 43.89 kiB 43.89 kiB Shape (11235,) (11235,) Dask graph 1 chunks in 2 graph layers Data type float32 numpy.ndarray

• Data variables: (3)
  • EWCT
    (time)
    float32
    dask.array<chunksize=(11235,), 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 43.89 kiB 43.89 kiB Shape (11235,) (11235,) Dask graph 1 chunks in 4 graph layers Data type float32 numpy.ndarray

  • NSCT
    (time)
    float32
    dask.array<chunksize=(11235,), 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 43.89 kiB 43.89 kiB Shape (11235,) (11235,) Dask graph 1 chunks in 4 graph layers Data type float32 numpy.ndarray

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

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

• Indexes: (1)
  • time
    PandasIndex

    PandasIndex(DatetimeIndex(['2019-01-01 00:00:00', '2019-01-01 00:00:00',
                   '2019-01-01 00:00:00', '2019-01-01 00:00:00',
                   '2019-01-01 00:00:00', '2019-01-01 00:00:00',
                   '2019-01-01 00:00:00', '2019-01-01 00:00:00',
                   '2019-01-01 00:00:00', '2019-01-01 00:00:00',
                   ...
                   '2019-12-30 18:00:00', '2019-12-30 18:00:00',
                   '2019-12-30 18:00:00', '2019-12-30 18:00:00',
                   '2019-12-30 18: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=11235, 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

Sea Surface current maps to evaluate

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

<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 currents maps onto drifters positions

ds_interp = run_interpolation_drifters(ds_maps, ds_drifter, time_min, time_max)
ds_interp = ds_interp.dropna('time')
ds_interp = ds_interp.sortby('time')
ds_interp

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

<xarray.Dataset>
Dimensions:            (time: 11230)
Coordinates:
  * time               (time) datetime64[ns] 2019-01-01 ... 2019-12-31
Data variables:
    EWCT               (time) float32 -0.1839 -0.2353 0.2235 ... -0.05235 -0.391
    NSCT               (time) float32 -0.3496 -0.2228 ... -0.3979 -0.4466
    sensor_id          (time) float64 6.188e+07 6.322e+07 ... 6.548e+07 6.64e+07
    latitude           (time) float32 35.95 36.12 36.39 ... 35.21 33.88 37.28
    longitude          (time) float32 -55.56 -57.22 -56.21 ... -56.45 -61.2
    ugos_interpolated  (time) float64 -0.2183 -0.2837 0.2636 ... -0.2945 -0.5087
    vgos_interpolated  (time) float64 -0.2328 -0.09264 ... -0.209 -0.6738

xarray.Dataset

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

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

• Data variables: (7)
  • EWCT
    (time)
    float32
    -0.1839 -0.2353 ... -0.05235 -0.391

    array([-0.18394, -0.23534,  0.22352, ...,  0.01319, -0.05235, -0.39097],
          dtype=float32)

  • NSCT
    (time)
    float32
    -0.3496 -0.2228 ... -0.3979 -0.4466

    array([-0.34964, -0.22283, -0.4015 , ..., -0.23825, -0.39785, -0.44663],
          dtype=float32)

  • sensor_id
    (time)
    float64
    6.188e+07 6.322e+07 ... 6.64e+07

    array([61879850., 63221590., 63225600., ..., 63940890., 65480940.,
           66398370.])

  • latitude
    (time)
    float32
    35.95 36.12 36.39 ... 33.88 37.28

    array([35.953, 36.121, 36.386, ..., 35.212, 33.875, 37.284], dtype=float32)

  • longitude
    (time)
    float32
    -55.56 -57.22 ... -56.45 -61.2

    array([-55.559   , -57.218994, -56.208008, ..., -61.596985, -56.44699 ,
           -61.195007], dtype=float32)

  • ugos_interpolated
    (time)
    float64
    -0.2183 -0.2837 ... -0.2945 -0.5087

    array([-0.21833286, -0.28365047,  0.26361407, ..., -0.08105203,
           -0.29453977, -0.50869766])

  • vgos_interpolated
    (time)
    float64
    -0.2328 -0.09264 ... -0.209 -0.6738

    array([-0.23281991, -0.0926385 , -0.149127  , ..., -0.21333821,
           -0.20904642, -0.67383731])

• Indexes: (1)
  • time
    PandasIndex

    PandasIndex(DatetimeIndex(['2019-01-01 00:00:00', '2019-01-01 00:00:00',
                   '2019-01-01 00:00:00', '2019-01-01 00:00:00',
                   '2019-01-01 00:00:00', '2019-01-01 00:00:00',
                   '2019-01-01 00:00:00', '2019-01-01 00:00:00',
                   '2019-01-01 00:00:00', '2019-01-01 00:00:00',
                   ...
                   '2019-12-30 18:00:00', '2019-12-30 18:00:00',
                   '2019-12-30 18:00:00', '2019-12-30 18:00:00',
                   '2019-12-30 18: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=11230, freq=None))

• Attributes: (0)

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

Once the surface currents maps have been interpolated to the position of the drifters, it is possible to calculate different statistics on the time series of zonal and meridional velocities.

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

# Compute gridded stats
compute_stat_scores_uv(ds_interp, stat_output_filename,method_name=method_name)

2023-07-05 16:56:38 INFO     Compute mapping error all scales
2023-07-05 16:56:38 INFO     Compute statistics
2023-07-05 16:56:41 INFO     Stat file saved as: ../results/stat_uv_duacs_geos_GS.nc
2023-07-05 16:56:41 INFO     Compute statistics by oceanic regime

# Plot gridded stats
# Hvplot
# plot_stat_score_map_uv(stat_output_filename)
# Matplotlib
plot_stat_score_map_uv_png(stat_output_filename,region=region,box_lonlat=box_lonlat_GS)

The figure shows that the maximum mapping errors are found in intense current systems, for example in the GulfStream, Kuroshio and Agulhas regions.

However, when considering the full scale of motion in the drifter database, the surface current maps capture up to 80% of the variability of drifter currents in the Western Boundary Currents and Antarctic Circumpolar Currents (ACC). The geostrophic signal dominates the ageostrophic signal in these regions. In regions with low ocean variability, only a few percent of the total drifter current variability is recovered in the maps, which may be associated with a larger ageostrophic signal in these regions.

plot_stat_uv_by_regimes(stat_output_filename)

	mapping_err_u_var [m²/s²]	mapping_err_v_var [m²/s²]	ugos_interpolated_var [m²/s²]	EWCT_var [m²/s²]	vgos_interpolated_var [m²/s²]	NSCT_var [m²/s²]	var_score_u_allscale	var_score_v_allscale
coastal	0.026696	0.028006	0.036708	0.057370	0.052822	0.075548	0.534664	0.629296
offshore_highvar	0.058719	0.056771	0.165690	0.201561	0.132635	0.167143	0.708681	0.660345
offshore_lowvar	0.025558	0.024275	0.033281	0.058642	0.027226	0.053325	0.564163	0.544767
equatorial_band	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN
arctic	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN
antarctic	NaN	NaN	NaN	NaN	NaN	NaN	NaN	NaN

2.4 Compute Spectral scores

# Compute PSD scores
compute_psd_scores_current(ds_interp, psd_output_filename, lenght_scale=segment_lenght, method_name=method_name)

2023-07-05 16:57:18 INFO     Segment computation...
2023-07-05 16:57:19 INFO     Spectral analysis...
2023-07-05 16:57:19 INFO     Write output...
2023-07-05 16:57:19 INFO     PSD file saved as: ../results/psd_uv_duacs_geos_GS.nc

# Plot Zonally averaged rotary spectra
# Hvplot
# plot_psd_scores_currents(psd_output_filename)
# Matplotlib
plot_psd_scores_currents_png(psd_output_filename,region=region)

# Plot Zonally averaged rotary spectra
plot_psd_scores_currents_1D(psd_output_filename)

The interactive plot above allows you to explore the spectral metrics by latitude band