R/spatial_datahandling.R

In Martin-Jung/naturepackage: Functions to run naturemap on a new computer

#### General handling scripts of spatial data ####
#' Convert Longitude / Latitude to UTM zone
#'
#' @author Martin Jung
#' @param lon A vector of longitude coordinates
#' @param lat A vector of latitude coordinates
#' @return Gives the UTM zone number for a given lon-lat coordinate
#' @export
latlong2UTMzone <- function(lon,lat){
  assert_that(
    length(lon) == length(lat),msg = "Provide an equal number of observations for both!"
  )
  # Normal calculation
  ZoneNumber = floor((lon + 180)/6) + 1

  # Special case for higher longitude levels
  if( lat >= 56.0 && lat < 64.0 && lon >= 3.0 && lon < 12.0 ){
    ZoneNumber = 32
  }
  # Special cases for svalbard
  if( lat >= 72.0 && lat < 84.0 ) {
    if  ( lon >= 0.0  && lon <  9.0 ) ZoneNumber = 31
    else if( lon >= 9.0  && lon < 21.0 ) ZoneNumber = 33
    else if( lon >= 21.0 && lon < 33.0 ) ZoneNumber = 35
    else if( lon >= 33.0 && lon < 42.0 ) ZoneNumber = 37
  }

  #Return the result
  return(ZoneNumber)
}

#' Save a compressed Geotiff file
#'
#' @author Martin Jung
#' @param file A Raster* object
#' @param fname The output file name of the rasterlayer
#' @param dt The data type of the output (Default: 'INT2S')
#' @return returns the unified total extent
#' @import raster
#' @import assertthat
#' @export
writeGeoTiff <- function(file, fname,dt = "INT2S"){
  if(!has_extension(fname,"tif")) {fname <- paste0(fname,".tif")}
  writeRaster(file,fname,
              format='GTiff',
              datatype = dt,
              NAflag = -9999,
              options=c("COMPRESS=DEFLATE","PREDICTOR=2","ZLEVEL=9"),
              overwrite= TRUE )
}

#' Create a empty raster based on a provided layer
#' @param x An existing Raster object
#' @return A RasterLayer with the same dimensions as the input
#' @author Martin Jung
#' @import raster
#' @export
emptyraster <- function(x, ...) { # add name, filename,
  emptyraster <- raster(nrows=nrow(x), ncols=ncol(x),
                        crs=x@crs,
                        ext=extent(x), ...)

  return(emptyraster)
}

#' Unify extent of a rasterstack
#'
#' @author Martin Jung
#' @param rList a list of raster objects
#' @return returns the unified total extent
#' @import raster
#' @export
max_extent <- function(rlist){
  # given list of rasters
  # returns union of extent
  xmin=min(sapply(rl,FUN=function(x){extent(x)@xmin}))
  xmax=max(sapply(rl,FUN=function(x){extent(x)@xmax}))
  ymin=min(sapply(rl,FUN=function(x){extent(x)@ymin}))
  ymax=max(sapply(rl,FUN=function(x){extent(x)@ymax}))

  extent(c(xmin,xmax,ymin,ymax))
}

#' Expand an extent by a certain number
#'
#' @author Martin Jung
#' @param e an extent object
#' @param f value to increase the extent (Default = 0.1)
#' @return returns the unified total extent
#' @export
extent_expand <- function(e,f=0.1){

  xi <- (e@xmax-e@xmin)*(f/2)
  yi <- (e@ymax-e@ymin)*(f/2)

  xmin <- e@xmin-xi
  xmax <- e@xmax+xi
  ymin <- e@ymin-yi
  ymax <- e@ymax+yi

  return(extent(c(xmin,xmax,ymin,ymax)))
}

#' Normalize a raster by it's range
#'
#' @author Martin Jung
#' @param x A target raster
#' @return Returns a normalized raster layer with range 0-1
#' @import raster
#' @export
normalize_raster <- function(x) {
  min <- raster::minValue(x)
  max <- raster::maxValue(x)
  return((x - min) / (max - min))
}

#' Allign raster data by bringing it in the same geometry and extend.
#'
#' If the data is not in the same projection as the template, the allignment
#' will be computed by reprojection only. If the data has already the same
#' projection, the data set will be croped and aggregated prior to resampling
#' in order to reduce computation time.
#'
#' @author Martin Jung
#' @author Thomas Nauss
#' @import raster
#' @param data raster layer to be resampled
#' @param template raster or spatial data set from which geometry can be extracted
#' @param method method for resampling ("ngb" or "bilinear")
#' @param cl Should cluster computation be used (Default=T)
#' @import raster
#' @return raster layer containing geometrically alligned data
#' @export

alignRasters <- function(data, template, method = "bilinear",func = mean,cl = T,expand = TRUE){
  lib <- c("raster")
  sapply(lib, function(...) stopifnot(require(..., character.only = T)))
  if(cl) beginCluster(parallel::detectCores()-1)
  if(projection(data) == projection(template)){
    data <- raster::crop(data, template, snap = "near")
    if(class(template) == "RasterLayer"){
      if(data@ncols / template@ncols >= 2){
        factor <- floor(data@ncols/template@ncols)
        data <- raster::aggregate(data, fact = factor, fun = func,na.rm = TRUE,
                                  expand=expand)
      }
      data <- raster::resample(data, template, method = method)
    }
  } else {
    data <- projectRaster(data, template, method = method)
  }
  if(cl) endCluster()
  return(data)
}

#' Split a given raster layer into multiple parts
#' This is done in parallel
#'
#' @author Martin Jung
#' @param infile The path to an input raster layer
#' @param outpath The output directory
#' @param parts Into how many input parts should the layer be split
#' @export
split_raster <- function(infile, outpath, parts=3) {
  ## parts = division applied to each side of raster, i.e. parts = 2 gives 4 tiles, 3 gives 9, etc.
  lib <- c("gdalUtils","parallel","reshape2")
  sapply(lib, function(...) stopifnot(require(..., character.only = T)))
  # Check if infile exists
  stopifnot(file.exists(infile))

  filename <- strsplit(basename(infile),".tif")[[1]]
  # Get the dimensions of the file
  dims <- as.numeric(
    strsplit(gsub('Size is|\\s+', '', grep('Size is', gdalinfo(infile), value=TRUE)),
             ',')[[1]]
  )

  # Generate window frames
  xy <- list()
  # t is nr. of iterations per side
  t <- parts - 1
  for (i in 0:t) {
    for (j in 0:t) {
      # [-srcwin xoff yoff xsize ysize] src_dataset dst_dataset
      srcwin <- paste(i * dims[1]/parts, j * dims[2]/parts, dims[1]/parts, dims[2]/parts)
      xy[[paste0(i,"_",j)]] <- data.frame(infile = infile,srcwin = srcwin, file = paste0(outpath,"/",filename,"_",i,"_",j,".tif"))
    }
  }
  df <- reshape2::melt(xy)

  # Then process per src_win
  cat("Start splitting: ",filename)

  # Create a function to split the raster using gdalUtils::gdal_translate
  split <- function(input, outfile, srcwin) {
    gdalUtils::gdal_translate(input, outfile, srcwin=srcwin)
  }

  # Make a copy for export
  df_org <- df
  # Kick out files already existing
  df <- df[which(!file.exists(as.character(df$file))),]

  # Make cluster
  cl <- makeCluster(detectCores()-1)
  clusterExport(cl, c('split', 'df'))
  clusterEvalQ(cl,library(gdalUtils))

  system.time({
    parLapply(cl, seq_len(nrow(df)), function(i) {
      split(df$infile[i], df$file[i], df$srcwin[i])
    })
  })
  stopCluster(cl)
  cat("\n")
  cat("Done")
  return(df_org)
}

#' Project a given rasterlayer to Mollweide projection
#'
#' @param image The path to an input raster image
#' @param output The path to an output filename
#' @return The RasterLayer object just created
#' @author Martin Jung
#' @export
projectMollWeide <- function(image, output = "test.tif",s.crs = "+proj=moll +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +units=m +no_defs"){
  if(output == "" || is.null(output)) { stop('Specify an output filename path!')}
  if(s.crs == "" || is.null(output)) { stop('Specify an output projection!')}
  out <- gdalwarp(srcfile = image@file@name,
                  dstfile = output,
                  s_srs = proj4string(image),
                  t_srs = s.crs,
                  r = "near",
                  co = c("COMPRESS=DEFLATE","PREDICTOR=2","ZLEVEL=9"),
                  output_Raster=TRUE, multi = TRUE,
                  overwrite=TRUE,verbose=TRUE)
  return( out )
}