R/readHYDE32.R

In pik-piam/geodata: Landuse geodata library

#' Read HYDE 3.2 data
#' 
#' Reads the HYDE 3.2 baseline data in 5' resolution.
#' 
#' 
#' @usage readHYDE32(years=NULL, types=c("cropland", "grazing"), res=NULL,
#' unit="km2")
#' @param years Specific years. E.g. 1980, 1990, 2000, 2005
#' @param types Land-use types such as "cropland" or "grazing"
#' @param res if a resolution is specified, data is aggregated accordingly
#' @param unit either in 'km2' or 'Mha'
#' @return Raster stack of data
#' @author Ulrich Kreidenweis
#' @examples
#' 
#' \dontrun{readHYDE32()}
#' 
#' @export readHYDE32
#' @importFrom raster stack extent<- extent projection<-
readHYDE32 <- function(years=NULL, types=c("cropland", "grazing"), res=NULL, unit="km2") {
    
  dat <- stack()
  
  for (t in types) { 
    files <- as.list(list.files(paste0(geodata$config$mainfolder, "/HYDE/hyde32_baseline/"), pattern = t, full.names = T, ignore.case = F))
    filenames <- list.files(paste0(geodata$config$mainfolder, "/HYDE/hyde32_baseline/"), pattern = t, full.names = F, ignore.case = F)
    
    if(is.numeric(years)) {
      # only select files that contain years
      files <- files[gsub(".*([0-9]{4}).*","\\1",(filenames)) %in% years]
      if (!all(is.element(years, gsub(".*([0-9]{4}).*","\\1",(filenames))))) stop("Data for the specified years not available")  
    }
    
    tmp <- stack(files)
    names(tmp) <- paste(t, gsub(".*([0-9]{4}).*","y\\1",(names(tmp))), sep=".")
    extent(tmp) <- round(extent(tmp), digits=3) # round away small extent errors: e.g. 179.999
    projection(tmp) <- "+proj=longlat +tmpum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"
    
    dat <- stack(dat, tmp)
  }
  
  # aggregate
  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)
  }
  
  # convert unit
  if(unit=="km2") {
  attr(dat, "unit") <- "km2"
  } else if (unit=="Mha") {
    dat <- dat/10000
    attr(dat, "unit") <- "Mha"
  } else {
    print("specified unit not available, returned in km2")
  }
  
return(dat)
}