CALCUL AND SAVE SST STATS

import numpy as N
import xarray as xr
import pandas as pd
import matplotlib.pyplot as plt
import sys
import dask
from os import path
from dask.distributed import Client, LocalCluster, progress
import time
from tools import my_percentile
from scipy.optimize import least_squares,curve_fit
import dask_hpcconfig
from dask_jobqueue import PBSCluster
import glob

conseils dask

pour Datarmor il vaut mieux avoir 1 thread/worker il vaut mieux avoir des chunks assez fa

#pour l'ensemble des pas de temps
#cluster = PBSCluster(processes=3,cores=9) #16min 44
cluster = PBSCluster(processes=3,cores=3) #16min 44

cluster.scale(jobs=4)
print(cluster.dashboard_link)

/home1/datahome/mcaillau/conda-env/pydask3/lib/python3.10/site-packages/dask_jobqueue/core.py:255: FutureWarning: job_extra has been renamed to job_extra_directives. You are still using it (even if only set to []; please also check config files). If you did not set job_extra_directives yet, job_extra will be respected for now, but it will be removed in a future release. If you already set job_extra_directives, job_extra is ignored and you can remove it.
  warnings.warn(warn, FutureWarning)
/home1/datahome/mcaillau/conda-env/pydask3/lib/python3.10/site-packages/dask_jobqueue/pbs.py:82: FutureWarning: project has been renamed to account as this kwarg was used wit -A option. You are still using it (please also check config files). If you did not set account yet, project will be respected for now, but it will be removed in a future release. If you already set account, project is ignored and you can remove it.
  warnings.warn(warn, FutureWarning)

http://10.148.1.73:8787/status

/home1/datahome/mcaillau/conda-env/pydask3/lib/python3.10/site-packages/dask_jobqueue/core.py:255: FutureWarning: job_extra has been renamed to job_extra_directives. You are still using it (even if only set to []; please also check config files). If you did not set job_extra_directives yet, job_extra will be respected for now, but it will be removed in a future release. If you already set job_extra_directives, job_extra is ignored and you can remove it.
  warnings.warn(warn, FutureWarning)
/home1/datahome/mcaillau/conda-env/pydask3/lib/python3.10/site-packages/dask_jobqueue/pbs.py:82: FutureWarning: project has been renamed to account as this kwarg was used wit -A option. You are still using it (please also check config files). If you did not set account yet, project will be respected for now, but it will be removed in a future release. If you already set account, project is ignored and you can remove it.
  warnings.warn(warn, FutureWarning)

# explicitly connect to the cluster we just created
client = Client(cluster)

#read grid
ds_grid=xr.open_dataset('/home/shom_simuref/CROCO/ODC/CONFIGS/MEDITERRANEE_GLOBALE/CROCO_FILES/test2.nc')

path="/home/shom_simuref/CROCO/ODC/SIMU-RESULT/HINDCAST_2012_2013/OUTPUTS_201207_201307/"

list_model=glob.glob(path+"croco_his_surf_2*.nc")

chunks={'time':121,'xi_rho':315,'eta_rho':878}
chunks={'time':800,'s_rho':-1,'xi_rho':865,'eta_rho':936}
date_start="2012-06-01"
date_end="2013-08-01"
time_range=slice(date_start,date_end)

#store data dir
stat_dir='/home/shom_simuref/CROCO/ODC/POSTPROC/SST/'

ds=xr.open_mfdataset(list_model, parallel=False,chunks="auto",
                        concat_dim="time", combine="nested",
                        data_vars='minimal', coords='minimal', compat='override')
#                        preprocess=partial_func,data_vars='minimal', coords='minimal', compat='override')

print('data size in GB {:0.2f}\n'.format(ds.nbytes / 1e9))

data size in GB 425.54

#ds=ds.isel(time=slice(0,-1,12))

#because of MPINOLAND
for var in ["xi_rho","eta_rho","xi_u","eta_v","lon_rho","lat_rho","lon_u","lat_v"]:
    ds[var]=ds_grid[var]

#remove duplicates
_,index=N.unique(ds.time,return_index=True)
ds=ds.isel(time=index)

#select time range
ds=ds.sel(time=time_range)

ds=ds.chunk(chunks)
ds.temp

<xarray.DataArray 'temp' (time: 8760, eta_rho: 936, xi_rho: 2595)>
dask.array<rechunk-merge, shape=(8760, 936, 2595), dtype=float32, chunksize=(800, 936, 865), chunktype=numpy.ndarray>
Coordinates:
  * xi_rho   (xi_rho) int64 0 1 2 3 4 5 6 ... 2588 2589 2590 2591 2592 2593 2594
  * eta_rho  (eta_rho) int64 0 1 2 3 4 5 6 7 ... 928 929 930 931 932 933 934 935
  * time     (time) datetime64[ns] 2012-07-01T13:00:00 ... 2013-07-01T12:00:00
    lon_rho  (eta_rho, xi_rho) float64 dask.array<chunksize=(936, 865), meta=np.ndarray>
    lat_rho  (eta_rho, xi_rho) float64 dask.array<chunksize=(936, 865), meta=np.ndarray>
Attributes:
    long_name:  SST
    units:      Celsius
    field:      

xarray.DataArray

'temp'

• time: 8760
• eta_rho: 936
• xi_rho: 2595

• dask.array<chunksize=(800, 936, 865), meta=np.ndarray>

  Array Chunk Bytes 79.26 GiB 2.41 GiB Shape (8760, 936, 2595) (800, 936, 865) Dask graph 33 chunks in 81 graph layers Data type float32 numpy.ndarray

• Coordinates: (5)
  • xi_rho
    (xi_rho)
    int64
    0 1 2 3 4 ... 2591 2592 2593 2594

    array([   0,    1,    2, ..., 2592, 2593, 2594])

  • eta_rho
    (eta_rho)
    int64
    0 1 2 3 4 5 ... 931 932 933 934 935

    array([  0,   1,   2, ..., 933, 934, 935])

  • time
    (time)
    datetime64[ns]
    2012-07-01T13:00:00 ... 2013-07-...

    long_name :

    time since initialization

    field :

    time, scalar, series

    standard_name :

    time

    axis :

    T

    array(['2012-07-01T13:00:00.000000000', '2012-07-01T14:00:00.000000000',
           '2012-07-01T15:00:00.000000000', ..., '2013-07-01T10:00:00.000000000',
           '2013-07-01T11:00:00.000000000', '2013-07-01T12:00:00.000000000'],
          dtype='datetime64[ns]')

  • lon_rho
    (eta_rho, xi_rho)
    float64
    dask.array<chunksize=(936, 865), meta=np.ndarray>

    Array Chunk Bytes 18.53 MiB 6.18 MiB Shape (936, 2595) (936, 865) Dask graph 3 chunks in 2 graph layers Data type float64 numpy.ndarray

  • lat_rho
    (eta_rho, xi_rho)
    float64
    dask.array<chunksize=(936, 865), meta=np.ndarray>

    Array Chunk Bytes 18.53 MiB 6.18 MiB Shape (936, 2595) (936, 865) Dask graph 3 chunks in 2 graph layers Data type float64 numpy.ndarray

• Indexes: (3)
  • xi_rho
    PandasIndex

    PandasIndex(RangeIndex(start=0, stop=2595, step=1, name='xi_rho'))

  • eta_rho
    PandasIndex

    PandasIndex(RangeIndex(start=0, stop=936, step=1, name='eta_rho'))

  • time
    PandasIndex

    PandasIndex(DatetimeIndex(['2012-07-01 13:00:00', '2012-07-01 14:00:00',
                   '2012-07-01 15:00:00', '2012-07-01 16:00:00',
                   '2012-07-01 17:00:00', '2012-07-01 18:00:00',
                   '2012-07-01 19:00:00', '2012-07-01 20:00:00',
                   '2012-07-01 21:00:00', '2012-07-01 22:00:00',
                   ...
                   '2013-07-01 03:00:00', '2013-07-01 04:00:00',
                   '2013-07-01 05:00:00', '2013-07-01 06:00:00',
                   '2013-07-01 07:00:00', '2013-07-01 08:00:00',
                   '2013-07-01 09:00:00', '2013-07-01 10:00:00',
                   '2013-07-01 11:00:00', '2013-07-01 12:00:00'],
                  dtype='datetime64[ns]', name='time', length=8760, freq=None))

• Attributes: (3)

  long_name :

  SST

  units :

  Celsius

  field :

coord_dict={"xi_rho":"X","eta_rho":"Y","xi_u":"X_U","eta_v":"Y_V"}
ds=ds.assign_coords({"X":ds_grid.lon_rho[0,:], "Y":ds_grid.lat_rho[:,0]})
ds2=ds.swap_dims(coord_dict)

sst=ds2.temp

sst

<xarray.DataArray 'temp' (time: 8760, Y: 936, X: 2595)>
dask.array<rechunk-merge, shape=(8760, 936, 2595), dtype=float32, chunksize=(800, 936, 865), chunktype=numpy.ndarray>
Coordinates:
    xi_rho   (X) int64 0 1 2 3 4 5 6 7 ... 2588 2589 2590 2591 2592 2593 2594
    eta_rho  (Y) int64 0 1 2 3 4 5 6 7 8 ... 927 928 929 930 931 932 933 934 935
  * time     (time) datetime64[ns] 2012-07-01T13:00:00 ... 2013-07-01T12:00:00
    lon_rho  (Y, X) float64 dask.array<chunksize=(936, 865), meta=np.ndarray>
    lat_rho  (Y, X) float64 dask.array<chunksize=(936, 865), meta=np.ndarray>
  * X        (X) float64 -7.0 -6.983 -6.967 -6.95 ... 36.18 36.2 36.22 36.23
  * Y        (Y) float64 30.23 30.25 30.27 30.28 30.3 ... 45.77 45.78 45.8 45.82
Attributes:
    long_name:  SST
    units:      Celsius
    field:      

xarray.DataArray

'temp'

• time: 8760
• Y: 936
• X: 2595

• dask.array<chunksize=(800, 936, 865), meta=np.ndarray>

  Array Chunk Bytes 79.26 GiB 2.41 GiB Shape (8760, 936, 2595) (800, 936, 865) Dask graph 33 chunks in 81 graph layers Data type float32 numpy.ndarray

• Coordinates: (7)
  • xi_rho
    (X)
    int64
    0 1 2 3 4 ... 2591 2592 2593 2594

    array([   0,    1,    2, ..., 2592, 2593, 2594])

  • eta_rho
    (Y)
    int64
    0 1 2 3 4 5 ... 931 932 933 934 935

    array([  0,   1,   2, ..., 933, 934, 935])

  • time
    (time)
    datetime64[ns]
    2012-07-01T13:00:00 ... 2013-07-...

    long_name :

    time since initialization

    field :

    time, scalar, series

    standard_name :

    time

    axis :

    T

    array(['2012-07-01T13:00:00.000000000', '2012-07-01T14:00:00.000000000',
           '2012-07-01T15:00:00.000000000', ..., '2013-07-01T10:00:00.000000000',
           '2013-07-01T11:00:00.000000000', '2013-07-01T12:00:00.000000000'],
          dtype='datetime64[ns]')

  • lon_rho
    (Y, X)
    float64
    dask.array<chunksize=(936, 865), meta=np.ndarray>

    Array Chunk Bytes 18.53 MiB 6.18 MiB Shape (936, 2595) (936, 865) Dask graph 3 chunks in 2 graph layers Data type float64 numpy.ndarray

  • lat_rho
    (Y, X)
    float64
    dask.array<chunksize=(936, 865), meta=np.ndarray>

    Array Chunk Bytes 18.53 MiB 6.18 MiB Shape (936, 2595) (936, 865) Dask graph 3 chunks in 2 graph layers Data type float64 numpy.ndarray

  • X
    (X)
    float64
    -7.0 -6.983 -6.967 ... 36.22 36.23

    array([-7.      , -6.983333, -6.966667, ..., 36.199827, 36.216494, 36.23316 ])

  • Y
    (Y)
    float64
    30.23 30.25 30.27 ... 45.8 45.82

    array([30.233333, 30.25    , 30.266666, ..., 45.783271, 45.799938, 45.816604])

• Indexes: (3)
  • time
    PandasIndex

    PandasIndex(DatetimeIndex(['2012-07-01 13:00:00', '2012-07-01 14:00:00',
                   '2012-07-01 15:00:00', '2012-07-01 16:00:00',
                   '2012-07-01 17:00:00', '2012-07-01 18:00:00',
                   '2012-07-01 19:00:00', '2012-07-01 20:00:00',
                   '2012-07-01 21:00:00', '2012-07-01 22:00:00',
                   ...
                   '2013-07-01 03:00:00', '2013-07-01 04:00:00',
                   '2013-07-01 05:00:00', '2013-07-01 06:00:00',
                   '2013-07-01 07:00:00', '2013-07-01 08:00:00',
                   '2013-07-01 09:00:00', '2013-07-01 10:00:00',
                   '2013-07-01 11:00:00', '2013-07-01 12:00:00'],
                  dtype='datetime64[ns]', name='time', length=8760, freq=None))

  • X
    PandasIndex

    PandasIndex(Float64Index([               -7.0,  -6.983333399999999,          -6.9666668,
                           -6.9500002,          -6.9333336, -6.9166669999999995,
                           -6.9000004,          -6.8833338,          -6.8666672,
                  -6.8500005999999996,
                  ...
                            36.083161,          36.0998276,          36.1164942,
                           36.1331608,          36.1498274,           36.166494,
                           36.1831606,          36.1998272,          36.2164938,
                           36.2331604],
                 dtype='float64', name='X', length=2595))

  • Y
    PandasIndex

    PandasIndex(Float64Index([        30.2333333, 30.249999900000002,         30.2666665,
                          30.2833331,         30.2999997,         30.3166663,
                  30.333332900000002, 30.349999500000003,         30.3666661,
                          30.3833327,
                  ...
                  45.666604899999996,         45.6832715,         45.6999381,
                          45.7166047,         45.7332713,         45.7499379,
                          45.7666045,         45.7832711,         45.7999377,
                  45.816604299999995],
                 dtype='float64', name='Y', length=936))

• Attributes: (3)

  long_name :

  SST

  units :

  Celsius

  field :

MONTHLY MEAN

%%time
#4min12s on 4 workers and 3 proc par worker
sst=ds2.temp
print('compute mean')
sst_mean = sst.resample(time="1MS").mean().compute()
sst_mean=sst_mean.assign_attrs({"start_date":date_start,"end_date":date_end})
#write results for faster treatment
sst_mean.to_netcdf(f'{stat_dir}/mean_sst_croco.nc',mode='w')

compute mean
CPU times: user 18.9 s, sys: 2.13 s, total: 21.1 s
Wall time: 4min 12s

DAILY MEAN

%%time
#4min12s on 4 workers and 3 proc par worker
sst=ds2.temp
print('compute daily mean')
sst_daily_mean = sst.resample(time="1D").mean().compute()
sst_daily_mean=sst_daily_mean.assign_attrs({"start_date":date_start,"end_date":date_end})
#write results for faster treatment
sst_daily_mean.to_netcdf(f'{stat_dir}/mean_daily_sst_croco.nc',mode='w')

compute daily mean
CPU times: user 39.2 s, sys: 9.66 s, total: 48.8 s
Wall time: 6min 31s

PERCENTILES

from tools import my_percentile

re chunks to percentile

#sst=sst.chunk({'time':-1,'Y':52,'X':173})

#%%time
#per95=my_percentile(sst,0.95,"time")
#per95=per95.compute()
#per95.to_netcdf(f'{stat_dir}/sst_percentile_95.nc',mode='w')

#%%time
#per5=my_percentile(sst,0.05,"time")
#per5=per5.compute()
#per5.to_netcdf(f'{stat_dir}/sst_percentile_5.nc',mode='w')

SEVIRI

seviri_path="/home/shom_simuref/CROCO/ODC/DATA/SST_SEVIRI/"

chunks_sevi={'time':444,"lon":-1,"lat":-1}

ds_sevi=xr.open_mfdataset(glob.glob(seviri_path+'201[2-3]*.nc'),chunks=chunks_sevi,
                         combine="by_coords",
                          data_vars='minimal', coords='minimal', compat='override')

reduce to croco time bnds

time_bnds=slice(ds.time[0],ds.time[-1])

ds_sevi=ds_sevi.sel(time=time_bnds)

sst_sevi=ds_sevi.sea_surface_temperature-273.15

monthly mean

sst_sevi_mean=sst_sevi.resample(time='1MS').mean()

%%time
sst_sevi_mean=sst_sevi_mean.compute()
sst_sevi=sst_sevi.assign_attrs({"start_date":date_start,"end_date":date_end})
sst_sevi_mean.to_netcdf(f'{stat_dir}/mean_sst_sevi.nc',mode="w")

CPU times: user 20 ms, sys: 0 ns, total: 20 ms
Wall time: 59.5 ms

daily mean

sst_sevi_daily_mean=sst_sevi.resample(time='1D').mean()
sst_sevi_daily_mean=sst_sevi_daily_mean.assign_attrs({"start_date":date_start,"end_date":date_end})
sst_sevi_daily_mean.to_netcdf(f'{stat_dir}/mean_daily_sst_sevi.nc',mode="w")

percentiles

re chunk for PERCENTILE

sst_sevi=sst_sevi.chunk({'time':-1,"lat":-1,"lon":-1})

sst_sevi2=sst_sevi.interpolate_na(method="linear",limit=None,max_gap=None,dim="time")

sst_sevi2=sst_sevi2.ffill('time')
sst_sevi2=sst_sevi2.bfill('time')

#%%time
##35 sec sur 4 jobs de 6 processes et 12 cores ( 2 threads/worker)
#sst_sevi_p95=my_percentile(sst_sevi2,0.95,"time").compute()
#sst_sevi_p95.to_netcdf(f'{stat_dir}/sst_sevi_p95.nc',mode="w")

#%%time
#sst_sevi_p5=my_percentile(sst_sevi2,0.05,"time").compute()
#sst_sevi_p5.to_netcdf(f'{stat_dir}/sst_sevi_p5.nc',mode="w")

RMSE CROCO-SEVIRI

First compute seviri

%%time
#1min 21s sur 4 workers avec 3 procs par worker
sst_sevi2=sst_sevi2.compute()

CPU times: user 5.1 s, sys: 3.95 s, total: 9.05 s
Wall time: 1min 21s

interp croco on seviri horizontal grid

sst_croco=sst.interp(X=sst_sevi.lon,Y=sst_sevi.lat)

CPU times: user 36 ms, sys: 0 ns, total: 36 ms
Wall time: 115 ms

sst_croco=sst_croco.chunk({'time':2000})

%%time
#6 min 23s sur 4 workers avec 3 procs par worker
sst_croco=sst_croco.compute()

CPU times: user 41.7 s, sys: 35.6 s, total: 1min 17s
Wall time: 6min 23s

interp croco on seviri time grid

%%time
#40 s sur 4 workers avec 3 procs par worker
sst_croco2=sst_croco.interp(time=sst_sevi2.time).compute()

CPU times: user 16.2 s, sys: 25.8 s, total: 42 s
Wall time: 40.1 s

Compute the RMSE

import xskillscore as xs

calcul of monthly rmse

ds_merge=xr.merge([sst_croco2.rename('croco'),sst_sevi2.rename('seviri')])

#custom function to use with resample in xarray
def calc_rmse(data):
    return xs.rmse(data.croco,data.seviri,dim="time",skipna=True)

%%time
#22s sur 4 workers avec 3 procs par worker
rmse_month=ds_merge.resample(time='1MS').apply(calc_rmse).compute()

CPU times: user 13.1 s, sys: 9.64 s, total: 22.7 s
Wall time: 22 s

rmse_month=rmse_month.assign_attrs({"start_date":date_start,"end_date":date_end})
rmse_month=rmse_month.rename('rmse')
rmse_month.to_netcdf(f'{stat_dir}/rmse_croco_sevi_monthly.nc',mode='w')

calcul of annual rmse

%%time
#31s sur 4 workers avec 3 procs par worker
rmse=xs.rmse(sst_croco2,sst_sevi2,dim="time",skipna=True)
rmse=rmse.compute()
rmse=rmse.rename('rmse')
rmse=rmse.assign_attrs({"start_date":date_start,"end_date":date_end})
rmse.to_netcdf(f'{stat_dir}/rmse_croco_sevi_annual.nc',mode='w')

CPU times: user 23.8 s, sys: 9.01 s, total: 32.8 s
Wall time: 31.8 s