R/readEsaCci.R

Defines functions readEsaCci

Documented in readEsaCci

### function to read ESA CCI: European Space Agency, Climate Change Initiative, Land_Cover_CCI



#' Reads ESA CCI Land Cover Maps
#' 
#' This function reads ESA Climate Change Initiative Land Cover Maps for
#' different years. It is also possible to select a certain type of land-cover
#' that is then extracted and transferred to percentage values based on a
#' substitution list.
#' 
#' 
#' @usage readEsaCci(years=NULL, types=NULL, rCache=T, wCache=T)
#' @param years Specify the years that you want the data for. Currently
#' available: 2000 (1998-2002), 2005 (2003-2007), 2010 (2008-2012)
#' @param types If specified a certain land-cover type is extracted and
#' transferred to percentage values in original resolution. Available types
#' see: geodata$class2area$EsaCci, e.g. "tree.avg", "water", ...
#' @param rCache TRUE: read the data from cache if available
#' @param wCache write the data to cache if it has been newly created
#' @return If one or more types are specified a list of raster stacks (of
#' rasters of different years) with different types. Otherwise a raster stack
#' @author Ulrich Kreidenweis
#' @examples
#' 
#' \dontrun{readEsaCci()}
#' 
#' @export readEsaCci
readEsaCci <- function(years=NULL, types=NULL, rCache=T, wCache=T) {
    
  dat <- NULL
  if(rCache) dat <- readCache(); fromcache=T
  if(is.null(dat)) {
  
    filenames <- list.files(paste0(geodata$config$mainfolder, "CCI/"), pattern = "ESACCI-LC-L4-LCCS-Map-300m-P5Y-[0-9]{4}-v1\\.6\\.1\\.tif$", full.names = F, ignore.case = F, recursive=F)  
    
    ## only select specified years
    if(is.numeric(years)) {
      filenames <- filenames[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")  
    }
    
    ## rename stack dimensions to years
    dat <- raster::stack(paste0(geodata$config$mainfolder, "CCI/",filenames))
    names(dat) <- gsub(".*([0-9]{4}).*","y\\1", filenames) 
    raster::extent(dat) <- round(raster::extent(dat), digits=3)
    
    ## if a certain type is selected reclassify it according to the rules
    if (!is.null(types)){ 
      out <- list()
        
      for (t in types){
        rclm <-  geodata$read$EsaCci[,c("class",t)]
        colnames(rclm) <- c("is", "becomes")
        rclm <- as.matrix(rclm)
        print(paste("ESA CCI: Reclassify", t, "now"))
        out[t] <- raster::reclassify(dat, rcl=rclm)
  
        ## version of reclassify based on several cores
  #      raster::beginCluster(n=geodata$config$default$ncluster)
  #      out[t] <- raster::clusterR(dat, fun=reclassify, args = list(rcl=rclm))
        
  #      out[t] <- raster::subs(dat, y=rcl, by=1, which=2) # according to help subs should be faster, but failed locally
        attr(out, "unit") <- "share"
      }
    } else {
      out <- dat
      attr(out, "unit") <- "class"
    } 
    
    if (length(types)==1) out <- out[[1]]
  }
  
  if (wCache & !fromcache) writeCache(out)
  
  return(out)
}
pik-piam/geodata documentation built on Nov. 5, 2019, 12:21 a.m.