R/climatoKrige.R

In rxlacroix/FoRCA: French observations for R Climate Analysis

#'@title climatoKrige
#'@description Krigeage selon un parametre climato, une bounding box, une taille de cellule et une periode
#'@param d une variable climatique (character) à choisir parmi celles proposées : "NORTAV","NORTNAV","NORTXAV","NORTRAV","NORTXQ90","NORTXQ10","NORTNQ10","NORTNQ90","NORTXND","NORTNND","NORTNHT","NORTXHWD","NORTNCWD","NORTNFD","NORTXFD","NORSD","NORTR","NORHDD","NORCDD"
#'@param xmin ... copier coller la bounding box créée depuis http://bboxfinder.com (choisir les coordonnées EPSG:2154)
#'@param taille  taille du côté d'une cellule de la grille (en metres)
#'@param periode un nombre entre 1 et 12 pour choisir un mois précis
#'@export
#'@examples climatoKrige("NORTAV", 806149.1821,6294916.9104,882619.2089,6382762.4764,150,7)

climatoKrige <- function(d,xmin,ymin,xmax,ymax,taille, periode){

  #VALIDATION INPUTS

  cndata <- colnames(data_climat)

  if (paste0(d,"_",1) %in% cndata)
  {
    dataset <- data_climat
    print("Dataset : data_climat")
  } else {
    stop("Mauvais paramètre climatologique")
  }

  cat("Paramètre : ", d, "\n")

  if(xmin < 90000){
    stop("xmin trop à l'ouest, vous pouvez vous aider de http://bboxfinder.com et choisir EPSG:2154 pour calibrer vos points")
  } else if(xmax > 1278778) {
    stop("xmax trop à l'est, vous pouvez vous aider de http://bboxfinder.com et choisir EPSG:2154 pour calibrer vos points")
  } else if(ymin <6040000){
    stop("ymin trop au sud, vous pouvez vous aider de http://bboxfinder.com et choisir EPSG:2154 pour calibrer vos points")
  } else if (ymax > 7150000){
    stop("ymax trop au nord, vous pouvez vous aider de http://bboxfinder.com et choisir EPSG:2154 pour calibrer vos points")
  } else {
    # print("Grille tous les "+taille+"m : xmin = "+xmin+ "  xmax = " + xmax + "  ymin = "+ymin+ "  ymax = " + ymax)
  }

  if (is.numeric(periode)){

    data <-
      dataset %>%
      dplyr::select(Point, Latitude, Longitude,paste0(d,"_",periode))

    print(head(data))
  } else if (is.vector(periode)){

  } else if (periode == "annee"){

  } else {
    stop("Mauvais paramètre de période")
  }

  #GRRRRID
  grid <- expand.grid(x=seq(from=xmin, to=xmax, by=taille), y=seq(from=ymin, to=ymax, by=taille))
  coordinates(grid) <- ~ x+y
  proj4string(grid) <- sp::CRS("+init=epsg:2154")
  gridded(grid) <- TRUE
  cat("Grille créée tous les", taille,"m\n")

  # GEODATA
  myTheme <- rasterTheme(region = viridis(n = 10)); myTheme$axis.list$col <- 'transparent'

  data <- subset(data, data$Longitude < 1.1*xmax & data$Longitude > 0.9*xmin & data$Latitude < 1.1*ymax & data$Latitude > 0.9*ymin)
  coordinates(data) <- ~Longitude+Latitude
  proj4string(data) <- CRS("+init=epsg:2154")
  variogram <- autofitVariogram(data@data[,paste0(d,"_",periode)]~1,data, model = c("Exp","Sph","Gau","Cir","Ste","Mat","Exc", "Leg","Lin","Bes","Pen","Per","Wav","Hol","Log","Spl"), verbose = TRUE, fix.values = c(NA,NA,NA))
  plot(variogram)
  k_i <- krige(data@data[,paste0(d,"_",periode)]~1, data, grid, model=variogram$var_model, nmax =10, nmin = 2)
  k_i <- as.data.frame(k_i)
  k_i <- k_i[,-c(4)]
  r_k_i <- rasterFromXYZ(k_i)
  plot <- levelplot(r_k_i, par.settings= myTheme, main = paste0(d,"_",periode))
  print(plot)
  }