Time series

[malerany]

Integrate functions to generate and plot time series into the methods of the DataDiagnostics class defined in src/pyhanami/diags/Diagnostics.py and the function time_series_plot defined in src/pyhanami/utils/plot.py, respectively.

Assignees: @malerany

Migrated from GitLab: https://gitlab.earth.bsc.es/ces/hanami/pyhanami/-/issues/10

[malerany]

Hi @kellekai, I believe you had already implemented the functions related to the generation of time series plots. Could you please point me to the scripts where they are defined so I can integrate them into the package?

[malerany]

In GitLab by @kellekai on May 15, 2025, 15:26 UTC:

Hi @malerany, yes, I think I like this plot most so far. I am just posting the full code here, it should be possible to run it as an example in a notebook. And then modify it to our needs.

from scipy.stats import bootstrap
import xarray as xr
import numpy as np
rng = np.random.default_rng()
import matplotlib.pyplot as plt

def compute_area_weighted_dataset(dr):
    return dr.weighted(np.cos(np.deg2rad(dr.lat)))

def compute_timeseries_annual(dr):
    dr_weighted = compute_area_weighted_dataset(dr)
    dr_annual = dr_weighted.mean(['lat','lon']).coarsen(time=12).mean()
    return dr_annual

atm2d_units_dct = {
    'tas'           : 'K',
    'psl'           : 'Pa',
    'pr'            : 'kg m**-2 s**-1',
    'tauv'          : 'N m**-2',
    'tauu'          : 'N m**-2',
    'net_sfc'       : 'W m**-2',
}

atm3d_units_dct = {
    'ta200'         : 'K', 
    'ta850'         : 'K', 
    'ua200'         : 'm s**-1', 
    'ua850'         : 'm s**-1', 
    'va200'         : 'm s**-1', 
    'va850'         : 'm s**-1', 
    'hus300'        : 'kg kg**-1',
    'hus850'        : 'kg kg**-1'
 }

units_dct = {
    'atm2d' : atm2d_units_dct,
    'atm3d' : atm3d_units_dct
}

realm = {
    'tas'           : 'atm2d',
    'psl'           : 'atm2d',
    'pr'            : 'atm2d',
    'tauv'          : 'atm2d',
    'tauu'          : 'atm2d',
    'net_sfc'       : 'atm2d',
    'ta200'         : 'atm3d',
    'ta850'         : 'atm3d',
    'ua200'         : 'atm3d',
    'ua850'         : 'atm3d',
    'va200'         : 'atm3d',
    'va850'         : 'atm3d',
    'hus300'        : 'atm3d',
    'hus850'        : 'atm3d',
}

vars = ['tas', 'psl', 'pr', 'tauv', 'tauu', 'net_sfc', 'ta200', 'ta850', 'ua200', 'ua850', 'va200', 'va850', 'hus300', 'hus850']

# a26d -> double precision (50 members)
members_a26d = list(range(25))
members_a26d.remove(1)
members_a26d.remove(17)

# a26j -> single precision (25 members - 2 failed [fc1, fc17])
members_a26j = list(range(25))
members_a26j.remove(1)
members_a26j.remove(17)

fig, ax = plt.subplots(4,4, figsize=(32,20))

for id, var in enumerate(vars):
    
    #ens_full = xr.open_dataset(f'/gpfs/scratch/bsc32/bsc032655/lab/DestinE/single_double_ensemble_25km_replicability_test/data/control_data/a26d_{realm[var]}_M50_{var}.nc')
    #ens = ens_full.isel(realization=range(25,50,1))
    ens_cont = xr.open_dataset(f'/gpfs/scratch/bsc32/bsc032655/lab/DestinE/single_double_ensemble_25km_replicability_test/data/30_years_atm/a26d_cont_{realm[var]}_M23_Y30_{var}.nc') 
    ens_very = xr.open_dataset(f'/gpfs/scratch/bsc32/bsc032655/lab/DestinE/single_double_ensemble_25km_replicability_test/data/30_years_atm/a26d_{realm[var]}_M23_Y30_{var}.nc') 
    ens_test = xr.open_dataset(f'/gpfs/scratch/bsc32/bsc032655/lab/DestinE/single_double_ensemble_25km_replicability_test/data/30_years_atm/a26j_{realm[var]}_M23_Y30_{var}.nc') 
    
    ens_cont_ts = compute_timeseries_annual(ens_cont)
    ens_cont_ts_mean = ens_cont_ts.mean('realization')
    
    ens_very_ts = compute_timeseries_annual(ens_very)
    ens_very_ts_mean = ens_very_ts.mean('realization')
    
    ens_test_ts = compute_timeseries_annual(ens_test)
    ens_test_ts_mean = ens_test_ts.mean('realization')
    
    print(var)
    
    row = id//4
    col = id%4
    
    q025 = []
    q975 = []
    for i in range(30):
        # compute 95% confidence interval for the ensemble mean (dark blue in plots)
        bs_res = bootstrap((ens_cont_ts[var].isel(time=i).values,), np.mean, confidence_level=0.95)
        q025.append(bs_res.confidence_interval.low)
        q975.append(bs_res.confidence_interval.high)
    
    # compute 5 and 95 percent quantiles for the ensemble spread (light blue in plots)
    quantiles_data = np.quantile(ens_cont_ts[var].values, [0.05, 0.95], axis=0)
    
    ax[row,col].plot(range(1990,2020,1), ens_very_ts_mean[var].values, lw=2, ls='-.', label=f'{var} - DP', color='tab:green')
    ax[row,col].plot(range(1990,2020,1), ens_test_ts_mean[var].values, lw=2, ls='-.', label=f'{var} - SP', color='tab:red')
    
    ax[row,col].plot(range(1990,2020,1), ens_cont_ts_mean[var].values, lw=2, label=f'{var} - Control', color='tab:blue')
    for i in range(23):
        ax[row,col].plot(range(1990,2020,1), ens_cont_ts[var].isel(realization=i).values, color='tab:blue', alpha=0.2)

    ax[row,col].fill_between(range(1990,2020,1), 
         q025, 
         q975,
         facecolor='C0', alpha=0.6)

    ax[row,col].fill_between(range(1990,2020,1), 
         quantiles_data[0], 
         quantiles_data[1],
         facecolor='C0', alpha=0.2)
     
    ax[row,col].set_title(f'Timeseries of annual mean {var}', fontsize=18)
    ax[row,col].legend(fontsize=16)
    ax[row,col].set_xlabel("Year", fontsize=14)
    ax[row,col].set_ylabel(f"{var} [{units_dct[realm[var]][var]}]", fontsize=16)
    ax[row,col].set_xticks(range(1990,2020,1))
    ax[row,col].tick_params(axis='y', labelsize=12)
    ax[row,col].set_xticklabels(range(1990,2020,1), rotation=45, ha='right', fontsize=12)
    ax[row,col].grid()
    
    
ax[3,2].remove()
ax[3,3].remove()

plt.tight_layout()

[malerany]

Great, thank you!

[malerany]

I implemented the time series methods within the DataDiagnostics() class and the corresponding function in utils/plot.py.

They still need to be fully tested.

[malerany]

In GitLab by @kellekai on Jul 25, 2025, 09:57 UTC:

Hi @malerany I noticed that the time series is stored always using the same name. Perhaps we can name it after the reference name + variable? Something like /path/to/output/timeseries.<reference1_name>.<reference2_name>..png ?

[malerany]

In GitLab by @kellekai on Jul 25, 2025, 10:00 UTC:

@malerany, also, I think we should add observations to the plots. I would say for a start, we can provide ERA5 and JRA55 with monthly means. We could store the annual means directly, in that way, we can easily store them in the repository. WDYT?

At some point, some CMIP6 time series would be great too. in that case we can provide timeseries annual mean data from 1850 - 2100.

[malerany]

Hi @kellekai, I just changed the default names for the plots (time series, spatial plots and matrix plot) to <plot_type>_<var_name>_<sim_name_1>-<sim_name_2>.png.

[malerany]

Hi @kellekai, I think it is a good idea to make the observations available through the repository itself. However, the period for which JRA55 is available is already covered by ERA5, right? Perhaps we could provide ERA5 and NOAA instead? This way we would cover a larger period (from 1850 onwards). What do you think?

[malerany]

In GitLab by @kellekai on Jul 28, 2025, 08:29 UTC:

I think that it would be good to provide as many different sources for observations as possible. JRA55 and ERA5 are great because they are easy to handle, and they are different, as they are generated by a different model. What would be available from NOAA? Do they also provide reanalysis / continuous observations? MERRA is also an option. Also, EN4 for the ocean. There is also ORAS5 for the ocean. Maybe we can make a list to collect obs/recons and corresponding periods, variables, frequency, and resolutions? We can do it together inside an online doc, for instance. WDYT?

[malerany]

Hi @kellekai, sounds like a good idea. I created a spreadsheet with a tag called ‘Obs_datasets’ with the characteristics of each observations and a tag called ‘Available_variables’ with the variables available in each one (for now I only included the variables that we used for the replicability test).

Here is the link: Observations_pyhanami

I added the info for ERA5 and will continue filling it in.

[malerany]

Hi, I found this website reanalysis.org which provides references to many reanalysis datasets.

[malerany]

I listed 11 general reanalyses and 2 focused on the ocean in the spreadsheet Observations_pyhanami. Only one of them (MERRA-2) contains all the variables that we used for the replicability test.

[malerany]

In GitLab by @kellekai on Aug 4, 2025, 12:16 UTC:

Thanks Marta, I gonna take a look.

[malerany]

I implemented the option to include observations in the time series plots. For now, apart from setting the observations to True, a path and a name for the observations have to be provided. Something like this:

diags.time_series_plots('tas', ['lens2_cmip6', 'lens2_smbb'], obs=True, obs_paths='/gpfs/scratch/bsc32/bsc032155/HANAMI/pyhanami_trials', obs_names='era5')

where obs_paths and obs_names can also be lists to plot more than one observations dataset.

I was thinking that, once the observations are downloaded together with the package, we could set an observations dataset as the default to be added to the plot if the user wants to include observations for comparison without having to select a specific observations dataset.

[malerany]

In GitLab by @kellekai on Aug 9, 2025, 14:30 UTC:

Hi Marta, yes, indeed. I was thinking that we can download the observations and store them as global mean timeseries. I think it would be great to have the option of plotting the timeseries for different time averages such as daily, monthly, and annual. Also, I think it might be good to have the option of plotting anomalies instead of absolute values. Because, sometimes observations are provided as anomalies.

Per default, we should enable JRA55 and ERA5. If the time series is out of the period for those, we can print a warning. Then, of course, we have to provide a mean of selecting the observations specifically.