R/oEFTA.R

In EhrmannS/rasterTools: Obtain and process earth observation data

#' Obtain data from the European Atlas of Forest Tree Species
#'
#' Obtain presence and habitat suitability
#' \href{http://forest.jrc.ec.europa.eu/european-atlas-of-forest-tree-species/atlas-data-and-metadata/}{data}
#' from the
#' \href{http://forest.jrc.ec.europa.eu/european-atlas-of-forest-tree-species/}{European
#' Atlas of Forest Tree Species}
#' @param mask [\code{geom} | \code{Spatial*} | \code{sf} | \code{raster}]\cr
#'   spatial object of which the extent is the area of interest.
#' @param species [\code{character(.)}]\cr name(s) of species to get data for.
#'   Can be abbreviated if a \code{\link{catalog}} is provided.
#' @param type [\code{character(1)}]\cr the dataset type, either \code{"rpp"}
#'   (relative probability of presence) or \code{"mhs"} (maximum habitat
#'   suitability).
#' @return bla
#' @references de Rigo, D., Caudullo, G., Houston Durrant, T., San-Miguel-Ayanz,
#'   J., 2016. The European Atlas of Forest Tree Species: modelling, data and
#'   information on forest tree species. In: San-Miguel-Ayanz, J., de Rigo, D.,
#'   Caudullo, G., Houston Durrant, T., Mauri, A. (Eds.), European Atlas of
#'   Forest Tree Species. Publ. Off. EU, Luxembourg, pp.
#' @details The values in this dataset are originally stored as values in the
#'   range of 0 and 1 (the probability of presence). Here they are transformed
#'   to integer factors from 0 to 100.
#'   
#'   The object provided in \code{mask} is treated as a single mask,
#'   irrespective of that object consisting of only one or several features. The
#'   extent comprising all features (point(s), line(s), polygon(s)) is used as
#'   area of interest. This is in contrast to \code{\link{obtain}}, where a mask
#'   may consist of several features, each of which are treated as seperate
#'   mask.
#' @family obtain operators (Europe)
#' @examples
#' \dontrun{
#'
#' myTrees <- oEFTA(mask = rtGeoms$mask,
#'                  species = c("Quercus robur", "Alnus incana", 
#'                              "Pinus sylvestris", "Betula sp"))
#' visualise(gridded = myTrees, trace = TRUE)
#'
#' # get the (updated) bibliography
#' reference(style = "bibtex")
#' }
#' @importFrom geometr getCRS gs_rectangle getExtent setCRS
#' @importFrom checkmate testClass
#' @importFrom raster crop unique stack
#' @export

oEFTA <- function(mask = NULL, species = NULL, type = "rpp"){
  
  # check arguments
  maskIsGeom <- testClass(mask, classes = "geom")
  maskIsSp <- testClass(mask, classes = "Spatial")
  maskIsSf <- testClass(mask, classes = "sf")
  assert(maskIsGeom, maskIsSp, maskIsSf)
  speciesIsDF <- testDataFrame(species, any.missing = FALSE, ncols = 2, min.rows = 1, col.names = "named")
  if(speciesIsDF){
    assertNames(names(species), must.include = c("original", "abbr"))
    species <- species$original
  } else{
    assertCharacter(species)
  }
  assertSubset(type, choices = c("rpp", "mhs"))
  
  if(type == "rpp"){
    steps <- c(5, 5, 20, 20, 20, 20, 10)
    labels <- rep(c("marginal", "low", "mid-low", "medium", "mid-high", "high", "very-high"), steps)
    efta_palette <- rtPalette(colors = c("#fcf0d400", "#e8f5c3ff", "#bfe361ff", "#7abd2aff", "#438532ff", "#16301bff", "#050707ff", "#050707ff"), 
                              steps = steps)
    outCols <- c(efta_palette(100), rep("#000000", 155))
  } else{
    # steps <- c(, , , , , )
    # labels <- rep(c("negligible", "low", "mid-low", "medium", "mid-high", "high"), steps)
    # efta_palette <- rtPalette(colors = c("#", "#", "#", "#", "#", "#", "#", "#"), 
    #                           steps = steps)
    # outCols <- c(efta_palette(100), rep("#000000", 155))
  }
  
  # transform crs of the mask to the dataset crs
  targetCRS <- getCRS(x = mask)
  maskExtent <- getExtent(x = mask)
  if(targetCRS != projs$laea){
    targetMask <- setCRS(x = mask, crs = projs$laea)
  } else{
    targetMask <- mask
  } 
  maskGeom <- gs_rectangle(anchor = getExtent(x = targetMask))
  maskGeom <- setCRS(x = maskGeom, crs = projs$laea)
  targetExtent <- getExtent(maskGeom)

  # determine the species that need to be dealt with
  if(any(!species %in% meta_efta$species)){
    warning(paste0("species '", species[which(!species %in% meta_efta$species)], "' is not part of this dataset"))
    species <- species[species %in% meta_efta$species]
  }
  if(type == "rpp"){
    dataExists <- meta_efta$rpp[meta_efta$species %in% species]
  } else{
    dataExists <- meta_efta$mhs[meta_efta$species %in% species]
  }
  
  # go through 'species' to extract data
  efta_out <- stack()
  history <- list()
  for(i in seq_along(species)){
    thisSpecies <- species[i]
    fileName <- paste0(sub(thisSpecies, pattern = " ", replacement  = "-"), "_", type, ".tif")
    fileExists <- testFileExists(paste0(rtPaths$efta$local, "/", fileName))
    message(paste0("I am handling the tree species '", thisSpecies, "':"))
    
    if(any(!fileExists)){
      downloadEFTA(file = fileName,
                  localPath = rtPaths$gfc$local)
    }    
    
    temp <- strsplit(fileName, "_")[[1]]
    temp2 <- strsplit(temp[2], "[.]")[[1]]
    shortName <- paste0(temp[1], "_", temp2[1])
    tempObject <- gdalwarp(srcfile = paste0(rtPaths$efta$local, "/", fileName),
                           dstfile = paste0(rtPaths$project, "/", shortName, "_", paste0(round(maskExtent$x), collapse = "."), "_", paste0(round(maskExtent$y), collapse = "."), ".tif"),
                           s_srs = projs$laea,
                           t_srs = targetCRS,
                           te = c(maskExtent$x[1], maskExtent$y[1], maskExtent$x[2], maskExtent$y[2]),
                           overwrite = TRUE,
                           output_Raster = TRUE)
    tempObject <- round(tempObject*100)
    
    # manage the object history
    history <- c(history, paste0("object loaded from '", fileName, "'"))
    history <-  c(history, paste0("object cropped between points (x, y) '", targetExtent$x[1], ", ", targetExtent$y[1], "' and '", targetExtent$x[2], ", ", targetExtent$y[2], "'"))
    if(targetCRS != projs$laea){
      crs_name <- strsplit(targetCRS, " ")[[1]][1]
      history <- c(history, list(paste0("object reprojected to ", crs_name)))
    }
    
    # create and set RAT table
    tempObject@data@isfactor <- TRUE
    ids <- unique(tempObject)
    if(type == "rpp"){
      tempObject@data@attributes <- list(data.frame(id = ids, presence = labels[ids+1]))
    } else{
      tempObject@data@attributes <- list(data.frame(id = ids, suitability = labels[ids+1]))
    }
    tempObject@history <- history
    
    # set colortable
    tempObject@legend@colortable <- outCols
    
    names(tempObject) <- sub(thisSpecies, pattern = " ", replacement = "_")
    efta_out <- stack(efta_out, tempObject)
  }

  # manage the bibliography entry
  bib <- bibentry(bibtype = "InBook",
                  title = "The European Atlas of Forest Tree Species: modelling, data and information on forest tree species",
                  author = c(person(given = "Daniele", family = "de Rigo"),
                             person(given = "Giovanni", family = "Caudullo"),
                             person(given = "Tracy", family = "Houston Durrant"),
                             person(given = "Jes\uFAs", family = "San-Miguel-Ayanz")
                  ),
                  chapter = 2,
                  year = 2016,
                  booktitle = "European Atlas of Forest Tree Species",
                  editor = c(person(given = "Jes\uFAs", family = "San-Miguel-Ayanz"),
                             person(given = "Daniele", family = "de Rigo"),
                             person(given = "Giovanni", family = "Caudullo"),
                             person(given = "Tracy", family = "Houston Durrant"),
                             person(given = "Achille", family = "Mauri")
                  ),
                  publisher = "Publ. Off. EU",
                  address = "Luxembourg"
  )

  if(is.null(getOption("bibliography"))){
    options(bibliography = bib)
  } else{
    currentBib <- getOption("bibliography")
    if(!bib%in%currentBib){
      options(bibliography = c(currentBib, bib))
    }
  }
  
  return(efta_out)
  
}

#' @describeIn oEFTA function to download data related to the EFTA dataset
#' @param file [\code{character(1)}]\cr the name of the file to download.
#' @param localPath [\code{character(1)}]\cr the local path where files are
#'   supposed to be stored (i.e. on your harddisc).
#' @importFrom httr GET write_disk progress
#' @export

downloadEFTA <- function(file = NULL, localPath = NULL){
  
  assertCharacter(file, any.missing = FALSE, len = 1, null.ok = TRUE)
  if(!is.null(localPath)){
    assertDirectory(localPath, access = "rw")
  }
  
  if(!is.null(file) & !is.null(localPath)){
    
    # https://w3id.org/mtv/FISE/map-data-RPP/v0-3-2/internet/Abies-alba
    # https://w3id.org/mtv/FISE/map-data-MHS/v0-3-2/internet/Abies-alba
    
    onlinePath <- rtPaths$efta$online
    message(paste0("  ... downloading the file from '", onlinePath, "'"))
    
    GET(url = paste0(onlinePath, file),
        write_disk(paste0(localPath, "/", file)),
        progress())
    
    message(paste0("  ... unzipping the files of '", file, "'"))
    unzip(paste0(localPath, "/", file), exdir = localPath)
  }
}