#---------------------------------------------------------------------
# This script tells you how to load experimental and observational data in a 
# consistent way, facilating the following comparison.

# First, we load the experimental data. Because the latitude order of observation
# is opposite with experiment, and the sdate/time dimension is also different, we 
# use the attributes (sdate and latitude) of experimental data to define the 
# selectors for observation.

# You can see how to use parameter '*_across', 'merge_across_dims', and
# 'split_multiselected_dims' to create the consistent dimension as experiment.

#---------------------------------------------------------------------
library(startR)

# exp
repos_exp <- paste0('/esarchive/exp/ecearth/a1tr/cmorfiles/CMIP/EC-Earth-Consortium/',
                    'EC-Earth3/historical/r24i1p1f1/Amon/$var$/gr/v20190312/',
                    '$var$_Amon_EC-Earth3_historical_r24i1p1f1_gr_$sdate$01-$sdate$12.nc')

exp <- Start(dat = repos_exp,
             var = 'tas',
             sdate = as.character(c(2005:2008)),
             time = indices(1:3),
             lat = 'all',
             lon = 'all',
             synonims = list(lat = c('lat', 'latitude'),
                             lon = c('lon', 'longitude')),
             return_vars = list(lon = NULL,
                                lat = NULL,
                                time = 'sdate'),
             retrieve = FALSE)

# Retrieve attributes for the following observation.
# Because latitude order in experiment is [-90, 90] but in observation is [90, -90],
# latitude values need to be retrieved and used below.
lats <- attr(exp, 'Variables')$common$lat
# The 'time' attribute is dependent on 'sdate'. You can see the dimension below.
dates <- attr(exp, 'Variables')$common$time
# dim(dates)
#sdate time 
#    4    3 

#-------------------------------------------

# obs
# 1. For lat, use experiment attribute. For lon, it is not necessary because they have
# same values.
# 2. For dimension 'date', it is a vector involving the first 3 months (ftime) of the four years (sdate).
# 3. Dimension 'time' is assigned by the matrix, so we can seperate 'sdate' and 'time' 
# using 'split_multiselected_dims' later.
# 4. Because the 'time' is actually across all the files, so we need to specify 
# 'time_across'. Then, use 'merge_across_dims' to make dimension 'date' disappears. 
# At this moment, the dimension is 'time = 12'. 
# 5. However, we want to seperate year and month (which are 'sdate' and 'ftime' in 
# experimental data). So we use 'split_multiselected_dims' to split the two dimensions 
# of dimension 'time'.

repos_obs <- '/esarchive/recon/ecmwf/erainterim/monthly_mean/$var$_f6h/$var$_$date$.nc'

obs <- Start(dat = repos_obs,
             var = 'tas',
             date = unique(format(dates, '%Y%m')),
             time = values(dates),  #dim: [sdate = 4, time = 3]
             lat = values(lats),
             lon = 'all',
             time_across = 'date',
             merge_across_dims = TRUE,
             split_multiselected_dims = TRUE,
             synonims = list(lat = c('lat', 'latitude'),
                             lon = c('lon', 'longitude')),
             return_vars = list(lon = NULL,
                                lat = NULL,
                                time = 'date'),
             retrieve = FALSE)

#==========================
# Check attributes
#==========================

##-----dimension-----
print(attr(exp, 'Dimensions'))
#  dat   var sdate  time   lat   lon 
#    1     1     4     3   256   512 

print(attr(obs, 'Dimensions'))
#  dat   var sdate  time   lat   lon 
#    1     1     4     3   256   512 

##-----time-----
print(attr(exp, 'Variables')$common$time)
# [1] "2005-01-16 13:14:44 CET" "2006-01-16 13:14:44 CET"
# [3] "2007-01-16 13:14:44 CET" "2008-01-16 13:14:44 CET"
# [5] "2005-02-15 01:14:44 CET" "2006-02-15 01:14:44 CET"
# [7] "2007-02-15 01:14:44 CET" "2008-02-15 13:14:44 CET"
# [9] "2005-03-16 13:14:44 CET" "2006-03-16 13:14:44 CET"
#[11] "2007-03-16 13:14:44 CET" "2008-03-16 13:14:44 CET"

print(attr(obs, 'Variables')$common$time)
# [1] "2005-01-31 18:00:00 CET"  "2006-01-31 18:00:00 CET" 
# [3] "2007-01-31 18:00:00 CET"  "2008-01-31 18:00:00 CET" 
# [5] "2005-02-28 18:00:00 CET"  "2006-02-28 18:00:00 CET" 
# [7] "2007-02-28 18:00:00 CET"  "2008-02-29 18:00:00 CET" 
# [9] "2005-03-31 19:00:00 CEST" "2006-03-31 19:00:00 CEST"
#[11] "2007-03-31 19:00:00 CEST" "2008-03-31 19:00:00 CEST"

##-----lat-----
print(attr(exp, 'Variables')$common$lat[1:3])
#[1] -89.46282 -88.76695 -88.06697
print(attr(exp, 'Variables')$common$lat[256])
#[1] 89.46282

print(attr(obs, 'Variables')$common$lat[1:3])
#[1] -89.46282 -88.76695 -88.06697
print(attr(obs, 'Variables')$common$lat[256])
#[1] 89.46282