Newer
Older
save_metrics <- function(recipe,
data_cube,
agg = "global",
# This function adds metadata to the skill metrics in 'skill'
# and exports them to a netCDF file inside 'outdir'.
# Define grid dimensions and names
lalo <- c('longitude', 'latitude')
archive <- get_archive(recipe)
dictionary <- read_yaml("conf/variable-dictionary.yml")
global_attributes <- .get_global_attributes(recipe, archive)
## TODO: Sort out the logic once default behavior is decided
if ((!is.null(recipe$Analysis$Workflow$Anomalies$compute)) &&
(recipe$Analysis$Workflow$Anomalies$compute)) {
global_attributes <- c(list(from_anomalies = "Yes"),
global_attributes)
} else {
global_attributes <- c(list(from_anomalies = "No"),
global_attributes)
}
# Time indices and metadata
fcst.horizon <- tolower(recipe$Analysis$Horizon)
store.freq <- recipe$Analysis$Variables$freq
if (global_attributes$system == 'Multimodel'){
if (fcst.horizon == 'decadal'){
calendar <- archive$System[[recipe$Analysis$Datasets$System$models[[1]]$name]]$calendar
} else {
calendar <- archive$System[[recipe$Analysis$Datasets$System$models[[1]]]]$calendar
}
} else {
calendar <- archive$System[[global_attributes$system]]$calendar
}
# Generate vector containing leadtimes
dates <- as.PCICt(ClimProjDiags::Subset(data_cube$attrs$Dates, 'syear', 1),
cal = calendar)
if (fcst.horizon == 'decadal') {
# init_month <- archive$System[[recipe$Analysis$Datasets$System$name]]$initial_month
# init_date <- as.PCICt(paste0(recipe$Analysis$Time$hcst_start, '-',
# sprintf('%02d', init_month), '-01'),
# cal = calendar)
init_date <- as.PCICt(paste0(as.numeric(recipe$Analysis$Time$hcst_start)+1,
'-01-01'), cal = calendar)
} else {
init_date <- as.PCICt(paste0(recipe$Analysis$Time$hcst_start,
recipe$Analysis$Time$sdate),
format = '%Y%m%d', cal = calendar)
}
# Get time difference in hours
leadtimes <- as.numeric(dates - init_date)/3600
# Select start date
# If a fcst is provided, use that as the ref. year. Otherwise use 1970.
if (fcst.horizon == 'decadal') {
if (!is.null(recipe$Analysis$Time$fcst_year)) {
#PROBLEM: May be more than one fcst_year
# fcst.sdate <- paste0(recipe$Analysis$Time$fcst_year[1],
# sprintf('%02d', init_month), '01')
fcst.sdate <- paste0(recipe$Analysis$Time$fcst_year[1])
} else {
fcst.sdate <- paste0("1970", sprintf('%02d', init_month), '01')
}
} else {
if (!is.null(recipe$Analysis$Time$fcst_year)) {
fcst.sdate <- paste0(recipe$Analysis$Time$fcst_year,
recipe$Analysis$Time$sdate)
} else {
fcst.sdate <- paste0("1970", recipe$Analysis$Time$sdate)
}
}
times <- .get_times(store.freq, fcst.horizon, leadtimes, fcst.sdate, calendar)
time <- times$time
# Loop over variable dimension
for (var in 1:data_cube$dims[['var']]) {
# Subset skill arrays
subset_metric <- lapply(metrics, function(x) {
ClimProjDiags::Subset(x, along = 'var',
indices = var,
drop = 'selected')})
# Generate name of output file
variable <- data_cube$attrs$Variable$varName[[var]]
outdir <- get_dir(recipe = recipe, variable = variable)
if (!dir.exists(outdir)) {
dir.create(outdir, recursive = T)
}
if (tolower(recipe$Analysis$Output_format) == "scorecards") {
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
for (i in 1:length(subset_metric)) {
if (any('syear' %in% names(dim(subset_metric[[i]])))) {
sdate_dim_save = 'syear'
dates <- data_cube$attrs$Dates
} else {
sdate_dim_save = NULL
dates <- Subset(data_cube$attrs$Dates, along = 'syear', indices = 1)
}
## TODO: Maybe 'scorecards' condition could go here to further simplify
## the code
extra_string <- get_filename(NULL, recipe, variable,
fcst.sdate, agg, names(subset_metric)[[i]])
SaveExp(data = subset_metric[[i]], destination = outdir,
Dates = dates,
coords = c(data_cube$coords['longitude'],
data_cube$coords['latitude']),
varname = names(subset_metric)[[i]],
metadata = data_cube$attrs$Variable$metadata, Datasets = NULL,
startdates = NULL, dat_dim = NULL, sdate_dim = sdate_dim_save,
ftime_dim = 'time', var_dim = NULL, memb_dim = NULL,
drop_dims = NULL, single_file = TRUE,
extra_string = extra_string)
}
} else {
outfile <- get_filename(outdir, recipe, variable,
fcst.sdate, agg, module)
# Remove singleton dimensions and rearrange lon, lat and time dims
if (tolower(agg) == "global") {
subset_metric <- lapply(subset_metric, function(x) {
Reorder(x, c(lalo, 'time'))})
}
attr(subset_metric[[1]], 'global_attrs') <- global_attributes
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
for (i in 1:length(subset_metric)) {
metric <- names(subset_metric[i])
long_name <- dictionary$metrics[[metric]]$long_name
missing_val <- -9.e+33
subset_metric[[i]][is.na(subset_metric[[i]])] <- missing_val
if (tolower(agg) == "country") {
sdname <- paste0(metric, " region-aggregated metric")
dims <- c('Country', 'time')
} else if (tolower(agg) == "region") {
sdname <- paste0(metric, " region-aggregated metric")
dims <- c('region', 'time')
} else {
sdname <- paste0(metric)
dims <- c(lalo, 'time')
}
metadata <- list(metric = list(name = metric,
standard_name = sdname,
long_name = long_name,
missing_value = missing_val))
attr(subset_metric[[i]], 'variables') <- metadata
names(dim(subset_metric[[i]])) <- dims
}
# Get grid data and metadata and export to netCDF
if (tolower(agg) == "country") {
country <- get_countries(grid)
ArrayToNc(append(country, time, subset_metric), outfile)
region <- array(1:dim(metrics[[1]])['region'], c(dim(metrics[[1]])['region']))
# TODO: check metadata when more than 1 region is store in the data array
metadata <- list(region = list(long_name = data_cube$attrs$Variable$metadata$region$name))
attr(region, 'variables') <- metadata
vars <- list(region, time)
vars <- c(vars, subset_metric)
ArrayToNc(vars, outfile)
latitude <- data_cube$coords$lat[1:length(data_cube$coords$lat)]
longitude <- data_cube$coords$lon[1:length(data_cube$coords$lon)]
latlon <- .get_latlon(latitude, longitude)
# Compile variables into a list and export to netCDF
vars <- list(latlon$lat, latlon$lon, time)
vars <- c(vars, subset_metric)
ArrayToNc(vars, outfile)
}
info(recipe$Run$logger,
paste("#####", toupper(module), " METRICS SAVED TO NETCDF FILE #####"))
}