R/helpers.R

In SantanderMetGroup/transformeR: A climate4R package for general climate data manipulation and transformation

##     helpers.R Helper functions of the climate4R bundle
##
##     Copyright (C) 2018 Santander Meteorology Group (http://www.meteo.unican.es)
##
##     This program is free software: you can redistribute it and/or modify
##     it under the terms of the GNU General Public License as published by
##     the Free Software Foundation, either version 3 of the License, or
##     (at your option) any later version.
## 
##     This program is distributed in the hope that it will be useful,
##     but WITHOUT ANY WARRANTY; without even the implied warranty of
##     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
##     GNU General Public License for more details.
## 
##     You should have received a copy of the GNU General Public License
##     along with this program. If not, see <http://www.gnu.org/licenses/>.


#' @title Get regular grid definition 
#' @description Get the (regular) grid definition from an existing (gridded) dataset
#' @param gridData A grid data object coming from \code{loadGridData} (package \pkg{loadeR}) or \code{\link{interpGrid}}
#'  or the function \code{loadECOMS} of package \pkg{loadeR.ECOMS}.
#' @return A list of two named components, \code{x} and \code{y}, consisting of a vector of length two each one, defining
#' the x/y lower and upper bounds. The grid-cell resolution is given by the attributes \code{'resX'} and
#'  \code{'resY'} respectively.
#' @details In case of irregular grid definitions, the function forces the grid to regularity.
#' The returned grid object inherits the attributes from the input \code{xyCoords} definition.
#' @export
#' @importFrom utils tail
#' @family get.helpers
#' @author S. Herrera and J. Bedia 
#' @examples \donttest{ 
#' require(climate4R.datasets)
#' # Regular grid
#' data("NCEP_Iberia_hus850")
#' str(getGrid(NCEP_Iberia_hus850))
#' # Regular rotated grid with correpondence with a non regular grid in the 
#' # lat/lon domain.
#' data("CORDEX_Iberia_tas")
#' str(getGrid(CORDEX_Iberia_tas))
#' # Station data
#' data("VALUE_Iberia_tas")
#' str(getGrid(VALUE_Iberia_tas))
#' 
#' # Interpolate NCEP onto the System4 grid:
#' data("CFS_Iberia_pr")
#' NCEP_interpolated <- interpGrid(NCEP_Iberia_hus850, getGrid(CFS_Iberia_pr), "bilinear")
#' # Note the warnings because of the non-overlapping domain extents (longitudes)
#' require(visualizeR)
#' spatialPlot(climatology(NCEP_Iberia_hus850), backdrop.theme = "countries")
#' spatialPlot(climatology(NCEP_interpolated), backdrop.theme = "countries")
#' str(getGrid(NCEP_interpolated))
#' str(getGrid(CFS_Iberia_pr))
#' 
#' # The other way round, using nearest neighbour interpolation:
#' CFS_interpolated <- interpGrid(CFS_Iberia_pr, getGrid(NCEP_Iberia_hus850))
#' spatialPlot(climatology(CFS_interpolated), backdrop.theme = "countries")
#' # In this case, the mismatch in domain extent occurs in the longitudes and latitudes
#' }
 

getGrid <- function(gridData) {
      rot <- FALSE
      if (is.data.frame(gridData$xyCoords)) {
            out <- list(x = as.numeric(gridData$xyCoords[,1]), y = as.numeric(gridData$xyCoords[,2]))
            attr(out, "type") <- "irregular"
            if (!exists("resX", attributes(gridData$xyCoords))) {
                  attr(out, "resX") <- NULL
            } else {
                  attr(out, "resX") <- attr(gridData$xyCoords, "resX")
            }
            if (!exists("resY", attributes(gridData$xyCoords))) {
                  attr(out, "resY") <- NULL
            } else {
                  attr(out, "resY") <- attr(gridData$xyCoords, "resY")   
            }
      }else{
            if ("lon" %in% names(gridData$xyCoords)) rot <- TRUE
            if (length(gridData$xyCoords$x) > 1){
                  grid.x <- c(gridData$xyCoords$x[1], tail(gridData$xyCoords$x, 1))
            }else{
                  grid.x <- gridData$xyCoords$x
            }
            if (length(gridData$xyCoords$y) > 1){
                  grid.y <- c(gridData$xyCoords$y[1], tail(gridData$xyCoords$y, 1))
            }else{
                  grid.y <- gridData$xyCoords$y
            }
            ## grid.x <- c(gridData$xyCoords$x[1], tail(gridData$xyCoords$x, 1))
            ## grid.y <- c(gridData$xyCoords$y[1], tail(gridData$xyCoords$y, 1))
            out <- list(x = grid.x, y = grid.y)
            if (rot) {
                  out$lon <- gridData$xyCoords$lon 
                  out$lat <- gridData$xyCoords$lat
            }
            attributes(out) <- attributes(gridData$xyCoords)
            if (!exists("resX", attributes(gridData$xyCoords))) {
                  attr(out, "resX") <- (tail(gridData$xyCoords$x, 1) - gridData$xyCoords$x[1]) / (length(gridData$xyCoords$x) - 1)
            } else {
                  attr(out, "resX") <- attr(gridData$xyCoords, "resX")
            }
            if (!exists("resY", attributes(gridData$xyCoords))) {
                  attr(out, "resY") <- (tail(gridData$xyCoords$y, 1) - gridData$xyCoords$y[1]) / (length(gridData$xyCoords$y) - 1)
                  
            } else{ 
                  attr(out, "resY") <- attr(gridData$xyCoords, "resY")   
            }
            if (rot) {
                  attr(out, "resLON") <- NA 
                  attr(out, "resLAT") <- NA
            } 
      }
      return(out)
}
# End


#' @title Get geographical coordinates of a climate data object
#' @description Returns the coordinates of a climate data object, either stations or grid
#' @param obj Any object extending the station or grid classes
#' @return A list with x and y components
#' @family get.helpers
#' @author J. Bedia 
#' @export
#' @examples \donttest{ 
#' require(climate4R.datasets)
#' # Regular grid
#' data("NCEP_Iberia_hus850")
#' str(getCoordinates(NCEP_Iberia_hus850))
#' # Regular rotated grid with correpondence with a non regular grid in the 
#' # lat/lon domain.
#' data("CORDEX_Iberia_tas")
#' str(getCoordinates(CORDEX_Iberia_tas))
#' # Station data
#' data("VALUE_Iberia_tas")
#' str(getCoordinates(VALUE_Iberia_tas))
#' }


getCoordinates <- function(obj) {
    if (is.data.frame(obj$xyCoords)) {
        xy <- obj$xyCoords#[ , ,drop = FALSE]
        return(xy)
    } else {
        x <- obj$xyCoords$x
        y <- obj$xyCoords$y
        if (!is.data.frame(obj$xyCoords) && exists("lon", where = obj$xyCoords) && exists("lat", where = obj$xyCoords)) {
            lon <- obj$xyCoords$lon 
            lat <- obj$xyCoords$lat
            return(list("x" = x, "y" = y, "lon" = lon, "lat" = lat))
        } else {
            return(list("x" = x, "y" = y))
        }
    }
}
# End



#' @title Get season 
#' @description Retrieves the season encompassed by a station or grid object
#' @param obj Any object extending the station or grid classes
#' @return An integer vector with the season
#' @author J. Bedia 
#' @family get.helpers
#' @export
#' @examples \donttest{ 
#' require(climate4R.datasets)
#' data("NCEP_Iberia_hus850")
#' getSeason(NCEP_Iberia_hus850) # Boreal winter (DJF)
#' }

getSeason <- function(obj) {
      if ("season" %in% names(attributes(obj$Dates))) {
            attr(obj$Dates, "season")
      } else {
            aux <- if (is.null(names(obj$Dates))) {
                  as.integer(substr(obj$Dates[[1]]$start, 6,7))      
            } else {
                  as.integer(substr(obj$Dates$start, 6,7))      
            }
            unique(aux)
      }
}
# End


#' @title Get years as a factor
#' @description Extract the year as a factor (e.g. for computing annual statistics)
#' @param obj Any object extending the station or grid classes
#' @return A vector of years of the same length as the time dimension of the object, 
#' seasonally-adjusted in the case of year-crossing seasons (e.g. DJF). See details.
#' @details The function performs a very basic operation, extracting the year element from the 
#' dates previously converted to POSIXlt. The trick lies in the year-crossing seasons. For instance:
#'  by convention, winter 2001 encompasses December 2000 and January, February 2001. Therefore, in order to compute
#' annual statistics for a year-crossing season, it is necessary to modify first the vector of years, 
#' and assign year 2001 to the preceding December. Similarly, the next December 2001 belongs to winter 2002,
#'  and so on... The function is useful for computing and/or plotting annual statistics, seasonal climatologies ... 
#' @note Warning:
#' The function should no be used to extract the vector of actual date years
#' @family get.helpers
#' @author J. Bedia 
#' @export
#' @examples \donttest{ 
#' require(climate4R.datasets)
#' data("NCEP_Iberia_hus850")
#' getSeason(NCEP_Iberia_hus850)
#' # Winter 1983-2002
#' range(NCEP_Iberia_hus850$Dates$start)
#' ## Time series for the first point
#' # Dates vector
#' time <- as.POSIXlt(NCEP_Iberia_hus850$Dates$start, tz = "GMT")
#' hus850 <- NCEP_Iberia_hus850$Data[ ,1,1]
#' plot(time, hus850, ty = "l")
#' ## Computation of the annual series for winter specific humidity:
#' par(mfrow = c(2,1))
#' ## Wrong:
#' years <- as.POSIXlt(NCEP_Iberia_hus850$Dates$start)$year + 1900
#' x <- tapply(hus850, INDEX = list(years), FUN = mean)
#' plot(unique(years), x, ty = "b")
#' points(1990, x[1], col = "red", cex = 2, lwd = 2)
#' ## Correct:
#' years <- getYearsAsINDEX(NCEP_Iberia_hus850)
#' x <- tapply(hus850, INDEX = years, FUN = mean)
#' plot(unique(years), x, ty = "b")
#' par(mfrow = c(1,1))
#' 
#' #Station data:
#' data("VALUE_Iberia_tas")
#' getYearsAsINDEX(VALUE_Iberia_tas)
#' }


getYearsAsINDEX <- function(obj) {
      season <- getSeason(obj)
      aux.dates <- getRefDates(obj)
      yrs <- as.numeric(substr(aux.dates,1,4))
      mon <- as.numeric(substr(aux.dates,6,7))
      if (identical(yrs, unique(yrs))) {
            yrs
            if (!identical(season, sort(season))) {
                  yrs <- yrs + 1
            }
      } else {    
            if (!identical(season, sort(season))) {
                  yy <- unique(yrs)[-1]
                  aux <- match(mon, season)
                  brks <- c(1, which(diff(aux) < 0) + 1, length(aux) + 1)
                  l <- lapply(1:(length(brks) - 1), function(x) {
                        a <- yrs[brks[x]:(brks[x + 1] - 1)]
                        return(rep(yy[x], length(a)))
                  })
                  yrs  <- do.call("c", l)
            }
      }
      return(yrs)
}
# End

#' @title  Retrieve dimensions attribute
#' @description Retrieve dimensions attribute
#' @param obj A grid or station object
#' @return A character vector with the dimensions attribute of the object's \code{Data} component.
#' @keywords internal
#' @family get.helpers
#' @export
#' @author J. Bedia
#' @examples \donttest{ 
#' require(climate4R.datasets)
#' #Station data:
#' data("VALUE_Iberia_tas")
#' getDim(VALUE_Iberia_tas)
#' #Regular grid
#' data("NCEP_Iberia_hus850")
#' getDim(NCEP_Iberia_hus850)
#' }

getDim <- function(obj) {
      attr(obj[["Data"]], "dimensions")
}

#' @title  Retrieve array shape 
#' @description Retrieve array attributes 'dimensions' and 'dim'
#' @param obj A grid or station object
#' @return An integer vector with dim values, labelled with the \code{"dimension"} attribute names
#' @family get.helpers
#' @keywords internal
#' @export
#' @author J. Bedia
#' @examples \donttest{ 
#' require(climate4R.datasets)
#' #Station data:
#' data("VALUE_Iberia_tas")
#' getShape(VALUE_Iberia_tas)
#' #Regular grid
#' data("NCEP_Iberia_hus850")
#' getShape(NCEP_Iberia_hus850)
#' }

getShape <- function(obj, dimension = NULL) {
      dimNames <- getDim(obj)
      shape <- dim(obj[["Data"]])
      if (!is.null(dimension)) {
            ind <- match(dimension, dimNames)
            if (anyNA(ind)) message("NOTE: Some of the specified dimensions do not exist in the grid")
            shape <- shape[ind]
            dimNames <- dimNames[ind]
      }
      names(shape) <- dimNames
      return(shape)
}

#' @title Leap years
#' @description Return the indices of leap years in a vector of years
#' @param years a integer vector of (gregorian) years
#' @return a vector of indices of the position of leap years
#' @references \url{https://en.wikipedia.org/wiki/Leap_year#Algorithm}
#' @keywords internal
#' @export
#' @family get.helpers
#' @author J. Bedia 
#' @examples
#' leap.years <- which.leap(1885:1937)
#' (1885:1937)[leap.years]

which.leap <- function(years) {
      which((years %% 4 == 0) & ((years %% 100 != 0) | years %% 400 == 0))
}


#' @title Retrieve reference dates
#' @description Retrieves the \code{$start} component of the \code{$Dates} element of a grid or station.
#'  In case of multigrids, the dates from the first grid are returned (i.e.: \code{grid$Dates[[1]]$start}).
#' @param obj A grid or station object
#' @param which Character string. The element of the \code{$Dates} list to be extracted.
#'  Default to \code{"start"}, otherwise \code{"end"}.
#' @return A character vector of dates
#' @details This utility function should be used always that reference dates are internally required by any function.
#' It takes into account the particular case when \code{\link{redim}} is applied to create a singleton \code{"var"}
#'  dimension: in this situation, the \code{"Dates"} component is not changed accordingly, and thus the criterion to
#'  distinguish multigrids from redim-ed grids is the shape of the array.
#' @keywords internal
#' @family get.helpers
#' @export
#' @author J. Bedia
#' @examples \donttest{ 
#' require(climate4R.datasets)
#' #Station data:
#' data("VALUE_Iberia_tas")
#' str(getRefDates(VALUE_Iberia_tas))
#' #Regular grid
#' data("NCEP_Iberia_hus850")
#' str(getRefDates(NCEP_Iberia_hus850, which = "end"))
#' }

getRefDates <- function(obj, which = "start") {
    which <- match.arg(which, choices = c("start", "end"))
    ind <- switch(which, "start" = 1, "end" = 2)
    if (("var" %in% getDim(obj)) && (getShape(obj, "var") > 1)) {
        obj$Dates[[1]][[ind]]
    } else {
        obj$Dates[[ind]]
    }
}


#' @title Select an appropriate *pply function
#' @description Selects an appropriate *pply function depending on the parallel choice
#' @param parallel.pars The output from \code{\link{parallelCheck}}, passing the parallelization options.
#' @param .pplyFUN *pply function to be used. Current options are \code{"apply"}, \code{"lapply"}, \code{"sapply"} and \code{"mapply"}.
#' @details The output function will be either \code{parallel::parLapply} or \code{lapply} for \code{.pplyFUN = "lapply"},
#' depending of whether parallelization is enabled or not. Same for \code{apply} and \code{parallel::parApply} and so on.
#' @section Warning:
#' 
#' From the \code{\link[parallel]{makeCluster}} help: 
#' \dQuote{It is good practice to shut down the workers by calling \code{stopCluster}}. It is therefore recommended the
#' following line of code after using this function: \code{on.exit(parallel::stopCluster(parallel.pars$cl))}
#' 
#' @importFrom parallel parApply parLapply
#' @return A function
#' @keywords internal
#' @family parallel.helpers
#' @export
#' @author J. Bedia

selectPar.pplyFun <- function(parallel.pars, .pplyFUN = c("apply", "lapply", "sapply", "mapply")) {
    .pplyFUN <- match.arg(.pplyFUN, choices = c("apply", "lapply", "sapply", "mapply"))
    if (parallel.pars$hasparallel) {
        if (.pplyFUN == "apply") {
            fun <- function(...) {
                parallel::parApply(cl = parallel.pars$cl, ...)
            }
        } else if (.pplyFUN == "lapply") {
            fun <- function(...) {
                parallel::parLapply(cl = parallel.pars$cl, ...)
            }
        } else if (.pplyFUN == "sapply") {
            fun <- function(...) {
                parallel::parSapply(cl = parallel.pars$cl, ...)
            }
        } else if (.pplyFUN == "mapply") {
            fun <- function(...) {
                parallel::clusterMap(cl = parallel.pars$cl, ..., SIMPLIFY = TRUE)
            }
        }
    } else {
        fun <- switch(.pplyFUN,
                      "apply" = apply,
                      "lapply" = lapply,
                      "sapply" = sapply,
                      "mapply" = mapply)
    }
    return(fun)
}


#' @title Spatial check of input grids
#' @description Checks that the input grids have the same lon/lat shape
#' @param ... Input grids to be compared
#' @param dimension Character vector. Dimensions to check for consistency. Several dimensions can be tested at once.
#' Possible values are: \code{"var"}, \code{"member"}, \code{"time"}, \code{"lat"} and \code{"lon"}.  
#' @return In case of spatial inconsistency of any of the inputs grids, 
#' the function stops the execution of the current expression, with an error message.
#' @description The function simply compares the size of the dimension across input grids for 
#' structural consistency of their respective data arrays. It does not make any additional checks on the metadata,
#' so equal array shapes do not entail full temporal/spatial consistency.
#' @seealso \code{\link{getShape}}
#' @keywords internal 
#' @export
#' @family check.helpers
#' @author J Bedia
#' @examples \donttest{ 
#' require(climate4R.datasets)
#' data("VALUE_Iberia_tas")
#' data("VALUE_Iberia_pr")
#' checkDim(VALUE_Iberia_tas, VALUE_Iberia_pr) # ok, go on
#' data("EOBS_Iberia_tas")
#' data("EOBS_Iberia_pr")
#' checkDim(EOBS_Iberia_tas, EOBS_Iberia_pr) # ok, go on
#' data("NCEP_Iberia_psl")
#' try(checkDim(EOBS_Iberia_tas, EOBS_Iberia_pr, NCEP_Iberia_psl)) # inconsistent dimensions
#' data("CFS_Iberia_tas")
#' try(checkDim(NCEP_Iberia_psl, CFS_Iberia_tas, dimensions = "member")) 
#' checkDim(VALUE_Iberia_tas, VALUE_Iberia_pr) # ok, go on
#' }


checkDim <- function(..., dimensions = c("var", "member", "time", "lat", "lon")) {
    grid.list <- list(...)
    grid.list <- lapply(grid.list, "redim")
    dimensions <- match.arg(dimensions, choices = c("var", "member", "time", "lat", "lon"), several.ok = TRUE)
    dimlist <- lapply(dimensions, function(x) vapply(grid.list, "getShape", integer(1), x))
    oops <- dimensions[which(!sapply(dimlist, function(x) all(x == x[1])))]
    if (length(oops) > 0) {
        stop("Inconsistent sizes found for dimensions: ", paste(oops, collapse = ", "), call. = FALSE)
    }
}


#' @title Seasonal check
#' @description Check the consistency of seasons across the input grids
#' @param ... Input grids to be compared
#' @return In case of seasonal inconsistency of any of the inputs grids, the function stops the execution 
#' of the current expression, with an error message.
#' @seealso \code{\link{getSeason}}
#' @author J Bedia
#' @keywords internal
#' @family check.helpers
#' @export
#' @examples \donttest{ 
#' require(climate4R.datasets)
#' data("EOBS_Iberia_tas")
#' data("EOBS_Iberia_pr")
#' grid1 <- subsetGrid(EOBS_Iberia_tas, season = 1)
#' grid2 <- subsetGrid(EOBS_Iberia_pr, season = 1)
#' checkSeason(grid1, grid2) # ok, go on
#' try(checkSeason(grid1, grid2, EOBS_Iberia_tas)) # oops
#' # However note that it is insensitive to the level of temporal aggregation:
#' agg <- aggregateGrid(EOBS_Iberia_tas, aggr.m = list(FUN = "mean", na.rm = TRUE))
#' identical(getSeason(agg), getSeason(EOBS_Iberia_tas))
#' checkSeason(agg, EOBS_Iberia_tas)
#' 
#' #Station data
#' data("VALUE_Iberia_tas")
#' data("VALUE_Iberia_pr")
#' checkSeason(VALUE_Iberia_tas, VALUE_Iberia_pr)
#' }


checkSeason <- function(...) {
    grid.list <- list(...)
    sealist <- lapply(grid.list, "getSeason")
    oops <- sapply(sealist, function(x) all(x == sealist[[1]]))
    if (!all(oops)) stop("Inconsistent seasons among input grids")
}    


#' @title Variable name check
#' @description Check the consistency of variable (short)names across grids/multigrids
#' @param ... Input grids to be compared
#' @param check.order Check if the order is consistent across input grids. This is the default. When set to
#' \code{TRUE}, the function does not only check that the variable names match, but also that their order is the same
#' (in case of multigrids)
#' @return In case of inconsistency of any of the inputs grids, the function stops the execution 
#' of the current expression, with an error message.
#' @seealso \code{\link{getVarNames}}
#' @author J Bedia
#' @keywords internal
#' @family check.helpers
#' @export
#' @examples \donttest{ 
#' require(climate4R.datasets)
#' data("NCEP_Iberia_hus850", "NCEP_Iberia_psl", "NCEP_Iberia_ta850")
#' #' # Testing grids that contain different variables:
#' try(checkVarNames(NCEP_Iberia_hus850, NCEP_Iberia_psl))
#' # Testing multigrid contents:
#' mg <- makeMultiGrid(NCEP_Iberia_hus850, NCEP_Iberia_psl, NCEP_Iberia_ta850)
#' mg1 <- makeMultiGrid(NCEP_Iberia_hus850, NCEP_Iberia_ta850, NCEP_Iberia_psl)
#' mg2 <- makeMultiGrid(NCEP_Iberia_ta850, NCEP_Iberia_psl, NCEP_Iberia_hus850)
#' # The input grids contains the same variables, but their ordering is different:
#' checkVarNames(mg, mg1, mg2, check.order = FALSE) # OK
#' # If we do check.order = TRUE (the default), it gives an error:
#' try(checkVarNames(mg, mg1, mg2, check.order = TRUE)) # Error
#' }

checkVarNames <- function(..., check.order = TRUE) {
    grid.list <- list(...)
    if (length(grid.list) == 1) {
        message("NOTE: Only one input grid: 'checkVarNames' was skipped")   
    } else {
        out <- lapply(grid.list, FUN = "getVarNames")
        lengths <- sapply(out, "length")
        if (!all(lengths == lengths[[1]])) {
            stop("The number of variables stored in the input grids differ")
        }
        if (!all(sapply(2:length(out), function(x) all(out[[x]] %in% out[[1]])))) {
            stop("Variable names of the input grids differ")
        }
        if (check.order) {
            if (!all(sapply(2:length(out), function(x) identical(out[[x]], out[[1]])))) {
                stop("Variable ordering of the input grids differ")
            }
        }
    }
}

#' @title Temporal consistency check across input grids
#' @description Check the consistency of temporal aspects across input climate4R objects
#' @param ... Input grids to be compared. Either as a list, or comma sepparated
#' @return In case of inconsistency of any of the inputs grids, the function stops the execution 
#' of the current expression, with an error message.
#' @seealso \code{\link{checkSeason}}, for a more specific -less strict- type of check
#' @details The function takes into account the temporal resolution of the input grids via 
#' the \code{\link{getTimeResolution}} helper, and ensures its consistency. Then, it makes the 
#' comparison among inputs at that temporal level of resolution. For instance, if the input grids are all 
#' monthly, the function will intercompare only the consistency of years AND months, but will disregard
#' lower temporal units (days, hours). 
#' @author J Bedia
#' @keywords internal
#' @importFrom magrittr %>% %<>% 
#' @family check.helpers
#' @export

checkTemporalConsistency <- function(...) {
    grid.list <- list(...)
    if (length(grid.list) == 1) grid.list %<>% unlist(recursive = FALSE)
    if (length(grid.list) < 2) message("NOTE: Only one grid passed as input. Nothing was done")
    timeres <- sapply(grid.list, "getTimeResolution") %>% unique() 
    if (length(timeres) > 1) stop("Binding grids of different temporal resolutions is not allowed\nDid you mean \'makeMultiGrid\'?")
    refdates <- getRefDates(grid.list[[1]]) %>% as.character()
    refmon <- substr(refdates, 6, 7)
    refyr <- substr(refdates, 1, 4)
    refday <- substr(refdates, 9, 10)
    mssg <- "Input data are not temporally consistent"
    if (timeres == "MM") {
        for (i in 2:length(grid.list)) {
            testmon <- getRefDates(grid.list[[i]]) %>% as.character() %>% substr(start = 6, stop = 7)
            testyr <- getRefDates(grid.list[[i]]) %>% as.character() %>% substr(start = 1, stop = 4)
            if (!identical(refmon, testmon) | !identical(refyr, testyr)) stop(mssg)
        }
    } else if (timeres == "YY") {
        for (i in 2:length(grid.list)) {
            testyr <- getRefDates(grid.list[[i]]) %>% as.character() %>% substr(start = 1, stop = 4)
            if (!identical(refyr, testyr)) stop(mssg)
        }
    } else if (timeres == "DD") {
        for (i in 2:length(grid.list)) {
            testmon <- getRefDates(grid.list[[i]]) %>% as.character() %>% substr(start = 6, stop = 7)
            testyr <- getRefDates(grid.list[[i]]) %>% as.character() %>% substr(start = 1, stop = 4)
            testday <- getRefDates(grid.list[[i]]) %>% as.character() %>% substr(start = 9, stop = 10)
            if (!identical(refday, testday) | !identical(refyr, testyr) | !identical(refmon, testmon)) stop(mssg)
        }
    } else {
        for (i in 2:length(grid.list)) {
            if (!identical(as.POSIXlt(refdates), as.POSIXlt(getRefDates(grid.list[[i]])))) stop(mssg)
        }
    }   
}



#' @title Check if the input grid is regular
#' @description Check if the coordinates in the data are regular or irregular
#' @param grid a grid or station data, or the object returned by function \code{\link{getGrid}}
#' @return Logical.
#' @seealso \code{\link{getGrid}}
#' @author M Iturbide
#' @keywords internal
#' @family is.helpers
#' @seealso typeofGrid, for a character representation of the type of grid
#' @export
#' @examples \donttest{ 
#' require(climate4R.datasets)
#' data("EOBS_Iberia_tas")
#' isRegular(EOBS_Iberia_tas)
#' isRegular(getGrid(EOBS_Iberia_tas))
#' data("VALUE_Iberia_tas")
#' isRegular(VALUE_Iberia_tas)
#' }

isRegular <- function(grid) {
    gr <- tryCatch({getGrid(grid)}, error = function(err) {grid})
    x <- sort(gr$x)
    y <- sort(gr$y)
    if (!is.null(attr(gr, "resX")) && !is.null(attr(gr, "resY"))) {
        if (attr(gr, "resX") == 0 && attr(gr, "resY") == 0) {
            FALSE 
        } else {
            TRUE
        }
    } else {
        if (length(x) == 1 && length(y) == 1) {
            FALSE
        } else {
            xdists <- lapply(1:(length(x) - 1), function(l) {
                x[l + 1] - x[l]
            })
            ydists <- lapply(1:(length(y) - 1), function(l) {
                y[l + 1] - y[l]
            })
            xa <- sum(unlist(xdists) - unlist(xdists)[1])
            ya <- sum(unlist(ydists) - unlist(ydists)[1])
            if (any(abs(c(xa, ya)) > 1e-05)) {
                FALSE
            } else {
                TRUE
            }
        }
    }
}

#End

#' @title Get time resolution
#' @description Get the time resolution of the input grid as a character
#' @param grid a grid or station data
#' @return A character representation of the time resolution. Current possible values are:
#' \code{"1h", "3h", "6h", "12h", "DD", "MM", "YY"}. If none of these matches, \code{"unknown"} is returned 
#' @author J Bedia
#' @keywords internal
#' @family get.helpers dateutils
#' @export
#' @examples \donttest{ 
#' require(climate4R.datasets)
#' data("EOBS_Iberia_tas")
#' getTimeResolution(EOBS_Iberia_tas)
#' monthly.grid <- aggregateGrid(EOBS_Iberia_tas, aggr.m = list(FUN = "mean", na.rm = TRUE))
#' stopifnot(identical(getTimeResolution(monthly.grid), "MM"))
#' annual.grid <- aggregateGrid(monthly.grid, aggr.y = list(FUN = "mean", na.rm = TRUE))
#' stopifnot(identical(getTimeResolution(annual.grid), "YY"))
#' }

getTimeResolution <- function(grid) {
    aux <- getRefDates(grid) 
    if (length(aux) == 1) {
        message("The input grid is a climatology: temporal resolution set to \'unknown\'")
        out <- "unknown"
    } else {
        dft <- difftime(aux[2], aux[1], units = "hours") %>% as.numeric()
        out <- if (dft == 1) {
            "1h"
        } else if (dft == 3) {
            "3h"    
        } else if (dft == 6) {
            "6h"
        } else if (dft == 12) {
            "12h"
        } else if (dft >= 23 & dft <= 24) {
            "DD"
        } else if (dft >= 672 & dft <= 745) {
            "MM"
        } else if (dft >= 8640 & dft <= 8784) {
            "YY"
        } else {
            "unknown"
        }
    }
    return(out)
}


#' @title Conversion 2D matrix into a 3D array for station data
#' @description Converts a 2D matrix of the form [time, lonlat] to a 3D array of the form
#' [time,lat,lon], in this order. Mainly for PCA analysis and grid reconstruction.
#' @param mat2D A 2D matrix with time in rows and lonlat in columns, as returned 
#' by \code{\link{array3Dto2Dmat.stations}} 
#' @param x unique X coordinates of the points, in ascending order
#' @param y As argument \code{x}, for the Y coordinates
#' @return A 3-dimensional array with the dimensions ordered: [time,lat,lon]
#' @importFrom abind abind
#' @family internal.helpers
#' @details The function is the inverse of \code{\link{array3Dto2Dmat.stations}} 
#' @author M. Iturbide
#' @seealso \code{\link{array3Dto2Dmat}}, which performs the inverse operation.

mat2Dto3Darray.stations <- function(mat2D, x, y) {
    mat <- matrix(NA, ncol = length(x), nrow = length(y))
    aux.list <- lapply(1:nrow(mat2D), function(i) {
        diag(mat) <- mat2D[i, ]
        mat
    })
    arr <- unname(do.call("abind", c(aux.list, along = -1)))
    aux.list <- NULL
    arr <- aperm(arr, perm = c(1,3,2))
    attr(arr, "dimensions") <- c("time", "lat", "lon")
    return(arr)
}
#End

#' @title Conversion of a 3D array to a 2D matrix for station data
#' @description Converts 3D arrays of the form [lon,lat,time] -not strictly in this order-,
#'  to 2D matrices of the form [time, grid-point], in this order. Mainly for PCA analysis.
#' @param array3D A 3-dimensional array with longitude, latitude and time dimensions
#' @return A 2-dimensional matrix with time in rows and grid-points in columns.
#' @author M. Iturbide
#' @keywords internal
#' @family internal.helpers
#' @seealso \code{\link{mat2Dto3Darray}}, which performs the inverse operation
#' 
array3Dto2Dmat.stations <- function(array3D) {
    dims <- dim(array3D)
    aux.list <- lapply(1:dims[1], function(i) {
        if (!is.matrix(array3D[i, ,])) {
            array3D[i, ,]
        } else{
            diag(array3D[i, ,])      
        }
    })
    mat <- unname(do.call("abind", c(aux.list, along = -1)))
    aux.list <- NULL
    attr(mat, "dimensions") <- c("time", "loc")
    return(mat)
}

#End


#' @title Get grid variable names
#' @description Get the names of the variables of a (multi)grid
#' @param grid A grid or station data
#' @param type Character. Should either the \code{"short"} (default) or the \code{"long"} variable name(s) be returned?.
#' See Details.
#' @param append.level Append the vertical level to the short name string?. Default to \code{TRUE}, ignored when \code{type = "long"}).
#'  See details.
#' @details 
#' 
#' \strong{About short names}
#' 
#' Note that for variables with vertical levels, if the option \code{append.level} is set to \code{TRUE} (the default),
#'  the function will "construct" the shortname following the climate4R convention,
#' that defines the code of the variable plus its vertical level, separated by a \dQuote{@@} symbol. 
#' This way, when querying grid variable names, the behaviour of \code{getVarNames} is different than directly accessing
#' the \code{varName} element of the grid structure as in, e.g.: \code{grid$Variable$varName}. Using \code{getVarNames} 
#' with \code{append.level = TRUE} is recommended in most applications when the variable string is needed.
#' 
#' \strong{About long names}
#' 
#' The long name is an optional attribute of a grid (not so the short name), and it is sometimes undefined (particularly in station datasets).
#' In this case, the function returns \code{NULL}.
#' @return A character vector with the variable name(s) in the indicated \code{type} format.
#' @author J Bedia
#' @keywords internal
#' @family get.helpers
#' @export
#' @examples \donttest{ 
#' require(climate4R.datasets)
#' data("CFS_Iberia_pr")
#' getVarNames(CFS_Iberia_pr)
#' grid <- CFS_Iberia_pr
#' getVarNames(CFS_Iberia_pr, "long")
#' ## Example with a multigrid
#' data(NCEP_Iberia_hus850, NCEP_Iberia_psl, NCEP_Iberia_ta850)
#' mg <- makeMultiGrid(NCEP_Iberia_hus850, NCEP_Iberia_psl, NCEP_Iberia_ta850)
#' getVarNames(mg)
#' # The level can be removed id needed:
#' getVarNames(mg, append.level = FALSE)
#' data("CFS_Iberia_tas", "CFS_Iberia_pr")
#' mg <- makeMultiGrid(CFS_Iberia_tas, CFS_Iberia_pr)
#' getVarNames(mg)
#' getVarNames(mg, "long")
#' ## Example with station data
#' data("VALUE_Iberia_tas")
#' ## The long name is an optional attribute, and might be undefined:
#' getVarNames(VALUE_Iberia_tas, "long")
#' }

getVarNames <- function(grid, type = c("short", "long"), append.level = TRUE) {
    stopifnot(is.logical(append.level))
    type <- match.arg(type, c("short", "long"))
    if (type == "short") {
        vnames <- grid$Variable$varName
        if (append.level) {
            if (!any(grepl(pattern = "@", vnames))) {
                lev <- grid$Variable$level
                ind <- suppressWarnings(which(!is.na(lev) & !is.null(lev)))
                if (length(ind) > 0) {
                    vnames[ind] <- vapply(ind, FUN.VALUE = character(1L), FUN = function(x) paste(vnames[x], lev[x], sep = "@"))    
                }
            }
        } else {
            vnames <- gsub("@.*$", "", vnames)    
        }
    } else {
        vnames <- gsub("\\\"","", attr(grid$Variable, "longname"))
    }
    return(vnames)
}

#' @title Get grid vertical levels
#' @description A helper function that returns a vector of the variable(s) vertical levels
#' @param grid Input grid
#' @param var.index A vector of indices indicating the positions of the variables whose level will be returned.
#' This can be either an integer vector of positions, or a character vector with the short names of the variables contained
#'   in the multigrid (see the \code{\link{getVarNames}} helper for more details). Default to \code{NULL} (levels of all variables)
#' @return A (short-)named vector with the vertical levels of the (multi)grid variable(s)
#' @author J Bedia
#' @keywords internal
#' @export
#' @family get.helpers
#' @examples \donttest{ 
#' require(climate4R.datasets)
#' data("CFS_Iberia_hus850")
#' getGridVerticalLevels(CFS_Iberia_hus850)
#' # Surface variables usually have an undefined vertical level <NA>
#' data("EOBS_Iberia_tas")
#' getGridVerticalLevels(EOBS_Iberia_tas)
#' 
#' data("NCEP_Iberia_hus850")
#' data("NCEP_Iberia_psl")
#' data("NCEP_Iberia_ta850")
#' mg <- makeMultiGrid(NCEP_Iberia_hus850, NCEP_Iberia_psl, NCEP_Iberia_ta850)
#' getGridVerticalLevels(mg)
#' # Use of var.index to select some variables:
#' # either by their shortname...
#' getGridVerticalLevels(mg, var.index = c("hus@@850","ta@@850"))
#' # ... or by index position in the multigrid
#' getGridVerticalLevels(mg, var.index = c(1,3))
#' }

getGridVerticalLevels <- function(grid, var.index = NULL) {
    var.names <- getVarNames(grid)
    var.index <- if (is.null(var.index)) {
        1:length(var.names)
    } else if (is.character(var.index)) {
        match(var.index, var.names)
    } else {
        var.index
    }
    if (!all(var.index %in% 1:length(var.names))) {
        stop("Invalid 'var.index' definition", call. = FALSE)
    }
    levs <- grid$Variable$level[var.index]
    # NULL values (surface vars) are converted to NA
    if (any(is.null(levs))) {
        levs[which(is.null(levs))] <- NA
    }
    names(levs) <- var.names[var.index]
    return(levs)
}


#' @title Type of grid
#' @description Describe the type of grid according to the spatial representation of the data
#' @param grid An input grid
#' @return A character vector decribing the type of spatial reference. Three disjoint values possible: \code{"station"}, \code{"rotated_grid"} and \code{"regular_grid"}
#' @keywords internal
#' @author J Bedia
#' @family get.helpers
#' @export

typeofGrid <- function(grid) {
    ref.grid <- getGrid(grid)
    if (attr(ref.grid, "resX") == 0) {
        "station"
    } else if (is.na(attr(ref.grid, "resX"))) {
        "rotated_grid"
    } else {
        "regular_grid"
    }
}    





#' @title Level depth in a list
#' @description Level depth in a list 
#' @param this list
#' @export
#' @return number of nesting lists
#' @keywords internal
#' @author M. Iturbide 

gridDepth <- function(this) {
    that <- this
    i <- 0
    while (is.list(that)) {
        i <- i + 1
        that <- that[[1]]
    }
    return(i)
}
#end



#' Reorder loc dimension in irregular grids
#' 
#' Retrieves a grid that is with the 'loc' dimension reordered.
#'  Multimember multigrids are supported. Subroutine of \code{\link{subsetGrid}}.
#'
#' @param grid Input irregular grid to be subset (possibly a multimember/multigrid).
#' @param lon Logical to order data according to longitudes.
#' @param lat Logical to order data according to latitudes.
#' @return An irregular grid (or multigrid).
#' @keywords internal
#' @importFrom magrittr %>% %<>%
#' @export
#' @author M. Iturbide

reorderStation <- function(grid, axis = c("x", "y")) {
    dimNames <- getDim(grid)
    axis <- match.arg(axis, choices = c("x", "y"))
    if (isRegular(grid)) stop("This function is applied only to irregular grids")
    indloc <- order(getCoordinates(grid)[[axis]])
    grid$Data <- asub(grid$Data, idx = indloc, dims = grep("loc", dimNames), drop = FALSE)
    attr(grid$Data, "dimensions") <- dimNames
    grid$xyCoords <- grid$xyCoords[indloc,]
    if ("Metadata" %in% names(grid)) {
        if ("station_id" %in% names(grid$Metadata)) grid$Metadata$station_id <- grid$Metadata$station_id[indloc]
        if ("name" %in% names(grid$Metadata)) grid$Metadata$name <- grid$Metadata$name[indloc]
        if ("altitude" %in% names(grid$Metadata)) grid$Metadata$altitude <- grid$Metadata$altitude[indloc]      
        if ("source" %in% names(grid$Metadata)) grid$Metadata$source <- grid$Metadata$source[indloc]
    }
    return(grid)
}
#end


#' @title Get grid coordinates as 2D matrix
#' @description Obtain grid coordinates as 2D matrix
#' @param grid An input grid
#' @return A 2D matrix of x-y coordinates (in this order)
#' @keywords internal
#' @author J. Bedia
#' @family get.helpers
#' @export

get2DmatCoordinates <- function(grid) {
    if (typeofGrid(grid) == "regular_grid") {
        coords <- getCoordinates(grid) 
        coords$lon <- NULL
        coords$lat <- NULL
        aux <- expand.grid(coords)
        aux[order(aux[,1]), ]
    } else if (typeofGrid(grid) == "station") {
        getCoordinates(grid)
    } else if (typeofGrid(grid) == "rotated_grid") {
        cbind(as.vector(getCoordinates(grid)$lon), as.vector(getCoordinates(grid)$lat))
    }
}


#' Check if object is a grid
#' 
#'
#' @param grid Input object.
#' @return Logical.
#' @keywords internal
#' @export
#' @author M. Iturbide

isGrid <- function(grid) {
    if (is.list(grid)) {
        all(c("Variable", "Data", "xyCoords", "Dates") %in% names(grid))
    } else {
        FALSE
    }
}


#' Check if object is a multigrid
#' 
#'
#' @param grid Input object.
#' @return Logical.
#' @keywords internal
#' @export
#' @author M. Iturbide

isMultigrid <- function(grid) {
    if (is.list(grid)) {
        if (all(c("Variable", "Data", "xyCoords", "Dates") %in% names(grid))) {
            gridDepth(grid$Dates) > 1
        } else {
            FALSE
        }
    } else {
        FALSE
    }
}



#' @title Get grid units
#' @description Get the \code{"units"} attribute of a grid
#' @param grid An input grid
#' @param var Character vector of variable length. Variable short name(s) whose units are returned.
#'  Only makes sense in the case of multigrids storing several variables.
#'   Otherwise ignored.
#' @return Returns the \code{"units"} attribute 
#' @keywords internal
#' @export
#' @author J Bedia
#' @family get.helpers unit.helpers
#' @examples \donttest{ 
#' require(climate4R.datasets)
#' data(NCEP_Iberia_ta850)
#' getGridUnits(NCEP_Iberia_ta850)
#' data(NCEP_Iberia_hus850)
#' getGridUnits(NCEP_Iberia_hus850)
#' data(NCEP_Iberia_psl)
#' getGridUnits(NCEP_Iberia_psl)
#' mf <- makeMultiGrid(NCEP_Iberia_hus850, NCEP_Iberia_psl, NCEP_Iberia_ta850)
#' getGridUnits(mf)
#' getVarNames(mf)
#' getGridUnits(mf, "hus@850")
#' getGridUnits(mf, var = c("hus@850", "ta@850"))
#' }

getGridUnits <- function(grid, var = NULL) {
    uds <- attr(grid$Variable, "units") %>% gsub(pattern = "\\\"", replacement = "")
    if (isMultigrid(grid)) {
        if (is.null(var)) {
            message("NOTE: The input is a multigrid: Units of all variables are shown.\nUse argument 'var' for displaying the units of a particular variable")
        } else {
            vn <- getVarNames(grid)
            if (!all(var %in% vn)) stop("At least one variable in \'var\' was not found. Use \'getVarNames\' for help")
            uds <- uds[match(var, vn)]
        }
    } else {
        if (!is.null(var)) warning("The input \'grid\' is not a multigrid: argument \'var\' was ignored")
    }
    return(uds)
}


#' @title Set grid units
#' @description Set the \code{"units"} attribute of a grid
#' @param grid An input grid
#' @param unit.string Character string: a udunits-parseable character string vector.
#'  See details.
#' @param var In case of multigrids, the names of the variables whose units attribute is
#'   to be updated (see examples).
#' @return Retunrs (invisible) the same input grid with the new \code{"units"} 
#' attribute in \code{"$Variable"} list element.
#' @details 
#' The length of the \code{unit.string} vector should match the number of variables
#'  within the grid (in case of \code{multiGrids}), i.e., that of 
#'  getVarNames(grid) or the length of \code{var}, in case the latter is used.
#' @export
#' @author J Bedia
#' @family get.helpers unit.helpers
#' @examples \donttest{ 
#' require(climate4R.datasets)
#' data(NCEP_Iberia_hus850, NCEP_Iberia_psl, NCEP_Iberia_ta850)
#' getGridUnits(NCEP_Iberia_hus850)
#' getGridUnits(NCEP_Iberia_psl)
#' mf <- makeMultiGrid(NCEP_Iberia_hus850, NCEP_Iberia_psl, NCEP_Iberia_ta850)
#' getGridUnits(mf)
#' mf2 <- setGridUnits(mf, unit.string = c("1", "Pa", "Kelvin"))
#' getGridUnits(mf2)
#' # Arbitrary subsets of variables within the multigrid can be updated:
#' getVarNames(mf)
#' mf3 <- setGridUnits(mf, unit.string = c("1", "Pa"), var = c("hus@850", "psl"))
#' getGridUnits(mf3)
#' }

setGridUnits <- function(grid, unit.string, var = NULL) {
    stopifnot(isGrid(grid))
    vn <- getVarNames(grid)
    if (isMultigrid(grid)) {
        if (is.null(var)) {
            if (length(vn) != length(unit.string)) {
                stop("The length of the \'unit.string\' vector does not match the number of variables in the grid")
            }
            ind <- 1:length(vn)
        } else {
            if (!all(var %in% vn)) stop("At least one variable in \'var\' was not found. Use \'getVarNames\' for help")
            if (length(unit.string) != length(var)) stop("Inconsistent \'unit.string\' and \'var\' vector lengths")
            ind <- match(var, vn)
        }
    } else {
        if (length(unit.string) > 1) stop("\'unit.string\' vector should have length 1")
        ind <- 1L
    }
    attr(grid$Variable, "units")[ind] <- unit.string
    invisible(grid)
}


#' @title Set grid projection
#' @description Set the \code{"projection"} attribute of a climate4R grid
#' @param grid An input grid
#' @param proj Character string. A string specifying the projection attribute value. See details.
#' @return Returns (invisible) the same input grid with the new \code{"projection"} 
#' attribute in \code{"$xyCoords"} list element.
#' @details 
#' The projection string format is free, although usually a valid epsg code or proj4 string is introduced here.
#' @export
#' @author J Bedia
#' @family get.helpers 
#' @examples \donttest{ 
#' require(climate4R.datasets)
#' data(NCEP_Iberia_hus850)
#' getGridProj(NCEP_Iberia_hus850)
#' # We introduce the standard definition of WGS84 lat-lon projection
#' NCEP_Iberia_hus850 <- setGridProj(NCEP_Iberia_hus850,
#'                                   proj= "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs")
#' getGridProj(NCEP_Iberia_hus850)                 
#' }                  

setGridProj <- function(grid, proj = NULL) {
  stopifnot(isGrid(grid))
  attr(grid[["xyCoords"]], "projection") <- proj
  invisible(grid)
}

#' @title Obtain station names and IDs
#' @description Obtain station names and IDs
#' @param grid An input station data grid
#' @return A 2D matrix of name-ID (in this order)
#' @author M. Iturbide
#' @family get.helpers
#' @export

getStationID <- function(grid) {
      if (typeofGrid(grid) != "station") stop("This is not station data")
      if (!is.null(grid$Metadata$name) & !is.null(grid$Metadata$station_id)) {
            cbind(grid$Metadata$name, grid$Metadata$station_id)      
      } else {
            stop("Incomplete metadata")
      }
}

#' @title Check the structure of the grid
#' @description Check the structure of the grid
#' @param ... Input grids (including station data)
#' @return Diagnostic message
#' @author M. Iturbide
#' @family check.helpers
#' @export


checkGrid <- function(...) {
      grid.list <- list(...)
      if (is.null(names(grid.list))) {
            nmes <- as.character(as.list(substitute(list(...)))[-1L])
            if (length(nmes) < length(grid.list)) nmes <- paste0(nmes, 1:length(grid.list))
            names(grid.list) <- nmes
      }
      end.dates.order <- lapply(grid.list, function(x) 
            identical(1:length(x$Dates$end), order(x$Dates$end))
      )
      start.dates.order <- lapply(grid.list, function(x) 
            identical(1:length(x$Dates$start), order(x$Dates$start))
      )
      start.end.dates.length <- lapply(grid.list, function(x) 
            length(x$Dates$end) == length(x$Dates$start) & length(x$Dates$end) == getShape(redim(x))[["time"]]
      )
      member.info <- lapply(grid.list, function(x) 
            if (getShape(redim(x))[["member"]] > 1) {
                  if ("Members" %in% names(x)) {
                        mem <- length(x[["Members"]]) == getShape(redim(x))[["member"]]
                        if (!mem) warning("Length of member dimension in '$Data' and length of '$Members' do not match.") 
                  } else {
                        mem <- FALSE
                        warning("There is no '$Members' in grid.") 
                  }
                  if ("InitializationDates" %in% names(x) & mem) {
                        length(x[["InitializationDates"]]) == getShape(redim(x))[["member"]]
                  } else if (mem) {
                        warning("There is no '$InitializationDates' in grid.")
                        "WARNING"
                  } else {
                        FALSE
                  }
            } else {
                  TRUE
            }
      )
      message.df <- rbind(start.dates.order, start.dates.order, start.end.dates.length, member.info)
      if (any(message.df == FALSE, na.rm = TRUE)) warning("Errors found in grid structure.")
      if (any(message.df == "WARNING", na.rm = TRUE)) warning("WARNING: There is no element 'InitializationDates in grid'.")
      message.df[which(message.df == TRUE)] <- "OK"
      message.df[which(message.df == FALSE)] <- "ERROR"
      message.df
}


#' @title Obtain projection information of a grid
#' @description Obtain projection information of a grid
#' @param grid An input climate4R grid
#' @return A character string with the projection definition
#' @details The projection may be defined loosely as e.g., \dQuote{LatLon Projection},
#' or more formally as a proj4 string. This function just retrieves the projection attribute
#' as included within the target grid.
#' Grid projection definition (and/or re-projection operations) can be achieved by the function \code{projectGrid}
#' from the climate4R package \pkg{geoprocessor} (\url{https://github.com/SantanderMetGroup/geoprocessoR})
#' @author J. Bedia
#' @family get.helpers
#' @export

getGridProj <- function(grid) {
  if ("projection" %in% names(attributes(grid$xyCoords))) {
    return(attr(grid$xyCoords, "projection"))
  } else {
    warning("Undefined projection.\nConsider using projectGrid from package geoprocessoR for setting a valid proj definition\n<https://github.com/SantanderMetGroup/geoprocessoR>")
  }
}


#' @title Obtain Weather Typing (WT) classification
#' @description Obtain the classification attribute from \code{\link{clusterGrid}} outputs
#' @param grid An input climate4R grid. 
#' @return The attribute vector \code{"wt.index"}, when present. Otherwise \code{NULL}. See \code{\link{clusterGrid}} for details.
#' @author J. Bedia
#' @family get.helpers
#' @export

getWT <- function(grid) attr(x = grid, which = "wt.index")