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

Jacob van Etten jacobvanetten at yahoo.com
Thu Jun 2 19:16:40 CEST 2016


There is a contradiction in that you create zero conductance for the sea, but then use this very long formula to convert from resistance (cost) to conductance, whereas you already have conductance values for the sea. The formula gives this result for a value of zero.
# value zerocost <- c(0)adj <- 1
# your complicated formula1 / ((((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)))[1] 0.6097561
It usually works better if you think about the conversion you need and write a single transitionFunction and test this with some values before you plug it into the transition() function. The transition function does nothing special with it, so it is just a matter of clear thinking.
Actually, there is no need to work with adjacency() in your case, I think. Try to create a single transitionFunction...
Jacob
 

    On Thursday, 2 June 2016, 11:00, chris english <englishchristophera at gmail.com> wrote:
 

 Toby,
I've only just installed Jacob's gdistance package and don't know as yet the underlying assumptions nor the, what I will term 'services', his package offers. Were I in your position, having my DEM, and knowing what you already know about using gdistance, I would swap out my sea NA(s) for 0.0(s), in my case probably using GDAL, and see how things went thereafter. 
Jacob will have a much better suggestion, but basically it seems like 'package x doesn't like my data, so I'll make my data more consumable for package x'. Wish I could be of more help but I'm playing in a different distance space than Hausdorf and probably shouldn't comment further, except as to what I'd try to finagle as to my DEM.
Chris 
On Thu, Jun 2, 2016 at 7:28 PM, T C Wilkinson <open at tobywilkinson.co.uk> wrote:

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
>





  
	[[alternative HTML version deleted]]



More information about the R-sig-Geo mailing list