R/spatial.R

In martindurocher/floodRFA: A R-package for regional frequency of floods

##################################################################
#' Corelation model
#'
#' Returns a vector or matrix of the correlation associated to
#' a given spatial correlation model.
#'
#' @param p Vector of row matrix of the model parameters. If it is
#'   a matrix, each row is set parameter.
#'
#' @param h Vector or matrix of the distance.
#'
#' @param type Correlation model
#'
#' @details
#'
#' The correlation model are Strings. The available model are:
#' exponential (exp), gaussian (gau), rational quadratic (rqd),
#' spherical (sph), circular (cir), power exponential (pow) and
#' matern (mat).
#'
#' @examples
#'
#' coord <- replicate(2,runif(30))
#' h <- as.matrix(dist(coord))
#'
#' Sexp <- corModel(c(.6,0), h,'exp')
#' Srq <- corModel(c(.3,.2), h,'rqd')
#' Smat <- corModel(c(.6,.1, 50), h,'mat')
#'
#' plot(as.vector(h),as.vector(Sexp))
#' points(as.vector(h),as.vector(Srq),col = 2)
#' points(as.vector(h),as.vector(Smat),col = 4)
#'
#' @export
##################################################################

corModel <- function(p, ...) UseMethod('corModel', p)

#' @export
#' @rdname corModel
corModel.numeric <- function(p, h, type = 'exp'){


  ## Scale the distance by the range
  if(p[1]<=0){
    stop('Invalid range parameter')
  } else {
    h <- h/p[1]
  }

  ## Calculate the model correlation
  if(type == 'exp'){
    ans <- exp(-3 * h)

  } else if(type =='gau'){
    ans <- exp(-3 * h * h)

  } else if(type == 'pow'){

    if(p[3] <= 0 | p[3]>2){
      stop('Invalid smooth parameter')
    }

    ans <- exp(-3 * h^p[3])

  } else if(type == 'rqd'){
    ans <-  1 / (1 + 19 * h^2 )

  } else if(type =='lin'){
    ans <- 1-pmax(0,pmin(1,h))

  } else if(type =='sph'){
    h[h>1] <- 1;
    ans <- 1 - 1.5 * h + .5 * h^3

  } else if(type =='cir'){
    h[h>1] <- 1
    ans <- 2/pi*(acos(h)- h * sqrt( 1 - h * h))

  } else if(type =='mat'){
    if(p[3]<= 0){
      stop('Invalid smooth parameter')
    }

    nu <- sqrt(2*p[3]) *h * 3
    ans <- 2^(1-p[3]) / gamma(p[3]) * nu^(p[3]) * besselK(nu,p[3])

  } else{
    stop('Invalid correlation model')
  }

  ## Add a nugget effect
  if(p[2] >1 | p[2] < 0){
    stop('Invalid nugget parameter')
  } else {
    ans <- (1-p[2]) * ans
    ans[h==0] <- 1
  }

  return(ans)
}

#' @export
corModel.matrix <- function(p, h, type = 'exp'){
  return(apply(p,1,corModel, h = h, type = type))
}

#' @export
corModel.data.frame <- function(p, h, type = 'exp'){
  return(corModel(as.matrix(p), h = h, type = type))
}

################################################################
#' Fitting a correlation model by least-squares
#'
#' Returns the parameter of the correlation model
#'
#' @param form Formula indication the correlation coefficient
#'     and the distance. For instance \code{x~h}
#'
#' @param h Vector of distance
#'
#' @param x Vector of correlation coefficients or Nx2 matrix
#'    of the distance and correlation coefficients.
#'
#' @param w Vector of weights.
#'
#' @param type Correlation model. See \code{\link{corModel}}
#'
#' @param nugget The nugget parameter. If \code{NULL} the
#'   parameter is optimized.
#'
#' @param smooth The smooth parameter for three parameters model.
#'
#' @param hPred A vector of distance for prediction (optional)
#'
#' @param startp The starting point for optimization. May be used
#'   when the nugget effect is estimated.
#'
#' @param ... Other parameters pass to the funtion optim.
#'
#' @export
#' @examples
#'
#' tt = 1:500
#' rr = corModel(c(200,.1), tt, type = 'exp') +
#'               rnorm(length(tt),0,.2)
#'
#' df <- data.frame(rho = rr, dist = tt)
#'
#' fitCorModel(tt, rr, type = 'sph')
#'
#' sol <- fitCorModel(rho ~ dist, df, nugget = NULL,
#'                    type = 'exp', hPred = tt, startp = c(200,.1))
#'
#' plot(tt,rr, ylab = 'rho', xlab = 'distance')
#' lines(sol$pred$dist, sol$pred$rho, col = 2, lwd = 2)
#'
################################################################


fitCorModel <- function(h, ...) UseMethod('fitCorModel',h)

## NOTE: range parameter is optimized in the logarithm scale
##    and nugget parameter is optimized in the logit scale
##    in order to respect bounderies.

#' @export
#' @rdname fitCorModel
fitCorModel.numeric <- function(h, x, type = 'exp', w = NULL,
                        nugget = 0, smooth = 1, hPred = NULL,
                        startp = NULL, bound = NULL, ...){

  if(is.null(w)){
    w <- 1/length(h)
  } else{
    w <- w / sum(w)
  }

  ## Case nugget unknown,
  if(is.null(nugget)){
    ofun <- function(p){
      rg <- exp(p[1])
      ng <- exp(p[2])/(1+exp(p[2]))
      return(sum(w * (x -
             corModel(c(rg,ng,smooth),h, type = type))^2))
    }

    ostart <- c(log(mean(h)),-2)

  ## Case nugget known
  } else if(!is.null(nugget)){
    ofun <- function(p){
      rg <- exp(p[1])
      return(sum(w * (x -
             corModel(c(rg,nugget,smooth),h, type = type))^2))
    }

    ostart <- log(mean(h))
  }

  ## Does a starting values are provided
  if(!is.null(startp)){
    ostart <- startp
    ostart[1] <- log(ostart[1])

    if(is.null(nugget)){
      ostart[2] <- log(ostart[2])-log(1-ostart[2])
    }
  }

  if(is.null(bound)) bound <- c(1e-4,5*max(h))

  ## Optimize accordin to number of parameter
  if(length(ostart) == 1){
    ans <- optim(ostart, ofun, method = 'Brent',
                 lower = log(bound[1]),
                 upper = log(bound[2]), ...)
  } else {
      ans <- optim(ostart, ofun, ...)

  }

  ## transform back parameter
  ans$par[1] <- exp(ans$par[1])

  if(is.null(nugget)){
    ans$par[2] <- exp(ans$par[2])/(1+exp(ans$par[2]))
  }

  if(!is.null(hPred)){
    if(is.null(nugget)){
      p0 <- c(ans$p, smooth)
    } else {
      p0 <-c(ans$p, nugget, smooth)
    }

    ans$pred <- data.frame(dist = hPred,
                           rho = corModel(p0, hPred, type = type))
  }

  return(ans)
}

#' @export
#' @rdname fitCorModel
fitCorModel.matrix <- function(x, ...){

  ans <- fitCorModel(x[,1], x[,2], type = type,
                        nugget = nugget, smooth = smooth,
                        startp = startp, method = method,
                        bound = bound, ...)
  return(ans)
}

#' export
fitCorModel.as.data.frame <- function(x, type = 'exp',
                        nugget = 0, smooth = 1,
                        startp = NULL, method = 'BFGS',
                        bound = NULL, ...){

  ans <- fitCorModel(x[,1], x[,2], type = type,
                        nugget = nugget, smooth = smooth,
                        startp = startp, method = method,
                        bound = bound, ...)
  return(ans)
}

#' @export
#' @rdname fitCorModel
fitCorModel.formula <- function(form, x, ...){

  ans <- fitCorModel(x[,all.vars(form)[2]],
                     x[,all.vars(form)[1]], ...)
  return(ans)
}


#################################################################
#' Simulation of at-site data
#'
#' Returns a matrix of the simulation of at-site data
#' based on the choice of a marginal distribution and a copula model.
#'
#' @param n Number of simulation.
#'
#' @param distr The marginal distribution. See \link{vec2par}.
#'
#' @param marg The parameter of the marginal distribution. Either a
#'   vector or a matrix with parameters in rows.
#'
#' @param cop A string representing the copula familly. Possible
#'    choice: Normal ('norm'), t-copula ('t'),
#'    Chi-square ('chisq'), squared t-copula ('tsq'),
#'    reverse chi-square ('rchisq') and
#'    reverse squared t-copula ('rtsq').
#'
#' @param h A matrix of distances.  If \code{model == 'ind'}
#'   and \code{distr == 'unif'}, a number of
#'   sites must be provided instead.
#'
#' @param type A string representing a correlation model.
#'   See \code{\link{corModel}}
#'
#' @param p The parameters of the correlation model .
#'
#' @param sigma,demi Covariance matrix or its square root.
#'
#' @param nu Degree of freedom of a t-copula (if necessary).
#'
#' @param unpivot If \code{TRUE} the results is returns in the form
#'   of a table with columns : id, time, value.
#'
#' @details
#'
#' The Chi-squared and the squared t-copula are the copulas of the
#' squared variables coming from a multivariate Normal and
#' Student distribution respectively.
#' The "squared" copulas introduce radial asymmetry and
#' their reverse have asymmetry in the opposite direction.
#'
#' @export
#'
#' @seealso \code{\link{corModel}}
#'
#' @examples
#'
#' xx <- simAtSite(100, distr = 'gev', marg = c(100, 10, .1),
#'                 cop = 'ind', h = 10, unpivot = TRUE)
#'
#' coord <- replicate(2,runif(200,0,500))
#' h <- as.matrix(dist(coord))
#'
#' # normal copula
#' xx <- simAtSite(1000, marg = c(100,20,-.1), h = h, p = c(400,0))
#'
#' hist(xx[,2])
#'
#' # t-copula
#'
#' demi <- chol(corModel(c(400,0),h))
#' xx <- simAtSite(100, distr = 'unif',
#'                 cop = 't', demi = demi, nu = 5)
#'
#' hist(xx)
#################################################################


simAtSite <- function(n, distr = 'gev', marg = c(0,1,0),
                      cop = 'norm', h = NULL,
                      type = 'exp', p = NULL,
                      sigma = NULL, demi = NULL,
                      nu = Inf, unpivot = FALSE){

  ## Computing the square root of the covariance matrix
  if(cop != 'ind'){
    if(is.null(sigma) & is.null(demi)){
      sigma <- corModel(p, as.matrix(h), type = type)
    }

    if(is.null(demi)) demi <- chol(sigma)
  }


  ## simulate from the copula
  if(cop == 'ind'){
    ans <- t(replicate(n,runif(h)))

  } else if(cop == 'norm'){
    ans <- mnormt::rmnorm(n, sqrt = demi)
    ans <- pnorm(ans)

  } else if(cop == 't'){
    ans <- mnormt::rmt(n, sqrt = demi, df = nu)
    ans <- pt(ans, df = nu)

  } else if(cop == 'chisq'){
    ans <- abs(mnormt::rmnorm(n, sqrt = demi))
    ans <- 2*pnorm(ans) - 1

  } else if(cop == 'tsq'){
    ans <- abs(mnormt::rmt(n, sqrt = demi, df = nu))
    ans <- 2*pt(ans, df = nu) -1

  } else if(cop == 'rchisq'){
    ans <- abs(mnormt::rmnorm(n, sqrt = demi))
    ans <- 2-2*pnorm(ans)

  } else if(cop == 'rtsq'){
    ans <- abs(mnormt::rmt(n, sqrt = demi, df = nu))
    ans <- 2-2*pt(ans, df = nu)

  }

  ## End the procedure if uniform marginal are selected
  if(distr == 'unif') return(ans)

  ## Transform the margins
  for(ii in seq(ncol(ans))){

    ## If marg is a vector
    if(is.vector(marg)) iiPar <- marg
    else iiPar<- marg[ii,]

    ## Remove NA parameters
    iiPar <- iiPar[!is.na(iiPar)]

    ## if a unique distribution is used
    if(length(distr) == 1) iiDistr <- distr
    else iiDistr<- distr[ii]

    ans[,ii] <- lmomco::par2qua(ans[,ii],lmomco::vec2par(iiPar,iiDistr))

  }

  ## If the simulation are not return in the form of a matrix
  if(unpivot){
    colnames(ans) <- colnames(h)
    ans <- (reshape2::melt(ans))[,c(2,1,3)]
    names(ans) <- c('id','time','value')
  }

  return(ans)
}

##################################################################
#' Draw maps of Canada
#'
#' Draw a schematic map of differents regions of Canada. It use
#' the database of the \link{maps} package with rectangular
#' projection. Best use for quick vizualization of Canadian locations.
#'
#' @param region A sting providing the region to zoom in.
#'
#' @param theme String suggesting default color parameter.
#'
#' @param ... Additional plotting parameter for map.
#'            See \code{\link{maps}} and \code{\link{plot}}
#'
#' @details
#'
#' Here the list of the preset regions :
#' 1. 'NL' : Newfoundland-Labrador
#' 2. 'MR' : Maritimes
#' 3. 'AT' : Atlantic coast
#' 4. 'QC' : Quebec
#' 5. 'SQ' : Sountern Quebec
#' 6. 'ON' : Ontario
#' 7. 'SO' : Southern Ontario
#' 8. 'ET' : East
#' 9. 'PR' : Prairies
#' 10. 'BC' : British-Columbia
#' 11. 'PC' : Pacific coast
#' 12. 'NO' : North
#' 13. 'WT' : West
#'
#' @export
#'
#' @examples
#'
#'
#' mapCan('WT', main = 'WEST CANADA',  axes = FALSE)
#'
#' mapCan('ET', col = 'gray', fill = TRUE, bg = 'lightblue',
#' col.lake = 'lightblue')
#'

mapCan <- function(region = 'canada', main = NULL, sub = NULL,
                   ylim = NULL, xlim = NULL, axes = TRUE,
                   xlab = 'Longitude', ylab = 'Latitude',
                   theme = 'white', ...){

  ## necessary as maps::map function is not fully encapsulated
  suppressPackageStartupMessages(require(maps))

  ## Convert in character in case of a number
  region <- as.character(region)

  if(theme == 'white'){
    fill = TRUE
    col = 'white'
    water = 'white'
  } else if(theme == 'bg'){
    fill = TRUE
    col = 'wheat'
    water = 'deepskyblue'
  } else{
    fill = TRUE
    col = theme
    water = 'white'
  }

  ## Select default map limits for a given region
  if(region %in% c('1', 'NL','nl',
                   'newfoundland-labrador')){
    yl <- c(46,61)
    xl <- c(-70,-50)
  }
  else if(region %in% c('2','MR','mr','maritimes')){
    yl <- c(43,50)
    xl <- c(-70,-59)
  }
  else if(region %in% c('3','AT','atl','atlantic')){
    yl <- c(42,53)
    xl <- c(-70,-50)
  }
  else if(region %in% c('4','QC','qc','quebec')){
    yl <- c(43,63)
    xl <- c(-85,-52)
  }
  else if(region %in% c('5','SQ','sq','southern-quebec')){
    yl <- c(44,52)
    xl <- c(-80,-61)
  }
  else if(region %in% c('6','ON','on','ontario')){
    yl <- c(41,57)
    xl <- c(-97,-73)
  }
  else if(region %in% c('7','SO','so','souhthern-ontario')){
    yl <- c(41,51)
    xl <- c(-97,-73)
  }
  else if(region %in% c('8','ET','et','east')){
    yl <- c(41,63)
    xl <- c(-97,-51)
  }
  else if(region %in% c('9','PR','pr','prairies')){
    yl <- c(48,61)
    xl <- c(-122,-87)
  }
  else if(region %in% c('10','BC','bc','british-columbia')){
    yl <- c(47,61)
    xl <- c(-142,-112)
  }
  else if(region %in% c('11','PC','pc','pacific')){
    yl <- c(47,70)
    xl <- c(-142,-112)
  }
  else if(region %in% c('12','NO','nor','north')){
    yl <- c(59,71)
    xl <- c(-142,-87)
  }
  else if(region %in% c('13','WT','wt','west')){
    yl <- c(47,71)
    xl <- c(-142,-87)
  }
  else{
    yl <- c(40,70)
    xl <- c(-143,-51)
  }

  ## Overrule preset limits
  if(!is.null(xlim)) xl <- xlim
  if(!is.null(ylim)) yl <- ylim

  ## Draw the maps

  maps::map(regions=c('Canada'),
           ylim = yl, xlim = xl, resolution = 0,
           col = col, fill= fill, bg = water, ...)

  maps::map(regions=c('USA'), add = TRUE, col = col,
            fill = fill, ...)

  maps::map('lakes', add = TRUE, col = water, fill = fill, ...)

  if(axes){
    maps::map.axes()
  } else {
    xlab = NULL
    ylab = NULL
  }

  title(xlab = xlab, ylab = ylab,
        main = main, sub = sub)

}


###############################################################################
#' Nearest sites
#'
#' Return a vector or martrix(in row) of the nearest sites to a target.
#'
#' @param ii Index of the target site
#' @param n Numbers of neighbors
#' @param h Matrix of distance or a distance object
#'
#' @export
#'
#' @examples
#'
#' xcoord <- replicate(2,runif(50))
#' xdis <- as.matrix(dist(xcoord))
#'
#' nearest(1,10,xdis)
#' nearest(xdis,10)
#'
#' ## if you wanted the indices only
#' t(apply(nearest(xdis,10),2,which))
#'
#' ## if the target is included in the distance, but not wanted
#' nearest(1,10,xdis, rm.target = TRUE)
#' nearest(xdis[1:5,1:5],3, rm.target = TRUE)
#'
#' ## idem ii = 0 implies that the target is zero distance
#' nearest(0, 20, xdis[2,])
#' nearest(2, 20, xdis)
#'
nearest <- function(x,...) UseMethod('nearest',x)


#' @export
#' @rdname nearest
nearest.default <- function(h,n){
  nearest.matrix(as.matrix(h,n,rm.target))
}

#' @export
#' @rdname nearest
nearest.matrix <- function(h,n, rm.target = FALSE){
  t(sapply(1:nrow(h), function(ii) nearest(ii,n,h, rm.target)))
}

#' @export
#' @rdname nearest
nearest.numeric <- function(ii,n,h, rm.target = FALSE){
  h <- as.matrix(h)
  isVector <- (ncol(h) == 1 | nrow(h) == 1)

  if(ii == 0 & isVector) ii <- which(h == 0)

  if(rm.target) n <- n+1

  if(isVector){
    ans <- seq(ncol(h)*nrow(h)) %in% order(h)[seq(n)]

  } else if(ncol(h) != nrow(h)){
    stop('The Distance Matrix is not square')

  } else  ans <- 1:nrow(h) %in% order(h[ii,])[seq(n)]

  if(rm.target) ans[ii] <- FALSE

  return(ans)
}