setwd("YOUR-PATH")

raw_data <- read.table(file = "data/FS_UTM.csv",
                       header = TRUE,
                       sep = ";",
                       stringsAsFactors = FALSE,
                       quote = "")

## check for issues/problems
## str(raw_data)
## apply(X=raw_data[,4:20], MARGIN = 2, FUN=function(x){x[x>1]})
raw_data$settlement[raw_data$settlement > 1] <- 1

library(sp)
sites <- raw_data[(raw_data$settlement == 1 & raw_data$iron.age == 1),]
coordinates(sites) <- ~xUTM+yUTM
proj4string(sites) <- CRS("+init=epsg:32634")

1. straight/ direct

   • Euclidean distance
   • "As the crow flies"
   • "In a beeline"
   • geodesic distance
   • great-circle distance

2. not direct/ random

   • random walk
   • "Drunkard's Walk"

\[d(i,j) = \sqrt{(x_i-x_j)^2 + (y_i-y_j)^2}\]

Distance between points using Euclidean distance using package rgeos

library(rgeos)
ed <- gDistance(spgeom1 = sites,
                byid = TRUE)
str(ed)

##  num [1:66, 1:66] 0 41669 57871 73645 74194 ...
##  - attr(*, "dimnames")=List of 2
##   ..$ : chr [1:66] "2" "22" "39" "53" ...
##   ..$ : chr [1:66] "2" "22" "39" "53" ...

dim(ed)

## [1] 66 66

head(ed[1,])

##        2       22       39       53       54       63 
##     0.00 41668.60 57871.43 73644.50 74193.66 70688.26

Can you imagine the structure if the ed object?

min(ed)

## [1] 0

round(data.frame(min = min(ed[ed>0]),
                 max = max(ed[ed>0]),
                 mean = mean(ed[ed>0]),
                 median = median(ed[ed>0])
                 ),
      2)

##    min      max     mean   median
## 1 1.73 225805.2 83029.08 79254.88

hist(ed)
abline(v = mean(ed[ed>0]), col = "red")
abline(v = median(ed[ed>0]), col = "blue")

Some further questions that might be interesting

mp <- apply(X = ed,
            MARGIN = 2,
            FUN = mean)

tmp <- ed
tmp[tmp==0] <- 9999

cp <- apply(X = tmp,
            MARGIN = 2,
            FUN = min)

sites$mp <- mp
sites@data$cp <- cp

What is the mean euclidean distance from a point to its surrounding points (mp)?

What is the closest euclidean distance from a point to its surrounding points (cp)?

Question: Why changing 0 to 9999?

spplot(sites, "cp")

spplot(sites, "mp")

Examples how to use spplot: http://rspatial.r-forge.r-project.org/gallery/

Leaflet in R: https://rstudio.github.io/leaflet/

library(leaflet)
tmp <- spTransform(sites, CRSobj = CRS("+init=epsg:4326"))
leaflet(data=tmp$cp) %>%
    addTiles() %>%
    addCircleMarkers(lng = tmp@coords[,1],
                     lat = tmp@coords[,2],
                     radius = tmp$mp/1000,
                     )

Mapview [as alternative to leaflet]: http://environmentalinformatics-marburg.github.io/mapview/introduction.html

library(mapview)
mapview(sites, zcol = "mp", legend = TRUE)

Basic work about least-cost distances and path analyses (usefulness, applicability, problems,…) was/is conducted quite excessively by Irmela Herzog (https://bodendenkmalpflege-lvr.academia.edu/IrmelaHerzog)

Technically speaking, there are two different ways how least-cost paths are created

1. create a cost raster
2. accumulate the costs from one of your points of interest
3. "walk" from a target point along the accumulated raster to your starting point (like water flows into a sink)

In GrassGIS:

• cost raster
• r.cost/r.walk
• r.drain

Technically speaking, there are two different ways how least-cost paths are created

In R:

1. create cost raster
2. create matrix of cell connections
3. find "cheapest" path between points along cell connections

gdistance in R a package by Jacob van Etten (https://cran.r-project.org/web/packages/gdistance)

What are the necessary steps?

1. Create cost surface

(https://cran.r-project.org/web/packages/gdistance/vignettes/gdistance1.pdf)

What are the necessary steps?

1. Create cost surface
2. Create transition matrix, i.e. matrix recording connections between cells

(https://cran.r-project.org/web/packages/gdistance/vignettes/gdistance1.pdf)

What are the necessary steps?

1. Create cost surface
2. Create transition matrix, i.e. matrix recording connections between cells
3. Transition matrix is filled with conductance, rather than resistance values; hence unconnected cells have value 0 (= no conductance); memory efficient

What are the necessary steps?

1. Create cost surface
2. Create transition matrix, i.e. matrix recording connections between cells
3. Transition matrix is filled with conductance, rather than resistance values; hence unconnected cells have value 0 (= no conductance); memory efficient
4. Geocorrection of values

What are the necessary steps?

1. Create cost surface
2. Create transition matrix, i.e. matrix recording connections between cells
3. Transition matrix is filled with conductance, rather than resistance values; hence unconnected cells have value 0 (= no conductance); memory efficient
4. Geocorrection of values
5. Identify adjacent cells

What are the necessary steps?

1. Create cost surface
2. Create transition matrix, i.e. matrix recording connections between cells
3. Transition matrix is filled with conductance, rather than resistance values; hence unconnected cells have value 0 (= no conductance); memory efficient
4. Geocorrection of values
5. Identify adjacent cells
6. Calculate walking speed for adjacent cells

What are the necessary steps?

1. Create cost surface
2. Create transition matrix, i.e. matrix recording connections between cells
3. Transition matrix is filled with conductance, rather than resistance values; hence unconnected cells have value 0 (= no conductance); memory efficient
4. Geocorrection of values
5. Identify adjacent cells
6. Calculate walking speed for adjacent cells
7. Geocorrect speed values

What are the necessary steps?

1. Create cost surface
2. Create transition matrix, i.e. matrix recording connections between cells
3. Transition matrix is filled with conductance, rather than resistance values; hence unconnected cells have value 0 (= no conductance); memory efficient
4. Geocorrection of values
5. Identify adjacent cells
6. Calculate walking speed for adjacent cells
7. Geocorrect speed values
8. Find path using least-cost ("as the wolf runs") or random walk ("drunkard's walk") algorithm

Now we create the transition object and geocorrect it

library(gdistance)
slope.tran <- transition(x = slope,
                         transitionFunction = mean,
                         directions = 8,
                         symm = TRUE
                         )
slope.geo <- geoCorrection(x = slope.tran, scl = TRUE)

To see how the results look like we create a raster stack for easy plotting

slopes <- stack(slope,
                raster(slope.tran),
                raster(slope.geo))
names(slopes) <- c("Slope","Transition","Geocorrection")

To see how the results look like we create a raster stack for easy plotting

plot(slopes, nr = 1)

Now, the final step. Calculate Tobler's hiking speed, geocorrect again and calculate a shortest path

adj <- adjacent(x = slope,
                cells=1:ncell(slope),
                pairs=TRUE,
                directions=8)
speed <- slope.geo
speed[adj] <-  6 * exp(-3.5 * abs(slope.geo[adj] + 0.05))

speed.geo <- geoCorrection(x = speed, scl=TRUE)

sp1 <- shortestPath(x = speed.geo,
                    origin = sites@coords[4,],
                    goal = sites@coords[50,],
                    output = "SpatialLines")

Now, the final step. Calculate Tobler's hiking speed, geocorrect again and calculate a shortest path

plot(raster(speed.geo))
lines(sp1)

What about the idea of a prehistoric autobahn? Run the code and get some coffee (most probably)

xy <- data.frame(X = sites@coords[,1],
                 Y = sites@coords[,2])

rows.xy <- row.names(xy)
sa.rows.xy <- sample(row.names(xy), length(xy$X)/2)
starts <- subset(xy, rows.xy %in% sa.rows.xy)
goals <- subset(xy, !(rows.xy %in% sa.rows.xy))
sg.passages <- brick(slope)

for(i in 1:length(starts$X)) {
    s <- c(starts$X[i],starts$Y[i])
    z <- c(goals$X[i],goals$Y[i])
    sg.passages[[i]] <- passage(x = speed.geo, origin = s, goal = z)
    cat("iteration ", i, " of ", length(starts$X),"\n")
}

## iteration  1  of  33 
## iteration  2  of  33 
## iteration  3  of  33 
## iteration  4  of  33 
## iteration  5  of  33 
## iteration  6  of  33 
## iteration  7  of  33 
## iteration  8  of  33 
## iteration  9  of  33 
## iteration  10  of  33 
## iteration  11  of  33 
## iteration  12  of  33 
## iteration  13  of  33 
## iteration  14  of  33 
## iteration  15  of  33 
## iteration  16  of  33 
## iteration  17  of  33 
## iteration  18  of  33 
## iteration  19  of  33 
## iteration  20  of  33 
## iteration  21  of  33 
## iteration  22  of  33 
## iteration  23  of  33 
## iteration  24  of  33 
## iteration  25  of  33 
## iteration  26  of  33 
## iteration  27  of  33 
## iteration  28  of  33 
## iteration  29  of  33 
## iteration  30  of  33 
## iteration  31  of  33 
## iteration  32  of  33 
## iteration  33  of  33

sg.passages.sum <- sum(sg.passages)