R/readHWSD.R

In pik-piam/geodata: Landuse geodata library

## function to read in data from the Harmonized World Soil Database
#' @importFrom raster stack extent<- extent projection<- setZ nlayers

## the output should be the relative distribution or the area in Mha of each slope class per 0.5 deg grid
## this area could be aggregated then within MAgPIE and the model could first convert the areas with the lower slope (possibly unrealistic, one cluster has homogenous conditions)
## therefore it should contain the median or mean slope class value

readHWSD <- function(types="slope", res=NULL) {
  
  
  ## list all files that can potentially be loaded
  
  files <- as.list(list.files(paste0(geodata$config$mainfolder, "/Harmonized World Soil Database/v1.2/"), pattern = ".asc$", full.names = T, ignore.case = F, recursive = T))
  # names(files) <- gsub("([a-z])\\.asc$", "\\1", list.files(paste0(geodata$config$mainfolder, "/Harmonized World Soil Database/v1.2/"), pattern = ".asc$", full.names = F, ignore.case = F, recursive = T))
  
  
  if(types=="slope"){
    files <- as.list(grep("GloSlopes", files, value = T))
    dat <- stack(files)
    
    # renaming for the slope classes
    n <- c("slope0_0.5perc", "slope0.5_2perc", "slope2_5perc", "slope5_10perc", "slope10_15perc", "slope15_30perc", "slope30_45perc", "slope>45perc")
    names(n) <- c("GloSlopesCl1_5min","GloSlopesCl2_5min","GloSlopesCl3_5min","GloSlopesCl4_5min","GloSlopesCl5_5min","GloSlopesCl6_5min","GloSlopesCl7_5min","GloSlopesCl8_5min")
    names(dat) <- n[names(dat)]
  } else {
    cat("Select one or several type(s). Available types are: ", paste0(names(files),collapse = "''") )
  }
  
  
  extent(dat) <- round(extent(dat), digits=3) # round away small extent errors: e.g. 179.999
  projection(dat) <- "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"
  
  
  ## the function should consider the right aggregation method for the different types
  if(!is.null(res)) {
    fact <- res/res(dat)
    if(any(fact<1)) stop("resolution has to be bigger than 0.08333333, e.g. 0.5")
    dat <- raster::aggregate(dat, fact=fact, fun=sum, expand=F)
  }
  
  setZ(dat, rep(as.Date("2000-01-01"), nlayers(dat)))
  #setZ(dat, rep(2000, nlayers(dat)))
  
  #   attr(tmp, "unit") <- "km2"
  
  return(dat)
}

## which layer has the highest value, so the highest count of cells of one class
# maxl <- which.max(dat)
# plot(maxl)