-
Notifications
You must be signed in to change notification settings - Fork 0
/
SiteExploitationTerritories.R
149 lines (114 loc) · 6.74 KB
/
SiteExploitationTerritories.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
#' 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)
}
}