Thursday, 24 May 2018

raster - Integrating use of surface distance when calculating cumulative cost-surface based on walking pace?


Disclaimer


This question is somewhat related to a number of earlier questions of mine:




each representing a step in my struggle to implement in R the calculation of isochrones around a given location, using the Tobler's hiking function. In what follows, I am going to ask something new, albeit it will entail making reference to what I have achieved in an attempt to address the problem I described in my earlier threads. Long story short: I hope this new question will not be flagged as possible duplicate.


Background


As said above, I was trying to calculate isochrones around a start location using the Tobler's hiking function (solved for time). I believe I eventually got something right, building on the gdistance package's documentation (https://mran.microsoft.com/snapshot/2014-11-17/web/packages/gdistance/vignettes/gdistance.pdf).


Below is a reproducible example, before I ask my new question:


First, load the gdistance library and a sample dataset (DTM):


library(gdistance)    
dtm <- raster(system.file("external/maungawhau.grd", package="gdistance"))

Secondly, calculate the terrain slope:



altDiff <- function(x){x[2] - x[1]}
hd <- transition(dtm, altDiff, 8, symm=FALSE)
slope <- geoCorrection(hd)

Thirdly, apply the Tobler's function to only those slope cells that are adjacents:


adj <- adjacent(dtm, cells=1:ncell(dtm), pairs=TRUE, directions=8)
speed <- slope
speed[adj] <- (6 * exp(3.5 * abs(slope[adj] + 0.05)))^-1 #s/m

In the above step, I used the reciprocal of the speed in m/s. Rationale: if we use speed, the final accumulated values will be 1/travel time (according to gdistance pdf guide; link above); therefore, I use the reciprocal of speed to eventually get travel time/1; travel time will be in seconds, and will require to be divided by 3600 to get hours.



Fourthly, geocorrection of the speed values (well, we should call "pace" not "speed"):


Conductance <- geoCorrection(speed)

Finally, accumulating the speed (i.e., pace) from a location:


A <- c(2667670, 6479000)
accum_sec <- accCost(Conductance, A)
accum_hour <- accum_sec/3600

What we get using rasterVis is the image further below:


levels <- seq(min(values(accum_hour)), max(values(accum_hour)), 0.15)

rasterVis::levelplot(accum_hour, par.settings = RdBuTheme, margin=FALSE, contour=TRUE, at = levels, lty=3, labels = list(cex = 0.7), label.style = 'align')

enter image description here


Given the background, here is the question


While the result looks quite good, in order to best mimic part of what is done by arcGIS's Path Distance tool (http://desktop.arcgis.com/en/arcmap/10.3/tools/spatial-analyst-toolbox/how-the-path-distance-tools-work.htm), I am wondering if is it possible to take into account surface distance in the calculation of the accumulated cost-surface, not the planar distance?


In other words, when passing from cell A to B, the cost should be (as customary) the average of A and B (since the traveller experiences half of the cost of the origin cell and half of the cost of the destination cell), but multiplied by the surface (not planar) distance (which could be quite easily derived using the Pytagorean theorem knowing the elevation difference between A and B, and their distance).


While in theory the above sounds easy, I cannot figure out how to actually integrate that to the workflow I have described earlier in this thread.




No comments:

Post a Comment

arcpy - Changing output name when exporting data driven pages to JPG?

Is there a way to save the output JPG, changing the output file name to the page name, instead of page number? I mean changing the script fo...