
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)
Parameters
Parameters
[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
Input files
Input files
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>
Spectral Analysis
Spectral Analysis
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.021172.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)
[17]:
plot_psd_scores(psd_output_filename)
[17]:
The interactive plot above allows you to explore the spectral metrics by latitude / longitude box
[ ]: