AlongTrack Data

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
import matplotlib.pyplot as plt
import matplotlib.colors as colors
import matplotlib.ticker as ticker
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
import seaborn as sns

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

%load_ext autoreload
%autoreload 2

The autoreload extension is already loaded. To reload it, use:
  %reload_ext autoreload

Data Challenge - 2020a - OSSE

This data challenge looks at some simulated Altimetry data from simulations

!ls "/gpfswork/rech/yrf/commun/data_challenges/dc20a_osse/raw/dc_obs/"

2020a_SSH_mapping_NATL60_envisat.nc
2020a_SSH_mapping_NATL60_geosat2.nc
2020a_SSH_mapping_NATL60_jason1.nc
2020a_SSH_mapping_NATL60_karin_swot.nc
2020a_SSH_mapping_NATL60_nadir_swot.nc
2020a_SSH_mapping_NATL60_topex-poseidon_interleaved.nc

files_nadir_dc20a = [
    "/gpfswork/rech/yrf/commun/data_challenges/dc20a_osse/raw/dc_obs/2020a_SSH_mapping_NATL60_jason1.nc",
    "/gpfswork/rech/yrf/commun/data_challenges/dc20a_osse/raw/dc_obs/2020a_SSH_mapping_NATL60_envisat.nc",
    "/gpfswork/rech/yrf/commun/data_challenges/dc20a_osse/raw/dc_obs/2020a_SSH_mapping_NATL60_geosat2.nc",
    "/gpfswork/rech/yrf/commun/data_challenges/dc20a_osse/raw/dc_obs/2020a_SSH_mapping_NATL60_topex-poseidon_interleaved.nc"
]

def preprocess_nadir_dc20a(da):
    
    # reorganized
    da = da.sortby("time").compute()
    
    # validate coordinates
    da = validate_latlon(da)
    da = validate_time(da)
    
    # validate variables
    da = da.rename({"ssh_model": "ssh"})
    da = validate_ssh(da)
    
    # # slice time period
    # da = da.sel(time=slice("2012-10-22", "2012-12-03"))
    
    # # slice region
    # da = where_slice(da, "lon", -64.975, -55.007)
    # da = where_slice(da, "lat", 33.025, 42.9917)
    

    
    # assign coordinates
    da = da.set_coords(["time", "lat", "lon"])
    
    return da[["ssh"]]

ds_nadir = xr.open_mfdataset(
    files_nadir_dc20a, 
    preprocess=preprocess_nadir_dc20a,
    combine="nested",
    engine="netcdf4",
    concat_dim="time"
)

ds_nadir = ds_nadir.sortby("time").compute()

ds_nadir

<xarray.Dataset>
Dimensions:  (time: 205232)
Coordinates:
  * time     (time) datetime64[ns] 2012-10-01T00:45:50.676369 ... 2013-09-30T...
    lon      (time) float64 -57.81 -57.83 -57.85 -57.87 ... -58.29 -58.31 -58.33
    lat      (time) float64 42.96 42.9 42.84 42.78 ... 42.82 42.88 42.94 43.0
Data variables:
    ssh      (time) float64 -0.0028 -0.0097 -0.0153 ... 0.244 0.2186 0.1897
Attributes: (12/26)
    description:               SWOT fixed grid
    corresponding_grid:        
    title:                     Altimeter like data simulated by SWOT simulator
    keywords:                  check keywords
    Conventions:               CF-1.6
    summary:                   SWOT grid data produced
    ...                        ...
    geospatial_lon_units:      degrees_east
    project:                   SWOT
    date_created:              2018-10-12T12:39:50Z
    date_modified:             2018-10-12T12:39:50Z
    keywords_vocabulary:       NASA
    references:                Gaultier, L., C. Ubelmann, and L.-L. Fu, 2016:...

xarray.Dataset

• Dimensions:
  • time: 205232
• Coordinates: (3)
  • time
    (time)
    datetime64[ns]
    2012-10-01T00:45:50.676369 ... 2...

    array(['2012-10-01T00:45:50.676369000', '2012-10-01T00:45:51.665139000',
           '2012-10-01T00:45:52.489113000', ..., '2013-09-30T14:29:45.527443000',
           '2013-09-30T14:29:45.527443000', '2013-09-30T14:29:48.164162000'],
          dtype='datetime64[ns]')

  • lon
    (time)
    float64
    -57.81 -57.83 ... -58.31 -58.33

    units :

    deg

    standard_name :

    longitude

    long_name :

    Longitude

    valid_min :

    0

    valid_max :

    359999999

    array([-57.807968, -57.828224, -57.848416, ..., -58.289344, -58.309536,
           -58.329792])

  • lat
    (time)
    float64
    42.96 42.9 42.84 ... 42.94 43.0

    units :

    deg

    standard_name :

    latitude

    long_name :

    Latitude

    valid_min :

    -80000000

    valid_max :

    80000000

    array([42.956744, 42.899372, 42.841996, ..., 42.88094 , 42.938324,
           42.995704])

• Data variables: (1)
  • ssh
    (time)
    float64
    -0.0028 -0.0097 ... 0.2186 0.1897

    valid_min :

    -1500000

    valid_max :

    1500000

    units :

    m

    long_name :

    Sea Surface Height

    standard_name :

    sea_surface_height

    array([-0.0028, -0.0097, -0.0153, ...,  0.244 ,  0.2186,  0.1897])

• Indexes: (1)
  • time
    PandasIndex

    PandasIndex(DatetimeIndex(['2012-10-01 00:45:50.676369', '2012-10-01 00:45:51.665139',
                   '2012-10-01 00:45:52.489113', '2012-10-01 00:45:53.477883',
                   '2012-10-01 00:45:54.466652', '2012-10-01 00:45:55.455422',
                   '2012-10-01 00:45:56.444191', '2012-10-01 00:45:57.432961',
                   '2012-10-01 00:45:58.421730', '2012-10-01 00:45:59.410500',
                   ...
                   '2013-09-30 14:29:37.617287', '2013-09-30 14:29:40.254006',
                   '2013-09-30 14:29:40.254006', '2013-09-30 14:29:42.890725',
                   '2013-09-30 14:29:42.890725', '2013-09-30 14:29:42.890725',
                   '2013-09-30 14:29:45.527443', '2013-09-30 14:29:45.527443',
                   '2013-09-30 14:29:45.527443', '2013-09-30 14:29:48.164162'],
                  dtype='datetime64[ns]', name='time', length=205232, freq=None))

• Attributes: (26)

  description :

  SWOT fixed grid

  corresponding_grid :

  title :

  Altimeter like data simulated by SWOT simulator

  keywords :

  check keywords

  Conventions :

  CF-1.6

  summary :

  SWOT grid data produced

  Metadata_Conventions :

  Unidata Dataset Discovery v1.0

  history :

  Grid File created by swotsimulator version 2.21

  processing_level :

  L2

  standard_name_vocabulary :

  CF-1.6

  creator_name :

  Lucile Gaultier and Clement Ubelmann

  creator_email :

  lucile.gaultier@gmail.com

  publisher_url :

  github/SWOTSimulator/

  time_coverage_start :

  0.3383046

  time_coverage_end :

  9.895123

  geospatial_lat_min :

  26.46

  geospatial_lat_max :

  66.15

  geospatial_lat_units :

  degrees_north

  geospatial_lon_max :

  360.00

  geospatial_lon_min :

  0.00

  geospatial_lon_units :

  degrees_east

  project :

  SWOT

  date_created :

  2018-10-12T12:39:50Z

  date_modified :

  2018-10-12T12:39:50Z

  keywords_vocabulary :

  NASA

  references :

  Gaultier, L., C. Ubelmann, and L.-L. Fu, 2016: The Challenge of Using Future SWOT Data for Oceanic Field Reconstruction. J. Atmos. Oceanic Technol., 33, 119-126, doi:10.1175/jtech-d-15-0160.1. http://dx.doi.org/10.1175/JTECH-D-15-0160.1.

%matplotlib inline

fig, ax = plt.subplots()

sub_ds = ds_nadir.sel(time=slice("2012-10-22","2012-10-25"))
variable = "ssh"
pts = ax.scatter(sub_ds.lon, sub_ds.lat, c=sub_ds[variable], s=0.1)
ax.set(
    xlabel="Longitude",
    ylabel="Latitude",
    xlim=[ds_nadir.lon.min(), ds_nadir.lon.max()],
    ylim=[ds_nadir.lat.min(), ds_nadir.lat.max()],
)

plt.colorbar(pts, label="Sea Surface Height [m]")
plt.tight_layout()
plt.show()

Data Challenge - 2021a - OSE

This data challenge looks at some real altimetry data for sea surface height.

!ls /gpfswork/rech/yrf/commun/data_challenges/dc21a_ose/test/train

dt_gulfstream_alg_phy_l3_20161201-20180131_285-315_23-53.nc
dt_gulfstream_h2g_phy_l3_20161201-20180131_285-315_23-53.nc
dt_gulfstream_j2g_phy_l3_20161201-20180131_285-315_23-53.nc
dt_gulfstream_j2n_phy_l3_20161201-20180131_285-315_23-53.nc
dt_gulfstream_j3_phy_l3_20161201-20180131_285-315_23-53.nc
dt_gulfstream_s3a_phy_l3_20161201-20180131_285-315_23-53.nc

files_nadir_dc21a = [
    "/gpfswork/rech/yrf/commun/data_challenges/dc21a_ose/test/train/dt_gulfstream_alg_phy_l3_20161201-20180131_285-315_23-53.nc",
    "/gpfswork/rech/yrf/commun/data_challenges/dc21a_ose/test/train/dt_gulfstream_h2g_phy_l3_20161201-20180131_285-315_23-53.nc",
    "/gpfswork/rech/yrf/commun/data_challenges/dc21a_ose/test/train/dt_gulfstream_j2g_phy_l3_20161201-20180131_285-315_23-53.nc",
    "/gpfswork/rech/yrf/commun/data_challenges/dc21a_ose/test/train/dt_gulfstream_j2n_phy_l3_20161201-20180131_285-315_23-53.nc",
    "/gpfswork/rech/yrf/commun/data_challenges/dc21a_ose/test/train/dt_gulfstream_j3_phy_l3_20161201-20180131_285-315_23-53.nc",
    "/gpfswork/rech/yrf/commun/data_challenges/dc21a_ose/test/train/dt_gulfstream_s3a_phy_l3_20161201-20180131_285-315_23-53.nc"
]

def preprocess_nadir_dc21a(da):
    
    # reorganized
    da = da.sortby("time").compute()
    
    # validate coordinates
    da = da.rename({"longitude": "lon", "latitude": "lat"})
    da = validate_latlon(da)
    da = validate_time(da)
    
    # # slice region
    # da = where_slice(da, "lon", -64.975, -55.007)
    # da = where_slice(da, "lat", 33.025, 42.9917)
    
    # # slice time period
    # da = da.sel(time=slice("2017-01-01", "2017-12-31"))
    
    # calculate SSH directly
    da = alongtrack_ssh(da)
    
    return da

ds_nadir = xr.open_mfdataset(
    files_nadir_dc21a, 
    preprocess=preprocess_nadir_dc21a,
    combine="nested",
    engine="netcdf4",
    concat_dim="time"
)

ds_nadir = ds_nadir.sortby("time").compute()

ds_nadir

<xarray.Dataset>
Dimensions:         (time: 1793663)
Coordinates:
  * time            (time) datetime64[ns] 2016-12-01T01:07:08.844252928 ... 2...
    lon             (time) float64 -50.5 -50.51 -50.52 ... -74.96 -74.97 -74.99
    lat             (time) float64 23.03 23.08 23.14 23.2 ... 25.21 25.27 25.33
Data variables:
    cycle           (time) float64 11.0 11.0 11.0 11.0 ... 193.0 193.0 193.0
    track           (time) float64 575.0 575.0 575.0 575.0 ... 83.0 83.0 83.0
    dac             (time) float32 -0.0565 -0.0569 -0.0571 ... -0.0842 -0.0847
    lwe             (time) float32 0.003 0.003 0.003 ... -0.005 -0.005 -0.005
    mdt             (time) float32 0.52 0.52 0.521 0.522 ... 0.621 0.622 0.623
    ocean_tide      (time) float64 0.0269 0.0288 0.0307 ... 0.302 0.3033 0.3046
    sla_filtered    (time) float32 0.101 0.103 0.105 0.109 ... 0.039 0.048 0.059
    sla_unfiltered  (time) float32 0.1 0.098 0.104 0.115 ... -0.01 0.069 0.046
    ssh             (time) float32 0.617 0.615 0.622 0.634 ... 0.616 0.696 0.674
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:          P23H16M26.458585S
    time_coverage_end:               2016-01-01T23:06:57Z
    time_coverage_resolution:        P1S
    time_coverage_start:             2015-12-31T23:50:30Z
    title:                           DT Altika Drifting Phase Global Ocean Al...

xarray.Dataset

• Dimensions:
  • time: 1793663
• Coordinates: (3)
  • time
    (time)
    datetime64[ns]
    2016-12-01T01:07:08.844252928 .....

    array(['2016-12-01T01:07:08.844252928', '2016-12-01T01:07:09.844252928',
           '2016-12-01T01:07:10.844252672', ..., '2018-01-30T22:55:33.260000000',
           '2018-01-30T22:55:34.284000000', '2018-01-30T22:55:35.308000000'],
          dtype='datetime64[ns]')

  • lon
    (time)
    float64
    -50.5 -50.51 ... -74.97 -74.99

    units :

    degrees_east

    standard_name :

    longitude

    long_name :

    Longitude of measurement

    array([-50.49506 , -50.509751, -50.524451, ..., -74.958654, -74.974188,
           -74.989733])

  • lat
    (time)
    float64
    23.03 23.08 23.14 ... 25.27 25.33

    units :

    degrees_north

    standard_name :

    latitude

    long_name :

    Latitude of measurement

    array([23.026026, 23.084948, 23.143869, ..., 25.212322, 25.26908 ,
           25.325836])

• Data variables: (9)
  • cycle
    (time)
    float64
    11.0 11.0 11.0 ... 193.0 193.0

    long_name :

    Cycle the measurement belongs to

    units :

    1

    array([ 11.,  11.,  11., ..., 193., 193., 193.])

  • track
    (time)
    float64
    575.0 575.0 575.0 ... 83.0 83.0

    long_name :

    Track in cycle the measurement belongs to

    units :

    1

    array([575., 575., 575., ...,  83.,  83.,  83.])

  • dac
    (time)
    float32
    -0.0565 -0.0569 ... -0.0842 -0.0847

    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([-0.0565    , -0.0569    , -0.0571    , ..., -0.0835    ,
           -0.08419999, -0.0847    ], dtype=float32)

  • lwe
    (time)
    float32
    0.003 0.003 0.003 ... -0.005 -0.005

    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([ 0.003,  0.003,  0.003, ..., -0.005, -0.005, -0.005], dtype=float32)

  • mdt
    (time)
    float32
    0.52 0.52 0.521 ... 0.622 0.623

    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([0.52000004, 0.52000004, 0.521     , ..., 0.62100005, 0.62200004,
           0.623     ], dtype=float32)

  • ocean_tide
    (time)
    float64
    0.0269 0.0288 ... 0.3033 0.3046

    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([0.0269, 0.0288, 0.0307, ..., 0.302 , 0.3033, 0.3046])

  • sla_filtered
    (time)
    float32
    0.101 0.103 0.105 ... 0.048 0.059

    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]-[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([0.101, 0.103, 0.105, ..., 0.039, 0.048, 0.059], dtype=float32)

  • sla_unfiltered
    (time)
    float32
    0.1 0.098 0.104 ... 0.069 0.046

    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]-[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([ 0.1       ,  0.098     ,  0.104     , ..., -0.01      ,
            0.06900001,  0.046     ], dtype=float32)

  • ssh
    (time)
    float32
    0.617 0.615 0.622 ... 0.696 0.674

    array([0.61700004, 0.615     , 0.622     , ..., 0.61600006, 0.69600004,
           0.674     ], dtype=float32)

• Indexes: (1)
  • time
    PandasIndex

    PandasIndex(DatetimeIndex(['2016-12-01 01:07:08.844252928',
                   '2016-12-01 01:07:09.844252928',
                   '2016-12-01 01:07:10.844252672',
                   '2016-12-01 01:07:11.844252928',
                   '2016-12-01 01:07:12.844252928',
                   '2016-12-01 01:07:13.844252928',
                   '2016-12-01 01:07:14.844252672',
                   '2016-12-01 01:07:15.844252928',
                   '2016-12-01 01:07:16.844252928',
                   '2016-12-01 01:07:17.844252928',
                   ...
                      '2018-01-30 22:55:26.092000',
                   '2018-01-30 22:55:27.116000256',
                      '2018-01-30 22:55:28.140000',
                      '2018-01-30 22:55:29.164000',
                      '2018-01-30 22:55:30.188000',
                   '2018-01-30 22:55:31.212000256',
                   '2018-01-30 22:55:32.236000256',
                      '2018-01-30 22:55:33.260000',
                      '2018-01-30 22:55:34.284000',
                      '2018-01-30 22:55:35.308000'],
                  dtype='datetime64[ns]', name='time', length=1793663, 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 :

  2019-03-28T10:13:03Z

  date_issued :

  2019-03-28T10:13:03Z

  date_modified :

  2019-03-28T10:13:03Z

  geospatial_lat_max :

  81.506562

  geospatial_lat_min :

  -78.185707

  geospatial_lat_resolution :

  0.03436100000000408

  geospatial_lat_units :

  degrees_north

  geospatial_lon_max :

  359.99886699999996

  geospatial_lon_min :

  0.013321

  geospatial_lon_resolution :

  0.017669000000012147

  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 :

  2019-03-28T10:13:03Z: 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 :

  2019

  project :

  COPERNICUS MARINE ENVIRONMENT MONITORING SERVICE (CMEMS)

  references :

  http://marine.copernicus.eu

  software_version :

  6.2_DUACS_DT2018_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 :

  P23H16M26.458585S

  time_coverage_end :

  2016-01-01T23:06:57Z

  time_coverage_resolution :

  P1S

  time_coverage_start :

  2015-12-31T23:50:30Z

  title :

  DT Altika Drifting Phase Global Ocean Along track SSALTO/DUACS Sea Surface Height L3 product

%matplotlib inline

fig, ax = plt.subplots()

sub_ds = ds_nadir.sel(time=slice("2017-01-10","2017-01-15"))
variable = "ssh"
pts = ax.scatter(sub_ds.lon, sub_ds.lat, c=sub_ds[variable], s=0.1)
ax.set(
    xlabel="Longitude",
    ylabel="Latitude",
    xlim=[ds_nadir.lon.min(), ds_nadir.lon.max()],
    ylim=[ds_nadir.lat.min(), ds_nadir.lat.max()],
)

plt.colorbar(pts, label="Sea Surface Height [m]")
plt.tight_layout()
plt.show()

Data Challenge - 2023a - OSE

This data challenge looks at the Western Mediterranean.

Data Challenge - 2023b - OSE

This data challenge looks at the Global dataset.

Pangeo-Forge

See this tutorial

TODO: Need to get the intake package working.

# from intake import open_catalog

# cat = open_catalog("https://raw.githubusercontent.com/pangeo-data/pangeo-datastore/master/intake-catalogs/ocean/altimetry.yaml")
# list(cat)