Newer
Older
#'@rdname CST_Analogs
#'@title Downscaling using Analogs based on large scale fields.
#'
#'@author Carmen Alvarez-Castro, \email{carmen.alvarez-castro@cmcc.it}
#'@author Nuria Perez-Zanon \email{nuria.perez@bsc.es}
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
#'@description This function perform a downscaling using Analogs. To compute
#'the analogs, the function search for days with similar large scale conditions
#'to downscaled fields in the local scale.
#'The large scale and the local scale regions are defined by the user.
#'The large scale is usually given by atmospheric circulation as sea level
#'pressure or geopotential height (Yiou et al, 2013) but the function gives the
#' possibility to use another field. The local scale will be usually given by
#' precipitation or Temperature, but might be another variable.
#' The analogs function will find the best analogs based in three criterias:
#' (1) Minimal distance in the large scale pattern (i.e. SLP)
#' (2) Minimal distance in the large scale pattern (i.e. SLP) and minimal
#' distance in the local scale pattern (i.e. SLP).
#' (3) Minimal distance in the large scale pattern (i.e. SLP), minimal
#' distance in the local scale pattern (i.e. SLP) and maxima correlation in the
#' local variable to downscale (i.e Precipitation).
#' The search of analogs must be done in the longest dataset posible. This is
#' important since it is necessary to have a good representation of the
#' possible states of the field in the past, and therefore, to get better
#' analogs. Once the search of the analogs is complete, and in order to used the
#' three criterias the user can select a number of analogs nAnalogs to restrict
#' the selection of the best analogs in a short number of posibilities, the best
#' ones. By default this parameter will be 1.
#' This function has not constrains of specific regions, variables to downscale,
#' or data to be used (seasonal forecast data, climate projections data,
#' reanalyses data).
#' The regrid into a finner scale is done interpolating with CST_Load.
#' Then, this interpolation is corrected selecting the analogs in the large
#' and local scale in based of the observations.
#' The function is an adapted version of the method of Yiou et al 2013.
#'@references Yiou, P., T. Salameh, P. Drobinski, L. Menut, R. Vautard,
#' and M. Vrac, 2013 : Ensemble reconstruction of the atmospheric column
#' from surface pressure using analogues. Clim. Dyn., 41, 1419-1437.
#' \email{pascal.yiou@lsce.ipsl.fr}
#'@param criteria different criteria to be used for the selection of analogs
#'if criteria = "Large_dist"
#'if criteria ="Local_dist"
#'if criteria ="Local_cor"
#'@param expL variable for the Large scale in the model (same variable
#'might be used in the local scale for criteria 2)
#'@param obsL variable for the large scale in the observations
#'@param expVar variable for the local scale in the model usually different
#'to the variable in expL
#'@param obsVar variable for the local scale in the observations usually
#'different to the variable in obsL
#'@param lon_local longitude in the local scale
#'@param lat_local latitude in the local scale
#'@param region region for the local scale
#'@param nAnalogs number of Analogs to be selected to apply the criterias (this
#'is not the necessary the number of analogs that the user can get)
#'@import ClimProjDiags
#'@return Dowscaled values of the best analogs in the criteria selected
CST_Analogs <- function(expL, obsL, expVar, obsVar, criteria = 'Large_dist') {
#checks
timevector <- obsL$Dates$start
region <- c(min(expVar$lon), max(expVar$lon), min(expVar$lat),
max(expVar$lon))
result <- Analogs(expL$data, obsL$data, time_obsL = timevector,
expVar = expVar$data, obsVar = obsVar$data,
criteria = criteria,
lon_local = expVar$lon, lat_local = expVar$lat,
region = region, nAnalogs = 1, return_list = FALSE)
obsVar$data <- result$AnalogsFields
result(obsVar)
}
#'@rdname Analogs
#'@title Search for analogs based on large scale fields.
#'
#'@author Carmen Alvarez-Castro, \email{carmen.alvarez-castro@cmcc.it}
#'@author Nuria Perez-Zanon \email{nuria.perez@bsc.es}
#'
#'@description This function search for days with similar large scale
#'conditions or similar large and local scale conditions.
83
84
85
86
87
88
89
90
91
92
93
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
#'The large scale and the local scale regions are defined by the user.
#'The large scale is usually given by atmospheric circulation as sea level
#'pressure or geopotential height (Yiou et al, 2013) but the function gives the
#' possibility to use another field. For the local scale the user can select
#' any variable.
#' The analogs function will find the best analogs based in three criterias:
#' (1) Minimal distance in the large scale pattern (i.e. SLP)
#' (2) Minimal distance in the large scale pattern (i.e. SLP) and minimal
#' distance in the local scale pattern (i.e. SLP).
#' (3) Minimal distance in the large scale pattern (i.e. SLP), minimal
#' distance in the local scale pattern (i.e. SLP) and maxima correlation in the
#' local variable to find the analog (i.e Precipitation).
#' Once the search of the analogs is complete, and in order to used the
#' three criterias the user can select a number of analogs nAnalogs to restrict
#' the selection of the best analogs in a short number of posibilities, the best
#' ones. By default this parameter will be 1.
#' This function has not constrains of specific regions, variables to find the
#' analogs, or data to be used (seasonal forecast data, climate projections
#' data, reanalyses data).
#' The input data might be interpolated or not.
#' The function is an adapted version of the method of Yiou et al 2013.
#'@references Yiou, P., T. Salameh, P. Drobinski, L. Menut, R. Vautard,
#' and M. Vrac, 2013 : Ensemble reconstruction of the atmospheric column
#' from surface pressure using analogues. Clim. Dyn., 41, 1419-1437.
#' \email{pascal.yiou@lsce.ipsl.fr}
#'@param criteria different criteria to be used for the selection of analogs
#'if criteria = "Large_dist"
#'if criteria ="Local_dist"
#'if criteria ="Local_cor"
#'@param expL variable for the Large scale in the model (same variable
#'might be used in the local scale for criteria 2)
#'@param obsL variable for the large scale in the observations
#'@param expVar variable for the local scale in the model usually different
#'to the variable in expL
#'@param obsVar variable for the local scale in the observations usually
#'different to the variable in obsL
#'@param lon_local longitude in the local scale
#'@param lat_local latitude in the local scale
#'@param region region for the local scale
#'@param nAnalogs number of Analogs to be selected to apply the criterias (this
#'is not the necessary the number of analogs that the user can get)
#'@param return_list TRUE if you want to get a list with the best analogs FALSE
#'if not.
#'@import ClimProjDiags
#'@return list list with the best analogs (time, distance)
#'@return values values of a certain variable
Analogs <- function(expL, obsL, time_obsL, expVar = NULL, obsVar = NULL,
criteria = "Large_dist",
lon_local = NULL, lat_local = NULL, region = NULL,
nAnalogs = 1, return_list = FALSE) {
# checks
if (!all(c('lon', 'lat') %in% names(dim(expL)))) {
stop("Parameter 'expL' must have the dimensions 'lat' and 'lon'.")
}
if (!all(c('lat', 'lon') %in% names(dim(obsL)))) {
stop("Parameter 'obsL' must have the dimension 'sdate'.")
}
if (any(is.na(expL))) {
warning("Parameter 'exp' contains NA values.")
}
if (any(is.na(obsL))) {
warning("Parameter 'obs' contains NA values.")
}
position <- Select(expL = expL, obsL = obsL, expVar = expVar,
obsVar = obsVar, criteria = criteria, lon_local = lon_local,
lat_local = lat_local,region = region)$position
best <- Apply(list(position), target_dims = c('time', 'pos'),
fun = BestAnalog, criteria = criteria, return_list = return_list,
nAnalogs = nAnalogs)
Analogs_dates <- time_obsL[best]
dim(Analogs_dates) <- dim(best)
if (is.null(obsVar)) {
obsLocal <- SelBox(obsL, lon = lon_local, lat = lat_local, region = region)
Analogs_fields <- Subset(obsLocal,
along = which(names(dim(obsLocal)) == 'time'), indices = best)
} else {
obsVar <- SelBox(obsVar, lon = lon_local, lat = lat_local, region = region)
Analogs_fields <- Subset(obsVar,along = which(names(dim(obsVar)) == 'time'),
indices = best)
}
lon_dim <- which(names(dim(Analogs_fields)) == 'lon')
lat_dim <- which(names(dim(Analogs_fields)) == 'lat')
if (lon_dim < lat_dim) {
dim(Analogs_fields) <- c(dim(Analogs_fields)[c(lon_dim, lat_dim)],
dim(best))
} else if (lat_dim > lon_dim) {
dim(Analogs_fields) <- c(dim(Analogs_fields)[c(lat_dim, lon_dim)],
dim(best))
} else {
stop("Dimensions 'lat' and 'lon' not found.")
}
return(list(DatesAnalogs = Analogs_dates, AnalogsFields = Analogs_fields))
BestAnalog <- function(position, criteria = 'Large_dist', return_list = FALSE,
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
pos_dim <- which(names(dim(position)) == 'pos')
if (dim(position)[pos_dim] == 1) {
pos1 <- position
if (criteria != 'Large_dist') {
warning("Dimension 'pos' in parameter 'position' has length 1,",
" criteria 'Large_dist' will be used.")
criteria <- 'Large_dist'
}
} else if (dim(position)[pos_dim] == 2) {
pos1 <- Subset(position, along = pos_dim, indices = 1)
pos2 <- Subset(position, along = pos_dim, indices = 2)
if (criteria == 'Local_cor') {
warning("Dimension 'pos' in parameter 'position' has length 2,",
" criteria 'Local_dist' will be used.")
criteria <- 'Local_dist'
}
} else if (dim(position)[pos_dim] == 3) {
pos1 <- Subset(position, along = pos_dim, indices = 1)
pos2 <- Subset(position, along = pos_dim, indices = 2)
pos3 <- Subset(position, along = pos_dim, indices = 3)
if (criteria != 'Local_cor') {
warning("Parameter 'criteria' is set to", criteria, ".")
}
} else {
stop("Parameter 'position' has dimension 'pos' of different ",
"length than expected (from 1 to 3).")
}
if (criteria == 'Large_dist') {
if (return_list == FALSE) {
pos <- pos1[1]
} else {
pos <- pos1[1 : nAnalogs]
}
} else if (criteria== 'Local_dist') {
# pos1 <- c(7, 13, 5, 3, 6, 12, 10, 1, 8, 9, 11, 4, 2, 14)
# pos2 <- c(4, 8, 13, 6, 3, 1, 12, 5, 9, 7, 10, 2, 11, 14)
pos1 <- pos1[1 : nAnalogs]
pos2 <- pos2[1 : nAnalogs]
best <- match(pos1, pos2)
pos <- pos1[as.logical(best)]
pos <- pos[which(!is.na(pos))]
if (return_list == FALSE) {
pos <- pos[1]
}
} else if (criteria == 'Local_cor') {
pos1 <- pos1[1 : nAnalogs]
pos2 <- pos2[1 : nAnalogs]
best <- match(pos1, pos2)
pos <- pos1[as.logical(best)]
pos <- pos[which(!is.na(pos))]
# pos3 <- c(6, 11, 14, 3, 13, 7, 2, 5, 1, 12, 10, 9, 8, 4)
best <- match(pos, pos3)
pos <- pos[order(best, decreasing = F)]
pos <- pos[which(!is.na(pos))]
if (return_list == FALSE) {
pos[1]
}
Select <- function(expL, obsL, expVar = NULL, obsVar = NULL,
criteria = "Large_dist",
lon_local = NULL, lat_local = NULL, region = NULL) {
#check expL
#check obsL
#check obsVar
names(dim(expL)) <- replace_repeat_dimnames(names(dim(expL)),names(dim(obsL)))
metric1 <- Apply(list(obsL), target_dims = list(c('lat', 'lon')),
fun = .select, expL, metric = "dist")$output1
dim_time_obs <- which(names(dim(metric1)) == 'time')
margins <- c(1 : length(dim(metric1)))[-dim_time_obs]
pos1 <- apply(metric1, margins, order)
names(dim(pos1))[1] <- 'time'
metric1 <- apply(metric1, margins, sort)
names(dim(metric1))[1] <- 'time'
dim(metric1) <- c(dim(metric1), metric = 1)
dim(pos1) <- c(dim(pos1), pos = 1)
return(list(metric = metric1, position = pos1))
}
if (criteria == "Local_dist" | criteria == "Local_cor") {
obs <- SelBox(obsL, lon = lon_local, lat = lat_local, region = region)$data
exp <- SelBox(expL, lon = lon_local, lat = lat_local, region = region)$data
metric2 <- Apply(list(obs), target_dims = list(c('lat', 'lon')),
fun = .select, exp, metric = "dist")$output1
margins <- c(1 : length(dim(metric2)))[-dim_time_obs]
pos2 <- apply(metric2, margins, order)
names(dim(pos2))[1] <- 'time'
metric2 <- apply(metric2, margins, sort)
names(dim(metric2))[1] <- 'time'
if (criteria == "Local_dist") {
metric <- abind(metric1, metric2, along = length(dim(metric1)) + 1)
position <- abind(pos1, pos2, along = length(dim(pos1)) + 1)
names(dim(metric)) <- c(names(dim(metric1)), 'metric')
names(dim(position)) <- c(names(dim(pos1)), 'pos')
return(list(metric = metric, position = position))
}
}
obs <- SelBox(obsVar, lon = lon_local, lat = lat_local,region = region)$data
exp <- SelBox(expVar, lon = lon_local, lat = lat_local,region = region)$data
metric3 <- Apply(list(obs), target_dims = list(c('lat', 'lon')),
fun = .select, exp, metric = "cor")$output1
margins <- c(1 : length(dim(metric3)))[-dim_time_obs]
pos3 <- apply(metric3, margins, order, decreasing = TRUE)
names(dim(pos3))[1] <- 'time'
metric3 <- apply(metric3, margins, sort)
names(dim(metric3))[1] <- 'time'
metric <- abind(metric1, metric2, metric3, along = length(dim(metric1)) + 1)
position <- abind(pos1, pos2, pos3, along = length(dim(pos1)) + 1)
names(dim(metric)) <- c(names(dim(metric1)), 'metric')
names(dim(position)) <- c(names(dim(pos1)), 'pos')
return(list(metric = metric, position = position))
}
else {
stop("Parameter 'criteria' must to be one of the: 'Large_dist', ",
"'Local_dist','Local_cor'.")
}
}
.select <- function(exp, obs, metric = "dist") {
if (metric == "dist") {
#metric <- sum((obs - exp) ^ 2)
#metric <- apply(obs, "time", function(x) {sum((x - exp) ^ 2)})
result <- Apply(list(obs), target_dims = list(c('lat', 'lon')),
fun = function(x) {sum((x - exp) ^ 2)})$output1
} else if (metric == "cor") {
result <- Apply(list(obs), target_dims = list(c('lat', 'lon')),
fun = function(x) {
cor(as.vector(x), as.vector(exp))})$output1
replace_repeat_dimnames <- function(names_exp, names_obs, lat_name = 'lat',
lon_name = 'lon') {
if (!is.character(names_exp)) {
stop("Parameter 'names_exp' must be a vector of characters.")
}
if (!is.character(names_obs)) {
stop("Parameter 'names_obs' must be a vector of characters.")
}
latlon_dim_exp <- which(names_exp == lat_name | names_exp == lon_name)
latlon_dim_obs <- which(names_obs == lat_name | names_obs == lon_name)
if (any(unlist(lapply(names_exp[-latlon_dim_exp],
function(x){x == names_obs[-latlon_dim_obs]})))) {
original_pos <- lapply(names_exp, function(x)
which(x == names_obs[-latlon_dim_obs]))
original_pos <- lapply(original_pos, length) > 0
names_exp[original_pos] <- paste0(names_exp[original_pos], "_exp")
}
return(names_exp)
## Improvements: other dimensions to avoid replacement for more flexibility.