Hi @eexarchou, @aho,
I have fixed the bug in AbsBiasSS() mentioned in the issue #80 (closed).
The significance is now computed before the time averaging. I have used the code you provided in the issue to make sure that the bug is fixed. I had to change lonmax <- 180
to lonmax <- 179
because only one longitude was loaded otherwise. Also, I have assumed that a511.mjj
, a511.aso
, a511.ndj
and a511.fma
are seasonal means and that the experiments are initialised in May (please check the lines just before computing the skill score, where s2dv::Season is used):
# Delete previous data
rm(list=ls())
gc()
library(boot)
library(s2dv)
library(abind)
library(s2dverification)
library(startR)
library(ncdf4)
library(multiApply)
library(SpecsVerification)
#source('AbsBiasSS.R')
#source('Bias.R')
#source('RandomWalkTest.R')
# ATL
lonmin <- -180
lonmax <- 179
latmin <- -60
latmax <- 50
sdates <- paste0(c(1981:2002), '0501')
# exp
repos_a511 <- paste0('/esarchive/scratch/eexarcho/Eleftheria/Pace_maker/Data/a511/tos/r360x180/',
'$var$_Omon_EC-Earth3-CC_troppac-pacemaker_s$sdate$-$member$_r360x180_$chunk$.nc')
repos_a53d <- paste0('/esarchive/scratch/eexarcho/Eleftheria/Pace_maker/Data/a53d/tos/r360x180/',
'$var$_Omon_EC-Earth3-CC_troppac-pacemaker_s$sdate$-$member$_r360x180_$chunk$.nc')
a511 <- Start(dat = repos_a511,
var = 'tos',
member = 'all',
sdate = sdates,
chunk = 'all',
# time = 'all',
time = indices(1:12), #first time step per day
chunk_depends = 'sdate',
time_across = 'chunk',
merge_across_dims = TRUE,
lat = values(list(latmin, latmax)),
lat_reorder = Sort(), ##1
lon = values(list(lonmin, lonmax)),
lon_reorder = CircularSort(-180, 180), ##2
synonims = list(lat = c('lat', 'latitude'),
lon = c('lon', 'longitude')),
return_vars = list(lon = 'dat', ##3
lat = 'dat', ##3
time = 'sdate'),
retrieve = T)
a53d <- Start(dat = repos_a53d,
var = 'tos',
member = 'all',
sdate = sdates,
chunk = 'all',
# time = 'all',
time = indices(1:12), #first time step per day
chunk_depends = 'sdate',
time_across = 'chunk',
merge_across_dims = TRUE,
lat = values(list(latmin, latmax)),
lat_reorder = Sort(), ##1
lon = values(list(lonmin, lonmax)),
lon_reorder = CircularSort(-180, 180), ##2
synonims = list(lat = c('lat', 'latitude'),
lon = c('lon', 'longitude')),
return_vars = list(lon = 'dat', ##3
lat = 'dat', ##3
time = 'sdate'),
retrieve = T)
###lons <- attr(exp, 'Variables')$common$tos$dim[[1]]$vals
###lats <- attr(exp, 'Variables')$common$tos$dim[[2]]$vals
lons <- as.vector(attr(a511, 'Variables')$dat1$lon)
lats <- as.vector(attr(a511, 'Variables')$dat1$lat)
dates <- attr(a511, 'Variables')$common$time
dates_file <- sort(unique(gsub('-', '', sapply(as.character(dates), substr, 1, 7))))
##
## obs
repos_obs <- '/esarchive/obs/ukmo/hadisst_v1.1/monthly_mean/$var$/$var$_$date$.nc'
obs.data <- Start(dat = repos_obs,
var = 'tos',
date = dates_file,
# time = 'all',
# time_across = 'date',
# merge_across_dims = TRUE,
# merge_across_dims_narm = TRUE,
## -----------------------------------------------------------------------
time = values(dates), #dim: [sdate = 2, time = 12]
#because time is assigned by 'values', set the tolerance to avoid too distinct match
time_var = 'time',
# time_tolerance = as.difftime(372, units = 'hours'),
#time values are across all the files
time_across = 'date',
#combine time and file_date dims
merge_across_dims = TRUE,
#exclude the additional NAs generated by merge_across_dims
merge_across_dims_narm = TRUE,
#split time dim, because it is two-dimensional
split_multiselected_dims = TRUE,
## -----------------------------------------------------------------------
lat = values(lats),
lon = values(lons),
#---------transform to get identical lat and lon------------ ##4
transform = CDORemapper,
transform_extra_cells = 2,
transform_params = list(grid = 'r360x180',
method = 'conservative',
crop = c(lonmin, lonmax, latmin, latmax)),
transform_vars = c('lat', 'lon'),
#-----------------------------------------------------------
synonims = list(lat = c('lat', 'latitude'),
lon = c('lon', 'longitude')),
return_vars = list(lat = NULL, ##3
lon = NULL, ##3
time = 'date'),
retrieve = TRUE)
##### I have added this
a511.mjj <- s2dv::Season(data = a511, time_dim = 'time', monini = 5, moninf = 5, monsup = 7, method = mean, na.rm = F, ncores = 4)
a511.aso <- s2dv::Season(data = a511, time_dim = 'time', monini = 5, moninf = 8, monsup = 10, method = mean, na.rm = F, ncores = 4)
a511.ndj <- s2dv::Season(data = a511, time_dim = 'time', monini = 5, moninf = 11, monsup = 1, method = mean, na.rm = F, ncores = 4)
a511.fma <- s2dv::Season(data = a511, time_dim = 'time', monini = 5, moninf = 2, monsup = 4, method = mean, na.rm = F, ncores = 4)
a53d.mjj <- s2dv::Season(data = a53d, time_dim = 'time', monini = 5, moninf = 5, monsup = 7, method = mean, na.rm = F, ncores = 4)
a53d.aso <- s2dv::Season(data = a53d, time_dim = 'time', monini = 5, moninf = 8, monsup = 10, method = mean, na.rm = F, ncores = 4)
a53d.ndj <- s2dv::Season(data = a53d, time_dim = 'time', monini = 5, moninf = 11, monsup = 1, method = mean, na.rm = F, ncores = 4)
a53d.fma <- s2dv::Season(data = a53d, time_dim = 'time', monini = 5, moninf = 2, monsup = 4, method = mean, na.rm = F, ncores = 4)
obs.mjj <- s2dv::Season(data = obs.data, time_dim = 'time', monini = 5, moninf = 5, monsup = 7, method = mean, na.rm = F, ncores = 4)
obs.aso <- s2dv::Season(data = obs.data, time_dim = 'time', monini = 5, moninf = 8, monsup = 10, method = mean, na.rm = F, ncores = 4)
obs.ndj <- s2dv::Season(data = obs.data, time_dim = 'time', monini = 5, moninf = 11, monsup = 1, method = mean, na.rm = F, ncores = 4)
obs.fma <- s2dv::Season(data = obs.data, time_dim = 'time', monini = 5, moninf = 2, monsup = 4, method = mean, na.rm = F, ncores = 4)
source('https://earth.bsc.es/gitlab/es/s2dv/-/raw/develop-AbsBiasSS-sign/R/AbsBiasSS.R')
.Bias <- s2dv:::.Bias
.RandomWalkTest <- s2dv:::.RandomWalkTest
##### Until here
abs.bias.a511.mjj <- AbsBiasSS(a511.mjj, obs.mjj-273.15, a53d.mjj,
memb_dim = 'member', dat_dim='dat', ncores=16)
abs.bias.a511.aso <- AbsBiasSS(a511.aso, obs.aso-273.15, a53d.aso,
memb_dim = 'member', dat_dim='dat', ncores=16)
abs.bias.a511.ndj <- AbsBiasSS(a511.ndj, obs.ndj-273.15, a53d.ndj,
memb_dim = 'member', dat_dim='dat', ncores=16)
abs.bias.a511.fma <- AbsBiasSS(a511.fma, obs.fma-273.15, a53d.fma,
memb_dim = 'member', dat_dim='dat', ncores=16)
min2=-2.
max2=2.
int2=(max2-min2)/20
interval2=seq(min2,max2,int2)
clim.palette(palette = "bluered")
color=clim.colors(20, palette = "bluered")
pdf( "~/Desktop/SST_bias_SS_a511_a53d.pdf" )
nf <- layout(matrix(c(1, 2 , 3, 4, 5, 5, 5, 5),2,4,byrow=TRUE),
widths=1, heights=c(1,0.3), TRUE)
s2dv::PlotEquiMap( abs.bias.a511.mjj$biasSS[1, 1, 1, 1, , ],
lons, lats ,
toptitle= "SST bias SS, MJJ",
cols=color, brks=interval2, drawleg = F,
triangle_ends = c(T,T), col_inf = color[1], col_sup = color[20],
filled.continents = T,
dots = abs.bias.a511.mjj$sign[1, 1, 1, 1, ,],
dot_size = 0.5,
numbfig=4
)
s2dv::PlotEquiMap( abs.bias.a511.aso$biasSS[1, 1, 1, 1, , ],
lons, lats ,
toptitle= "SST bias SS, ASO",
cols=color, brks=interval2, drawleg = F,
triangle_ends = c(T,T), col_inf = color[1], col_sup = color[20],
filled.continents = T,
dots = abs.bias.a511.aso$sign[1, 1, 1, 1, ,],
dot_size = 0.5,
numbfig=4
)
s2dv::PlotEquiMap( abs.bias.a511.ndj$biasSS[1, 1, 1, 1, , ],
lons, lats ,
toptitle= "SST bias SS, NDJ",
cols=color, brks=interval2, drawleg = F,
triangle_ends = c(T,T), col_inf = color[1], col_sup = color[20],
filled.continents = T,
dots = abs.bias.a511.ndj$sign[1, 1, 1, 1, ,],
dot_size = 0.5,
numbfig=4
)
s2dv::PlotEquiMap( abs.bias.a511.fma$biasSS[1, 1, 1, 1, , ],
lons, lats ,
toptitle= "SST bias SS, FMA",
cols=color, brks=interval2, drawleg = F,
triangle_ends = c(T,T), col_inf = color[1], col_sup = color[20],
filled.continents = T,
dots = abs.bias.a511.fma$sign[1, 1, 1, 1, ,],
dot_size = 0.5,
numbfig=4
)
ColorBar(brks=interval2,
triangle_ends = c(T,T), col_inf = color[1], col_sup = color[20],
cols = color, vert = FALSE)
dev.off()
Please let me know if you think the new version is correct or if any modifications are needed.
A suggestion: for plotting skill scores, I think a more suitable colorbar would be one from -0.9 to 0.9 with triangle_ends = c(T,T)
, or from -1 to 1 with triangle_ends = c(T,F)
. This is because the skill scores range from -Inf to 1. But, of course, everyone can decide different bar limits :)
Best regards,
Carlos