R/readNLCD.R

In pik-piam/geodata: Landuse geodata library

### read landcover data for the US from the NLCD

#notes:
# especially linear elements like roads are lost if data is converted to tif with WGS projection in ArcGIS.


# TO DO:
# change in geodata$read$NCLD1992 pasture to past!!



#' Read National Land Cover Dataset of the USA
#' 
#' Reads the National Land Cover Dataset of the Mulit-resolution Land
#' Characteristics Consortium (MRLC) for the USA.
#' 
#' 
#' @usage readNLCD(years=NULL, types=NULL, reproject=F, multicore=T)
#' @param years Specify the year the data is to be read. Available are 1992,
#' 2001, 2006 and 2011. This does not necessarily equal the year the data is
#' representative for.
#' @param types Land-cover type. If specified only
#' @param reproject If TRUE, data is reprojected to WGS84 (lat/long)
#' @param multicore If TRUE, several cores are used to reclassify the data
#' @return If one or more types are specified a list of raster stacks.
#' Otherwise simply a raster stack
#' @author Ulrich Kreidenweis
#' @examples
#' 
#' \dontrun{readNLCD()}
#' 
#' @export readNLCD
#' @importFrom raster raster
#' @importFrom parallel makeCluster detectCores
readNLCD <- function(years=NULL, types=NULL, reproject=F, multicore=T) {
  
  if(length(years)>1 & any(years==1992)) stop("Currently only one year supported")
  
  filenames <- list.files(paste0(geodata$config$mainfolder, "/NLCD/"), pattern = ".img$", full.names = F, ignore.case = F, recursive=T)  
  
  ## only select specified years: This doesn't work for 1992
  if(years==1992) {
    # this is a version that has been manually reprojected with ArcGIS to WGS84. Reprojection with R on cluster failed several times
    # consider reprojecting all files with ArcGIS. However some detail is lost. Therefore better do this later
    dat <- raster(paste0(paste0(geodata$config$mainfolder, "/NLCD/1992/nlcd_1992_30meter_whole_orig.tif")))
  } else if(is.numeric(years)) {
      filenames <- filenames[gsub(".*nlcd_([0-9]{4})_landcover.*","\\1",(filenames)) %in% years]
     if (!all(is.element(years, gsub(".*nlcd_([0-9]{4})_landcover.*","\\1",(filenames))))) stop("Data for the specified years not available")  
     
     dat <- raster::stack(paste0(geodata$config$mainfolder, "/NLCD/", filenames))
     names(dat) <- gsub(".*nlcd_([0-9]{4})_landcover.*","y\\1", filenames) 
     
  } else stop("year not set correctly")
  
    
  ## reproject the raster: very slow: consider doing this only after everything has been processed
#   if (reproject) {
#     print("Reprojecting dataset now")
#     newproj <- "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"
#     dat <- raster::projectRaster(dat, crs=newproj)
#   }
  
  ## if a certain type is selected reclassify it according to the rules
  if (!is.null(types)){ 
    out <- list()
    
    for (t in types){
      
      if (years == 1992) {
        rclm <-  geodata$read$NLCD1992[,c("class",t)]  
      } else {
        rclm <-  geodata$read$NLCD[,c("class",t)] 
      }
      colnames(rclm) <- c("is", "becomes")
      rclm <- as.matrix(rclm)
      print(paste("ESA CCI: Reclassify", t, "now"))
          
      if (multicore) {
#        raster::beginCluster(n=geodata$config$default$ncluster)
        raster::beginCluster() # try without a limitation of nodes
        no_cores <- detectCores() - 1
        cl <- makeCluster(no_cores)
        
        out[t] <- raster::clusterR(dat, fun=reclassify, args = list(rcl=rclm))
        raster::endCluster()
#        stopCluster(cl)

      } else {
        out[t] <- raster::reclassify(dat, rcl=rclm)
      }
      attr(out, "unit") <- "share"
    }
    
  } else {
    out <- dat
    attr(out, "unit") <- "class"
  } 
  
  if (length(types)==1) out <- out[[1]]
  
  return(out)
}