[R-sig-Geo] gDistance transition objects: dealing with NAs in source cost raster

[T C Wilkinson]

Dear all,

I’m using two functions (copied below, based on Jacob van Etten’s own documentation for distance) in order to construct the required transition objects for use in gDistance functions such as shortestPath and accCost, using an input DEM (digital elevation model).

The particular DEM I am using includes an area of sea which is marked as NA and, for my particular analyses, should be ignored (i.e. should be modelled as 0 conductance) in the distance calculations.  My assumption was that these functions would ignore the area of NA-marked sea by default (in part because the equivalent Arc- models exclude NoData areas).  

However, a plot of the transition as raster `plot(raster(conductance))` suggests the sea is included as a low-level conductance (it appears as pink from the col.palette and not black as specified in colNA field).  Similarly, subsequently-created accCost rasters include the sea as part of the accumulate cost surface and not, as I hoped, as NA.

I tried a different altDiffFunction, in an attempt to exclude the NA values, but I suspect this is merely a less elegant method of producing the same result.

Whilst I wondered if some kind of transitionStack might solve the problem, my (albeit limited) understanding of the way gDistance transition layers work suggest that this is overkill for what is essentially a single-variable (slope) function.

Does the resolution of the DEM play any role here?

Any suggestions for how to solve this (i.e. how to exclude the NA area from the transition model) would be very gratefully received!

Many thanks,
Toby




**CODE EXTRACT**:

sw_slope <- function(dem, dirs=16, symm=FALSE){
  #calculate the altitudinal differences between cells
  #altDiffFunction <- function(x){x[2] - x[1]}
  altDiffFunction <- function(x){
    return(ifelse(((!is.na(x[2])) & (!is.na(x[1]))),(x[2] - x[1]),NA))
  }
  
  heightDiff <- transition(dem, altDiffFunction, directions=dirs, symm=symm)
  #divide by the distance between cells
  slope <- geoCorrection(heightDiff)
  return(slope)
} 



sw_conductance <- function(dem, dirs=16, symm=FALSE){
  #first find slope
  slope <- sw_slope(dem, dirs, symm)
  
  #create an index for adjacent cells (adj) with the function adjacent
  adj <- adjacent(dem, cells=1:ncell(dem), pairs=TRUE, directions=16)
  cost <- slope
  
  #extract and replace adjacent cells
  #note the 1/ to convert from resistivity (as cited in publication above) to conductivity
  cost[adj] <- 1 / ((((1337.8 * cost[adj]^6) + (278.19 * cost[adj]^5) - (517.39 * cost[adj]^4) - (78.199 * cost[adj]^3) + (93.419 * cost[adj]^2) + (19.825 * cost[adj]) + 1.64)))
  
  #correction to take account of distance between cell centres
  gc <- geoCorrection(cost)
  
  return(gc)
}


dem <- raster() # the dem, based on an reprojected aster gdem tile, is loaded from file

conductance <- sw_conductance(dem)

plot(raster(conductance),
       col = col.palette,
       colNA = "black",
  )