header

Alongtrack spectral metrics

2023-04-27 DUACS_SSH_BENCHMARK_DEMO


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

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

Sea Surface Height from Saral/AltiKa

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

Sea Surface Height maps to evaluate

[6]:
%%time
list_of_maps = sorted(glob('../data/maps/DUACS_global_allsat-alg/*.nc'))
ds_maps = xr.open_mfdataset(list_of_maps, combine='nested', concat_dim='time')
ds_maps = ds_maps.sel(time=slice(time_min, time_max))
ds_maps
CPU times: user 6.2 s, sys: 741 ms, total: 6.94 s
Wall time: 8.42 s
[6]:
<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>

2.1 Interpolate sea surface height maps onto along-track positions

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

2.2 Compute Spectral scores

[15]:
compute_psd_scores_v2(ds_interp, psd_output_filename, lenght_scale=segment_lenght,method_name=method_name)
2023-07-05 17:03:17 INFO     Segment computation...
2023-07-05 17:03:33 INFO     Spectral analysis...
2023-07-05 17:19:27 INFO     Saving ouput...
2023-07-05 17:19:35 INFO     PSD file saved as: ../results/psd_sla_duacs_geos.nc
[16]:
# Plot effective resolution
# Hvplot
# plot_effective_resolution(psd_output_filename)
# Matplotlib
plot_effective_resolution_png(psd_output_filename)
../_images/gallery_alongtrack_spectral_description_17_0.png
[17]:
plot_psd_scores(psd_output_filename)
[17]:

The interactive plot above allows you to explore the spectral metrics by latitude / longitude box

[ ]: