R/readHYDE31.R

In pik-piam/geodata: Landuse geodata library

## read HYDE 3.1 data

## to do:
# include check whether years are within the data
# make a selection of types possible



#' Read HYDE 3.1 cropland and grassland data
#' 
#' Reads the HYDE 3.1 cropland and grassland data in 5' resolution. Grasland
#' and cropland values are available. From 1961 on data should be consistent
#' with FAO data.
#' 
#' 
#' @usage readHYDE31(years=NULL, types=c("crop", "gras"), res=NULL, unit="km2")
#' @param years Specific years. E.g. 1980, 1990, 2000, 2005
#' @param types Land-use type: "gras" or "crop"
#' @param res if a resolution is specified, data is aggregated accordingly
#' @param unit either in 'km2' or 'Mha'
#' @return If one or more types are specified a list of RasterStacks (of
#' different years) with different types. Otherwise just a RasterStack
#' @author Ulrich Kreidenweis
#' @examples
#' 
#' \dontrun{readHYDE31()}
#' 
#' @export readHYDE31
#' @importFrom raster stack extent<- extent projection<-
readHYDE31 <- function(years=NULL, types=c("crop", "gras"), res=NULL, unit="km2") {
    
  tmp <- list()
  
  for (t in types) { 
    files <- as.list(list.files(paste0(geodata$config$mainfolder, "/HYDE/hyde31/"), pattern = t, full.names = T, ignore.case = F))
    filenames <- list.files(paste0(geodata$config$mainfolder, "/HYDE/hyde31/"), 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")  
    }
    
    dat <- stack(files)
    names(dat) <- gsub(".*([0-9]{4}).*","y\\1",(names(dat)))
    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"
    
    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)
    }

    tmp[t] <- dat
  }
  
  # convert unit
  if(unit=="km2") {
  attr(tmp, "unit") <- "km2"
  } else if (unit=="Mha") {
    tmp <- lapply(tmp, function(x) x/10000)
    attr(tmp, "unit") <- "Mha"
  } else {
    print("specified unit not available, returned in km2")
  }
  
  
  # only a list of RasterStacks if there is more than one type.
  if (length(types)==1) tmp <- tmp[[1]]
  
  
return(tmp)
}