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

[T C Wilkinson]

Dear Chris, dear Jacob,

Many thanks indeed for suggestions and help.

1. The input raster (dem) definitely shows NA in the areas of sea using click().

2. Tried the 0.0 in the ifelse: no luck, same result.

I tried plotting and clicking for the values of the initial `slope` TransitionLayer: here the area of the sea reports as a value of NaN, so perhaps the problem is not in the transition function, but in the later conversion of slope to cost (in sw_conductance).

I tried changing adjacent() to adjacencyFromTransition(), which does indeed result in the sea being plotted as NA (or presumably NaN), but the values of the conductance on land now don’t seem correct (i.e. some areas that I know are flat are shown as relatively lower conductance than areas which are sloping).

Any ideas?

Thank you again,
Toby



sw_slope <- function(dem, dirs=16, symm=FALSE){
  # calculate the altitudinal differences between cells:
  #altDiffFunction <- function(x){x[2] - x[1]}
  # with a little help from Jacob van Etten and Chris English:
  altDiffFunction <- function(x){ifelse(((!is.na(x[2])) & (!is.na(x[1]))), (x[2] - x[1]), 0.0)}
  
  heightDiff <- transition(dem, altDiffFunction, directions=dirs, symm=symm)
  #divide by the distance between cells
  slope <- geoCorrection(heightDiff)
  return(slope)
} 

slope <- sw_slope(dem)
plot(raster(slope))
click(raster(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)
  adj <- adjacencyFromTransition(slope)
  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)
}

plot(raster(sw_conductance(dem))



On 2 June 2016 at 14:29:42, chris english (englishchristophera at gmail.com) wrote:
> Toby,
>  
> This may seem silly, and a thousandth of something is not zero, but perhaps
> this is a kind of preallocation issue and using Jacob's if else with 0.0
> as else might do the trick.
>  
> Fingers crossed?
> Chris
> On Jun 2, 2016 16:07, "T C Wilkinson" wrote:
>  
> > Dear Jacob,
> >
> > Many thanks for fast reply. [Apologies for the typo and potential
> > confusion on gdistance, I am still relative new to R — and thus borrowing
> > camelCase habit from previous coding.]
> >
> > I’ve tried your suggested re-worked function to avoid NA in the
> > transition, but I still get the same result.
> >
> > When plotted, the `plot(raster(conductance))` still colours the sea (i.e.
> > doesn’t show it in the colNA colour) and the resultant accCosts etc include
> > the sea as part of the result surface.
> >
> > An interactive `click(raster(conductance))` reveals the sea to be a
> > uniform value of 0.001436461.
> >
> > Thanks again,
> > Toby
> >
> >
> >
> > On 2 June 2016 at 12:30:43, Jacob van Etten (jacobvanetten at yahoo.com)
> > wrote:
> > > Dear Toby,
> > > TransitionLayers should normally not hold NA values. If you want to set
> > cell connections
> > > to zero conductance, the transitionFunction should give a value zero if
> > it finds an NA.
> > > Your function gives an NA as result, which is not what you want:
> > > ifelse(((!is.na(x[2])) & (!is.na(x[1]))),(x[2] - x[1]),NA)
> > >
> > > You probably want to do:
> > > ifelse(((!is.na(x[2])) & (!is.na(x[1]))), (x[2] - x[1]), 0)
> > >
> > > Also, note that gDistance is not the name of the package.
> > > Best,
> > > Jacob
> >
> > _______________________________________________
> > R-sig-Geo mailing list
> > R-sig-Geo at r-project.org
> > https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>