CST_Analogs.R

    best <- match(pos1, pos2)
    if (length(best) == 1) {
      warning("Just 1 best analog matching Large_dist and ", 
              "Local_dist criteria")
    } 
    if (length(best) < 1 | is.na(best[1]) == TRUE) {
      stop("no best analogs matching Large_dist and Local_dist criterias, 
            please increase nAnalogs") 
    }
    pos <- pos2[as.logical(best)]
    pos <- pos[which(!is.na(pos))]
    if (AnalogsInfo == FALSE) { 
      pos <- pos[1]
    }else {
      pos <- pos}
  } else if (criteria == 'Local_cor') {
    pos1 <- pos1[1 : nAnalogs]
    pos2 <- pos2[1 : nAnalogs]
    best <- match(pos1, pos2)
    if (length(best) == 1) {
      warning("Just 1 best analog matching Large_dist and ", 
              "Local_dist criteria")
    } 
    if (length(best) < 1 | is.na(best[1]) == TRUE) {
      stop("no best analogs matching Large_dist and Local_dist criterias, 
            please increase nAnalogs") 
    }
    pos <- pos1[as.logical(best)]
    pos <- pos[which(!is.na(pos))]
    pos3 <- pos3[1 : nAnalogs]
    best <- match(pos, pos3)
    if (length(best) == 1) {
      warning("Just 1 best analog matching Large_dist, Local_dist and ", 
              "Local_cor criteria")
    } 
    if (length(best) < 1 | is.na(best[1]) == TRUE) {
      stop("no best analogs matching Large_dist, Local_dist and Local_cor 
            criterias, please increase nAnalogs") 
    }
    pos <- pos[order(best, decreasing = F)]
    pos <- pos[which(!is.na(pos))]
    if (AnalogsInfo == FALSE) { 
      pos <- pos[1]
    } else{
      pos <- pos
    }
    return(pos)
  }
}
Select <- function(expL, obsL,  expVar = NULL, obsVar = NULL, 
                   criteria = "Large_dist", lonL = NULL, latL = NULL,
                   lonVar = NULL, latVar = NULL, region = NULL, 
                   lon_name = 'lon', lat_name = 'lat') {
  names(dim(expL)) <- replace_repeat_dimnames(names(dim(expL)), 
                                              names(dim(obsL)),
                                              lon_name = lon_name,
                                              lat_name = lat_name)
  metric1 <- Apply(list(obsL), target_dims = list(c(lat_name, lon_name)), 
                   fun = .select, expL, metric = "dist",
                   lon_name = lon_name, lat_name = lat_name)$output1
  metric1.original = metric1
  if (length(dim(metric1)) > 1) {
    dim_time_obs <- which(names(dim(metric1)) == 'time' | 
                          names(dim(metric1)) == 'ftime')
    dim(metric1) <- c(dim(metric1), metric=1)
    margins <- c(1 : (length(dim(metric1))))[-dim_time_obs]
    pos1 <- apply(metric1, margins, order)      
    names(dim(pos1))[1] <- 'time'
    metric1.original = metric1
    metric1 <-  apply(metric1, margins, sort)
    names(dim(metric1))[1] <- 'time'
    names(dim(metric1.original)) = names(dim(metric1))
  } else {
    pos1 <- order(metric1)
    dim(pos1) <- c(time = length(pos1))
    metric1 <- sort(metric1)
    dim(metric1) <- c(time = length(metric1))
    dim(metric1.original) = dim(metric1)
    dim_time_obs = 1
  }
  if (criteria == "Large_dist") {
    dim(metric1) <- c(dim(metric1), metric = 1)
    dim(pos1) <- c(dim(pos1), pos = 1)
    dim(metric1.original) = dim(metric1)
    return(list(metric = metric1, metric.original = metric1.original,
                position = pos1))
  }
  if (criteria == "Local_dist" | criteria == "Local_cor") {
    obs <- SelBox(obsL, lon = lonL, lat = latL, region = region)$data
    exp <- SelBox(expL, lon = lonL, lat = latL, region = region)$data
    metric2 <- Apply(list(obs), target_dims = list(c(lat_name, lon_name)), 
                     fun = .select, exp, metric = "dist", 
                     lon_name = lon_name, lat_name = lat_name)$output1
    metric2.original = metric2
    dim(metric2) <- c(dim(metric2), metric=1)
    margins <- c(1 : (length(dim(metric2))))[-dim_time_obs]
    pos2 <- apply(metric2, margins, order)
    dim(pos2) <- dim(pos1)
    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)
      metric.original <- abind(metric1.original,metric2.original,
                               along = length(dim(metric1))+1)
      position <- abind(pos1, pos2, along = length(dim(pos1))+1) 
      names(dim(metric)) <- c(names(dim(pos1)), 'metric')
      names(dim(position)) <- c(names(dim(pos1)), 'pos')
      names(dim(metric.original)) = names(dim(metric)) 
      return(list(metric = metric, metric.original = metric.original,
                  position = position))
    }   
  }
  if (criteria == "Local_cor") {
    obs <- SelBox(obsVar, lon = lonVar, lat = latVar, region = region)$data
    exp <- SelBox(expVar, lon = lonVar, lat = latVar, region = region)$data
    metric3 <- Apply(list(obs), target_dims = list(c(lat_name, lon_name)), 
                     fun = .select, exp, metric = "cor",
                     lon_name = lon_name, lat_name = lat_name)$output1
    metric3.original = metric3
    dim(metric3) <- c(dim(metric3), metric=1)
    margins <- c(1 : (length(dim(metric3))))[-dim_time_obs]
    pos3 <- apply(abs(metric3), margins, order, decreasing = TRUE)
    names(dim(pos3))[1] <- 'time'
    metricsort <- metric3[pos3]
    dim(metricsort) = dim(metric3)
    names(dim(metricsort))[1] <- 'time'
    metric <- abind(metric1, metric2, metricsort, 
                    along = length(dim(metric1)) + 1)
    metric.original <- abind(metric1.original, metric2.original, 
                             metric3.original, 
                             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')
    names(dim(metric.original)) = names(dim(metric)) 
    return(list(metric = metric, metric.original=metric.original,
                position = position))
  }
  else {
    stop("Parameter 'criteria' must to be one of the: 'Large_dist', ",
         "'Local_dist','Local_cor'.")
  }
}
.select <- function(exp, obs, metric = "dist", 
                    lon_name = 'lon', lat_name = 'lat') {
  if (metric == "dist") {
    result <- Apply(list(obs), target_dims = list(c(lat_name, lon_name)), 
            fun = function(x) {sqrt(sum((x - exp) ^ 2, na.rm = TRUE))})$output1
  } else if (metric == "cor") {
    result <- Apply(list(obs), target_dims = list(c(lat_name, lon_name)), 
                    fun = function(x) {cor(as.vector(x), 
                                           as.vector(exp),
                                           method = "spearman")})$output1
  } 
  result
}
.time_ref <- function(time_obsL,time_expL,excludeTime) {
  sameTime = which(time_obsL %in% time_expL)
  result<- c(time_obsL[1:(sameTime - excludeTime - 1)],
             time_obsL[(sameTime + excludeTime + 1):length(time_obsL)])
  result
}

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)
}

replace_time_dimnames <- function(dataL, time_name = 'time', 
                                  stdate_name = 'stdate', ftime_name='ftime') {
  names_obs = names(dim(dataL))
  if (!is.character(names_obs)) {
    stop("Parameter 'names_obs' must be a vector of characters.")
  }
  time_dim_obs <- which(names_obs == time_name | 
                        names_obs == stdate_name | names_obs == ftime_name)
  if (length(time_dim_obs) > 1) {
    stop ("more than 1 time dimension, please give just 1")
  }
  if (length(time_dim_obs) == 0) {
    warning ("name of time dimension is not 'ftime' or 'time' or 'stdate' 
              or time dimension is null")
  }
  if (length(time_dim_obs) != 0) {
    names_obs[time_dim_obs]= time_name
  }
  names(dim(dataL)) = names_obs
  return(dataL)
}