R/readSERVIR.R

In pik-piam/geodata: Landuse geodata library

## function to read the Servir data for different countries



#' Read SERVIR land cover data
#' 
#' Reads the SERVIR land cover data (http://apps.rcmrd.org/landcoverviewer/)
#' that is available for several African countries: "Botswana", "Ethiopia",
#' "Lesotho", "Malawi", "Namibia", "Rwanda", "Tanzania", "Uganda" and "Zambia"
#' 
#' 
#' @usage readSERVIR(years=NULL, types=NULL, country=NULL, rCache=T, wCache=T)
#' @param years the year for which to read the data
#' @param types land-use type
#' @param country specific country dataset
#' @param rCache TRUE: read the data from cache if available
#' @param wCache write the data to cache if it has been newly created
#' @return A raster of the SERVIR land-cover dataset
#' @author Ulrich Kreidenweis
#' @examples
#' 
#' \dontrun{readSERVIR()}
#' 
#' @export readSERVIR
readSERVIR <- function(years=NULL, types=NULL, country=NULL, rCache=T, wCache=T) {
  
  dat <- NULL
  if(rCache) dat <- readCache(); fromcache=T
  if(is.null(dat)) {
    
    name <-c("Botswana", "Ethiopia", "Lesotho", "Malawi", "Namibia", "Rwanda", "Tanzania", "Uganda", "Zambia")
    ISO <- c("BWA", "ETH", "LSO", "MWI", "NAM", "RWA", "TZA", "UGA", "ZMB")
    ISO_name <- data.frame(ISO,name)
    rownames(ISO_name) <- ISO_name$ISO
    
    if(is.null(country)) stop("One country has to be specified")
    if(length(country) > 1) stop("Only one country can be specified")
    if(nchar(country)==3) {
      country <- ISO_name[country,"name"]
      if(is.na(country)) stop(paste0("At least one country has to be specified. Available are: ",paste0(ISO_name$name,collapse = " ")))
    }
    
    filenames <- list.files(paste0(geodata$config$mainfolder, "SERVIR/"), pattern = "_Landcover_[0-9]{4}_Scheme_I.tif$", full.names = F, ignore.case = T, recursive=T)  
    
    ## only select specified country
    filenames <- filenames[gsub("^([A-Za-z]*)_.*","\\1",(filenames)) %in% country]
    
    ## 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::brick(paste0(geodata$config$mainfolder, "SERVIR/",filenames))
    names(dat) <- gsub(".*([0-9]{4}).*","y\\1", filenames) 

    ## 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$SERVIR[,c("class",t)]
        colnames(rclm) <- c("is", "becomes")
        rclm <- as.matrix(rclm)
        print(paste("SERVIR: 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)
}