1

Goal

I aim at setting to 0 (i.e., not traversable) the cells of a transition matrix (storing conductance values) that correspond to a (poly)line.

Where I am stuck

I am having hard time to figure out how to modify the values of the transition matrix.

Reproducible example

Hereunder, I am pasting a code that calculates the least-cost path (lcp) between two locations. The cost is conceived in terms of walking pace (i.e., reciprocal of speed), on the basis of the Tobler's hiking function. The code is after literature ('gdistance' package documentation).

library(gdistance)

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

heightDiff <- function(x){x[2] - x[1]}
hd <- transition(r,heightDiff,8,symm=FALSE)
slope <- geoCorrection(hd, scl=FALSE)
adj <- adjacent(r, cells=1:ncell(r), pairs=TRUE, directions=8)
speed <- slope
speed[adj] <- exp(-3.5 * abs(slope[adj] + 0.05))

x <- geoCorrection(speed, scl=FALSE)

A <- c(2667670,6479000)
B <- c(2667800,6479400)

AtoB <- shortestPath(x, A, B, output="SpatialLines")

plot(r)
plot(AtoB, add=T)

The code produces the following: enter image description here

Now, I would like to calculate the lcp between the two locations a second time, but this time setting to 0 the conductance value for the cells that correspond to the lcp previously calculated. In other words, I would like to come up with a "second best" lcp. As far as I understand, this would entail "editing" the transition matrix in order to set those cells to 0, but I am at loss of figuring out how to accomplish that.

7
  • The problem with setting these cells to 0, or 9999 maybe, is that it will create a barrier - you'll have a second path at the left or the right of the first one and not exactly a second best. An alternative I've been using is to use various neighbouring starting points and look at where/how the results switch from one path to another - which is not exactly the same either, I agree...
    – gvanhavre
    May 13 at 9:35
  • I see. But, should I want to give that a try, how can I create the barrier (i.e., setting to 0 the conductance of only those cells belonging to the 1st lcp)?
    – NewAtGis
    May 13 at 9:48
  • > library(spatialEco) > library(raster) > r2 <- rasterize(AtoB,r,background=0) > r2 <- raster.invert(r2) > r_alt <- r*r2
    – gvanhavre
    May 13 at 10:23
  • Now, that produces a clean flat zero path from point A to point B, where your next LCP will point to.
    – gvanhavre
    May 13 at 10:25
  • Maybe I am not understanding, but (using the raster produced by that procedure) the new lcp will pass exactly along the "earlier" path. The "earlier" path should be "barred" instead.
    – NewAtGis
    May 13 at 11:46

1 Answer 1

1

Transform your first AtoB path into a 0/1 raster, invert the values and create a new raster with 0 values on that path.

> library(spatialEco)
> r2 <- rasterize(AtoB,r,background=0)
> r2 <- raster.invert(r2) 
> r_alt <- r*r2

Now you can change your heightDiff function to introduce an ifelse statement considering these 0 on your path.

heightDiff_alt <- function(x){
if(x[1] == 0) {99}
else{x[2]-x[1]}
}

Then you can build a secondary path avoiding the first one.

hd_alt <- transition(r_alt,heightDiff_alt,8,symm=FALSE)

The result should show you a secondary path (on the image blow, a plain line - the first one is dotted).

Secondary path

6
  • Thanks for your reply, helpful. I have two follow-up questions: (1) can you elaborate on the rationale behind assigning 99 to the heightDiff values corresponding to the cells along the first calculated path? (2) I am wondering whether the ifelse statement that consider the 0s along the path is bound to create problems should the DTM contains terrain at/close to the sea level...not really sure about this...just thinking out loud...
    – NewAtGis
    May 15 at 16:49
  • To clarify my question (1) in the previous comment: I was wondering why 99 and not (say) 100, or 150, or 999....
    – NewAtGis
    May 15 at 16:56
  • (1) The rationale is just putting a two-digit high-enough value so that the computer won't ever dare to cross it. But it could be anything higher than your actual values, as shortestPath formula will always try to minimize the total cost. By using a very high value, you just make it absolutely sure: 99, 100, 150 or 999 are all fine. (2) Indeed, close to sea-level, there could be confusion between actual low areas and this new "no-trespassing" area. In this case, maybe one should change the r2 values, maybe by adding 999 or something...
    – gvanhavre
    May 15 at 18:30
  • But then again, what about an area with very high cliffs right on the shore? I'm afraid there won't be a straight universal answer - one has to know the area one is studying :D
    – gvanhavre
    May 15 at 18:34
  • This is why I originally thought to work on the transition layer, instead of "editing" the DTM values...To keep with your suggestion, I tried to set the "on-path" cells to NA and change the ifelse statement so that if x[1]==NA then 99, but didnt work...
    – NewAtGis
    May 16 at 8:22

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.