R/calcArea.R

In pik-piam/geodata: Landuse geodata library

## function that reads pre-calculated area information from disk
# it checks if resolution and projection agree
# if data is not available yet it will be calculated and stored
# x is a file whose res and projection are read and it is looked for the correct area

# maybe function could be much simpler: raster::area only available for latlong, and for the equal area calculation is much simpler



#' Calculate Cell Raster Area
#' 
#' Function calculates the area of a raster file. Either a raster file is
#' directly specified, or defined by it's resolution, projection and extent.
#' Method works for input in longlat or equal area projection with unit in
#' meters.
#' 
#' 
#' @usage calcArea(x=NULL, res=NULL, proj=NULL, ext=NULL, cache=T,
#' writecache=T)
#' @param x File for which the area should be calculated
#' @param res Resolution for which calculations should be done. If only one
#' number is specified, assumes square cells
#' @param proj projection as CRS. e.g. "+proj=longlat +datum=WGS84 +no_defs
#' +ellps=WGS84 +towgs84=0,0,0"
#' @param ext extent for calculations. e.g. c(-180,180,-90,90)
#' @param cache if TRUE function looks for pre-calculated area in cache
#' @param writecache if TRUE function writes newly calculated area to cache
#' folder
#' @return Either a raster or just the number for equal area projections
#' @author Ulrich Kreidenweis
#' @examples
#' 
#' \dontrun{calcArea(res=0.5, proj="+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0",
#'          ext=c(-180,180,-90,90))}
#' 
#' @export calcArea
#' @importFrom raster raster res<- extent<- projection<- area xres yres compareCRS projection writeRaster
calcArea <- function(x=NULL, res=NULL, proj=NULL, ext=NULL, cache=T, writecache=T) {
  
  if( !is.null(x) & (!is.null(res) |  !is.null(res) | !is.null(ext)) ) stop ("Either x or resulution and projection have to be specified, not both!")
  if( is.null(x) & (is.null(res) |  is.null(res) | is.null(ext)) ) stop ("Either x, or resolution and projection have to be specified")
  

  # read the data from x if it is specified
  if (!is.null(x)) {
    res <- raster::res(x)
    proj <- raster::projection(x)
    ext <- extent(x)
  }
  
  # assume the same resolution in both directions
  if (length(res)==1) {
    res <- c(res,res) 
  }
  
  # convert extent to extent object, or assume one if not specified
  if (is.null(ext)) {
    ext <- extent(c(-180,180,-90,90))
  } else if (is.vector(ext)) {
    ext <- extent(ext)
  }
  
  nameproj <- gsub("\\+proj=","",grep("+proj=", unlist(strsplit(proj, " ")), value = T))
  namedatum <- gsub("\\+datum=","",grep("+datum=", unlist(strsplit(proj, " ")), value = T))
  naming <- paste0(geodata$config$mainfolder, "/cache/calcArea/", res[1], "_", res[2], "_", nameproj, "_", namedatum, ".grd")
  
  nameunit <- gsub("\\+units=","",grep("+units=", unlist(strsplit(proj, " ")), value = T))
  
  
  ## calculate the area based on resolution, extent and projection
  .calcA <- function(res, ext, proj) {
    # for longlat use the area function from the raster package

    if (nameproj=="longlat") {
      r <- raster(crs=proj, ext=ext)
      res(r) <- res
      extent(r) <- ext
      projection(r) <- proj
      A <- area(r)
    # If data is projected in equal area and unit is meters, then calculation is easy. Just x*y.
    # In this case no raster is returned but just the value
    } else if ((nameproj=="aea" | nameproj=="laea") & nameunit=="m"){
      A <- xres(x)*yres(x)/1000000
    } else {
      stop ("No method available for this type of data")
    }
    return(A)
  }
  
  # read file from cache
  if(file.exists(naming) & cache) {
    A <- raster(naming)
    fromcache <- TRUE
    
    # if file from cache doesn't agree with requirements calculate the area newly
    if(!compareCRS(proj, projection(A)) | extent(A)!=ext) {
      cat("Projection or extent of file in cache don't aggree with input. Area is calculated")
      A <- .calcA(res, ext, proj)
      fromcache<-FALSE
    }
    
  } else {
    A <- .calcA(res, ext, proj)
    fromcache<-FALSE
  }
  
  if (writecache & !fromcache & !is.vector(A)) {
    # create the functions cache directory if it doesn't exist yet
    dir.create(file.path(geodata$config$mainfolder, "cache/calcArea"), showWarnings = FALSE)
    writeRaster(A, naming, overwrite=T)
  }
  
  attr(A, "unit") <- "km2"
  names(A) <- "km2"
  
  return(A)
}