R/SiteExploitationTerritories.R

In eScienceCenter/SpArchR: SiteExploitationTerritories

#' converting input df of points into SpatialPointsDataFrame()
#'
#' @keywords internal

convertDFtoSPDF <- function(pts) {
  coords = cbind(pts$X, pts$Y)
  data <- data.frame(pts$id)
  points = SpatialPointsDataFrame(coords, data)
  return(points)
}


#' Site Exploitation Territories
#'
#' \code{SiteExploitationTerritories} calculates non-isotropic spatial relationship by integrating human energy expenditure in terrain based estimations. The function is based on Technical Note 3: "Distance relationships or does distance matter – calculating a non-isotropic  spatial  relationship  by  integrating  human  energy  expenditure  in  terrain based estimations – A seamless workflow for defining archaeological Site Exploitation Territories (SET) using the open source (geo-)statistical language R" by Jan Johannes Ahlrichs, Philipp Gries and Karsten Schmidt; see: http://www.uni-tuebingen.de/index.php?eID=tx_nawsecuredl&u=0&g=0&t=1499250113&hash=af99ce857ef61c6df5a2b70a907cf91473fa6deb&file=fileadmin/Uni_Tuebingen/SFB/SFB_1070/dokumente/technical_notes/technical_notes_3/CRC1070_2016_Technical_Note_03_SET.pdf
#'
#' @param pts Points
#' @param dem DEM
#' @param slope slope raster
#' @param epsg EPSG
#' @param TheoreticalSlopes a sequence
#' @param numberOfNeighbors an integer, only if no SLOPE input raster available, select algorithm to calculate slope: 4 == 4 neighbours, smooth surfaces, Fleming and Hoffer (1979), Ritter (1987) (see raster package terrain); 8 == 8 neighbours, rough surfaces, Horn (1981) (see raster package terrain)
#' @param damping TRUE/FALSE; damp hiking speed if slope is steeper than the dampingFactor (slope in degrees)
#' @param dampingFactor Damping Factor
#' @param numberOfDirections Number of Directions; move cases for transition layer; select move case for calculating time-cost; cell connections, directions for transition: 4 == 4 directions, Rook case; 8 == 8 directions, Queen case; 16 == 16 directions, knight move
#' @param timeOfInterest Time interval of interest; select a time interval of interest; use integer > 0
#' @param numberOfIsochrones Select a number and interval of isochrones output; number of isochrones [hours], integer
#' @param intervallOfIsochrones Intervals of insochrones [hours], numeric
#' @return returns a list containing for each Site a list of the Accumulated Cost Surface and Isochrones; call by using variable[[a]][[b]], while a == 1 for the first site and b == 1 for rAccumulatedCostSurface.h and b == 2 for vContourLines
#' @export

SiteExploitationTerritories <- function(pts,
                                       dem,
                                       slope = NULL,
                                       epsg,
                                       TheoreticalSlopes = seq(-70,70,1),
                                       numberOfNeighbors = 8,
                                       damping = FALSE,
                                       dampingFactor = 16,
                                       numberOfDirections = 8,
                                       timeOfInterest = 2,
                                       numberOfIsochrones = 2,
                                       intervallOfIsochrones = 1){

  WlkSpeed <- ToblersHikingFunction(TheoreticalSlopes)

  points <- convertDFtoSPDF(pts)

  #rDEM <- raster(dem)
  rDEM <- dem

  SLOPE <- slope

  # SLOPE GRADIENT
  rSLOPE <- NULL # empty slope raster variable
  if (nchar(SLOPE) > 0) { # if SLOPE is set
    rSLOPE <- raster(SLOPE) # convert SLOPE to raster datatype
  } else { # if SLOPE is not set
    # calculate SLOPE from DEM
    if(!is.na(projection(rDEM))) {
      rSLOPE <- terrain(rDEM, opt='slope', unit='degrees', neighbors=numberOfNeighbors) # calculate slope from elevation model
    } else {
      # no projection of elevation model result in projection error. Create projection file.
      print("PROJECTION ERROR: no projection is set for ELEVATION input file.")
    }
  }

  # clip raster if time limit is set
  rSLOPE4TimeCost <- rSLOPE
  rDEM4Statistics <- rDEM

  foreach(i = 1:length(points),
          .packages = c('raster',
                        'sp')) %dopar% {

    SPATIALPOINT <- data.frame(x = points@coords[[i,1]],
                               y = points@coords[[i,2]])
    coordinates(SPATIALPOINT) <- ~ x+y
    projection(SPATIALPOINT) <- projection(rDEM)

    if (timeOfInterest > 0) {

      # maximal distance in limit + time buffer 0.25h
      maxHikingDistance <- round(max(WlkSpeed)*(timeOfInterest+0.25)*1000)

      # buffer hiking distance around point of interest
      bufferMaxHikingDistance <- buffer(SPATIALPOINT, maxHikingDistance)

      #clip raster and set new extent
      rDEM_clip <- crop(rDEM,extent(bufferMaxHikingDistance))
      rSLOPE_clip <- crop(rSLOPE,extent(bufferMaxHikingDistance))

      # DEM raster for statistic section below
      rDEM4Statistics <- rDEM_clip

      # slipped slope raster for time cost calculation
      rSLOPE4TimeCost <- rSLOPE_clip
    }

    rVelocity.kmh <- calc(rSLOPE4TimeCost, ToblersHikingFunction)

    rVelocity.ms <- calc(rVelocity.kmh, fun=function(x) { ((x*1000)/3600) })

    if (damping) {
      rDamping <- rSLOPE4TimeCost
      rDamping[rDamping >  dampingFactor]  = 1000
      rDamping[rDamping <= dampingFactor]  = 1
      rVelocity.ms <- rVelocity.ms/rDamping

    }

    lTransition <- transition(rVelocity.ms, transitionFunction=mean, directions=numberOfDirections)

    lGeoCorrection <- geoCorrection(lTransition, type="r")
    rAccumulatedCostSurface.s <- accCost(lGeoCorrection, SPATIALPOINT)

    rAccumulatedCostSurface.h <- rAccumulatedCostSurface.s/3600

    vContourLines <- rasterToContour(rAccumulatedCostSurface.h, levels=seq(0, numberOfIsochrones, intervallOfIsochrones))

    zonalStatistics   <- data.frame(matrix(0,2,5))
    statNames <- c('1st hour','min','max','mean','sd')
    names(zonalStatistics) <- statNames
    zonalStatistics[1,1] <- 'DEM'
    zonalStatistics[2,1] <- 'SLOPE'

    rasterZones <- rAccumulatedCostSurface.h
    rasterZones[rAccumulatedCostSurface.h <= 1] = 1 # pixels of first hour
    rasterZones[rAccumulatedCostSurface.h > 1]  = 2 # pixels more than first hour

    for(st in 2:length(statNames)) {
      zonalStatistics[1,st] <- zonal(rDEM4Statistics, rasterZones, statNames[st])[1,2]
      zonalStatistics[2,st] <- zonal(rSLOPE4TimeCost, rasterZones, statNames[st])[1,2]
    }

    TCR = paste("Site", i, "TimeCostRaster", sep = "_") # name of time cost raster
    SLG = paste("Site", i, "SlopeGradient", sep = "_") # name of slope gradient raster
    CTL = paste("Site", i, "ContourLines", sep = "_") # name of contour lines shapefile
    rdt = "ascii" # datatype of output raster files, for more information use 'help(writeRaster)'


    # rSLOPE4TimeCost ist in allen Iterationen gleich!

    results <- c(rAccumulatedCostSurface.h, vContourLines)

    return(results)
  }
}