R/makerasters.R

In g33k5p34k/PleistoDistR: Calculating Island Geomorphologies Over Pleistocene Time

#' Generate maps of Pleistocene island extents
#'
#' This function generates a series of map outputs based on the interval bins contained in the interval
#' file (see getintervals_time() and getintervals_sealvl()). This function takes a user-supplied input bathymetry raster and
#' generates a series of island extents based on the mean sea levels specified in the interval file. This function generates maps in
#' three formats: ESRI shapefile format, a flat raster format (with no topography), and a topographic raster format that preserves
#' the original bathymetric elevations of each island pixel, reprojected to the user-specified projected coordinate system.
#'
#' @param inputraster An input bathymetry raster in ASCII (.asc) format. Although PleistoDist should theoretically be able to
#' use any type of ASCII-formatted bathymetry grid as input, this tool has been tested specifically with data from the General
#' Bathymetric Chart of the Oceans (GEBCO: https://www.gebco.net). Locality-specific bathymetric maps can be downloaded from
#' https://download.gebco.net/.
#' @param epsg The projected coordinate system in EPSG code format. Because of the curvature of the Earth's surface, we need to apply a map
#' projection to accurately calculate straight-line distances between points instead of using the default WGS84 geographical coordinate system.
#' Users should specify a projected coordinate system appropriate to the geographic region being analysed using the projection's
#' associated EPSG code (https://epsg.org/home). Geographic coordinate system projections are not recommended as those will result
#' in distance matrices calculated in decimal degrees rather than in distance units.
#' @param intervalfile This is the master control file generated using either the getintervals_time() or getintervals_sealvl()
#' function. Users can also specify their own custom interval file (with nice round mean sea level values, for example)
#' although the file will need to have a column named MeanDepth and be in CSV format.
#' @param outdir The output directory for storing the maps generated by this function. If the specified output directory does not exist, PleistoDist will create
#' the folder and all necessary subsidiary folders. The makemaps() function creates three folders, one containing the output
#' maps in flat raster format (without preserving topography), one folder with maps in topographic raster format, and another folder with maps in ESRI shapefile format.
#' @param offset A correction factor to account for constant-rate tectonic uplift or subsidence in the study area. A positive offset value
#' results in PleistoDist assuming a constant rate of uplift from the past to the present (e.g. an offset of 1 causes the map for each subsequent interval to be corrected by
#' subtracting 1 * TimeInterval metres from all pixel values as PleistoDist moves backward in time). A negative offset therefore forces a constant rate of subsidence from the past to the present
#' (e.g. an offset of -1 adds 1 *TimeInterval metres to the map for each subsequent interval as PleistoDist moves backward in time). If the offset is set to the default value of 0,
#'  PleistoDist will not account for tectonic uplift or subsidence.
#' @return This function generates map outputs corresponding with the sea levels specified in the interval file. Maps are generated in flat ASCII raster (without topography),
#'  topographic ASCII raster, and ESRI shapefile formats.
#' @examples
#' \donttest{
#' #create temp directory
#' path <- file.path(tempdir())
#' #load bathymetry file
#' fiji <- system.file("extdata","FJ.asc",package="PleistoDist")
#' #generate interval file for 2 intervals and 20 kya cutoff time, binning by time
#' getintervals_time(time=20,intervals=2,outdir=path)
#' #generate maps of the Fiji Archipelago with the EPSG:3141 projection,
#' #with default interval file and a subsidence offset of 0 m per kya
#' makemaps(inputraster=fiji,epsg=3141,intervalfile=paste0(path,"/intervals.csv"),
#'     outdir=path,offset=0)
#' }
#' @export
makemaps <- function(inputraster,epsg,intervalfile,outdir,offset=0) {
  intervalfile = utils::read.csv(intervalfile)
  inputraster = terra::rast(inputraster)
  offset = as.numeric(offset)
  #create output folders
  if (base::dir.exists(outdir)==FALSE) {
    base::dir.create(outdir)
    base::dir.create(base::paste0(outdir,"/raster_topo"))
    base::dir.create(base::paste0(outdir,"/raster_flat"))
    base::dir.create(base::paste0(outdir,"/shapefile"))
  } else {
    #check to see if raster folder exists, create folder if it doesn't already exist
    if (base::dir.exists(base::paste0(outdir,"/raster_topo")) == FALSE) {
      base::dir.create(base::paste0(outdir,"/raster_topo"))
    }
    #check to see if raster_flat folder exists, create folder if it doesn't already exist
    if (base::dir.exists(base::paste0(outdir,"/raster_flat")) == FALSE) {
      base::dir.create(base::paste0(outdir,"/raster_flat"))
    }
    #check to see if shapefile folder exists, create folder if it doesn't already exist
    if (base::dir.exists(base::paste0(outdir,"/shapefile")) == FALSE) {
      base::dir.create(base::paste0(outdir,"/shapefile"))
    }
  }

  numintervals = max(intervalfile$Interval)

  for (x in 0:numintervals) {

    message("Generating maps for interval ",x,", sea level = ",intervalfile$MeanDepth[x+1])

    #apply offset value to inputraster
    inputraster_offset <- inputraster - (offset * sum(intervalfile$TimeInterval[1:(x+1)]))

    #reclassify values below interval sea level as NA
    convmat <- cbind(as.numeric(terra::global(inputraster_offset,fun="min")),intervalfile$MeanDepth[x+1],NA)
    outraster <- terra::classify(inputraster_offset,convmat,right=FALSE)

    #reproject outraster using bilinear method (since raw input is continuous)
    outraster_projected <- terra::project(outraster,base::paste0("EPSG:",epsg),method="bilinear")

    #write projected raster to raster folder
    terra::writeRaster(outraster_projected,paste0(outdir,"/raster_topo/interval",x,".tif"),filetype="GTiff",overwrite=TRUE)

    #reclassify all values as 1 to create flat raster
    convmat <- cbind(intervalfile$MeanDepth[x+1],as.numeric(terra::global(inputraster_offset,fun="max")),1)
    outraster_flat <- terra::classify(outraster_projected,convmat)

    #write flat raster to raster_flat folder
    terra::writeRaster(outraster_flat,paste0(outdir,"/raster_flat/interval",x,".tif"),filetype="GTiff",overwrite=TRUE)

    outvector <- terra::as.polygons(outraster_flat)
    outvector <- terra::disagg(outvector)

    #write reprojected vector to shapefile folder
    terra::writeVector(outvector,paste0(outdir,"/shapefile/interval",x,".shp"),filetype="ESRI Shapefile",overwrite=TRUE)
  }
  message("Done!")
}