Accumulation Cost or Cost distance feature

[pierrethiriet]

Thanks for this very valuable package. It is not a bug report but more of a feature question. Would it be feasible, or is there any interest in a cost distance analysis solution such as:

• Accumulated Cost from SAGA
• Cost distance from R Terra
  Thanks a lot

[rafaqz]

I have a fast algorithm for this already, but it relies in my not yet registered Neighborhoods.jl package.

Another package Neighborhood.jl was registered so I cant get the close name. I havent thought of another one yet!

When thats finaliaed I can add cost distance, slope, aspect, focal etc methods here.

[pierrethiriet]

Thanks for that fast answer!
Such an addition would be fantastic; let's hope the name collision will not be a definitive blocking issue...

[rafaqz]

It may be quire a few months until I get around to adding those functions here ether way.

[pierrethiriet]

Anyway, thanks for your work!

[rafaqz]

Neighborhoods.jl is registered now as Stencils.jl:
https://github.com/rafaqz/Stencils.jl

It should be fairly easy to add these things now, I just need a good test suite for it.

[pierrethiriet]

Thanks for your work !
I'll try Stencils ASAP.

[rafaqz]

That new PR has cost_distance in it, will get it tested and working in the next weeks.

Can I ping you for feedback on The PR at some point? It would be good to make sure the functionality is broad enough

[pierrethiriet]

Thanks for the proposition. I am neither a specialist on the cost_distance questions nor Julia's good programmer. But I surely can try to give feedback as I am also interested in the slope and focal function (I still have not check yet if your implementation allows using several input rasters and arbitrary focal shapes)

[rafaqz]

Im no specialist either, just a random ecologist.

But those kind of questions are exactly what Im looking for!

focal shapes are extremely arbitrary you can defone any offset patterns you like with a Positional stencil.

[rafaqz]

One question I have is what kind of cost functuons should be available out of the box? Ive got a walking time model there at the moment but there could be more and they could be more flexible.

[pierrethiriet]

Our use case might be walking distance / least path cost but with custom weighted cost to model some preferences based on path surroundings (vectorial and raster approaches are still on the table in this particular situation at this stage ). So I have not yet an opinion on how to design and the degree of flexibility of such an option.

[rafaqz]

Ok. Currently Im allowimg users to pass in a function:

f(elev1, elev2, distance, [resistance])

Where resistance is optional if you dont use a resistance matrix.

So you can use any function of those arguments and the result will be minimised over the grid.

Would that be flexible enough?

[pierrethiriet]

Sorry for the delay. I am not sure to understand the proposed function (or the args) properly.
In the tool I did use in the past, either you provide a point location or value cells as starting points only, and you get a "proximity" raster, or you provide an endpoint, and you get a cumulated distance.
In that perspective, I don't fully get the meaning of elev1 and elev2 args. Similarly, Does distance serve as a threshold?
But I may also misunderstand the purpose of the proposed function.

[rafaqz]

Sorry maybe that was confusing. This is the format for a custom cost function you can pass in, with walking_time being the default

Rasters.jl/ext/RastersStencilsExt/cost_distance.jl

Lines 120 to 131 in 43a34df

	"""
	walking_time(e1, e2, d, v)
	Calculate the cost of moving between elevation `e1` to elevation `e2`
	over distance `d` through resistance `r`
	cost = S * I * A * E
	"""
	function walking_time(e1, e2, d, r)
	s = walking_speed(ImhofTobler(), abs(e1 - e2), d)
	return d / (s * r)
	end

The idea is you can write whatever formula you like to define cost. You could have driving time here instead, using the same function arguments. Or energy expended by a laden camel ;)

Just something to return a cost based on elevations, distance, and optionally landscape resisance.

The main cost distanct fuction is here:

Rasters.jl/ext/RastersStencilsExt/cost_distance.jl

Lines 34 to 39 in 43a34df

	cost_distance(origins::AbstractRaster, elevation::AbstractRaster, resistance=nothing; kw...)
	cost_distance(time_taken, origins, elevation, resistance; kw...)
	function cost_distance(
	f::Function, origins::AbstractRaster, elevation::AbstractRaster, resistance=nothing;
	kw...
	)

It has pretty much what you describe - you provide origin points (currently an array but but could be geometries) and an elevation raster, and optionally a resistance surface, and it returns the cumulative cost over the whole area.

The function f in the second method here is where you would pass your own custom functions like driving_time or camel_energy_expenditure, following the format of walking_time

You didn't mention elevation or resistance layers - do you not use those?

[pierrethiriet]

Ok, I'll get a better picture now.
In my use case (if raster approaches remain...), we will definitively use a resistance layer to model the walking path conditions (question of home to waste bin distances), including infrastructure but also social perceptions of the path itself and surroundings. The elevation (and slope) is obviously and usual important parameter, but it may not play a crucial role here. I was thinking of combining the slope to other parameters (hedonist quality of a portion of a path, etc...) in a preprocessed cost (or resistance) raster.

[rafaqz]

Ah ok sounds like an interesting social study.

I hadn't thought about cases where elevation doesn't matter. I'll think about how to make that possible. It would be annoying for you to have to pass in an array of zeros for elevation and have two dummy elevation variables you ignore.

Maybe elevation and resistance should be generalised to behave in the same way, and each should get an argument that is a tuple of the source and destination cell values.

So if you pass in elevation and resistance rasters you would be passed these values in your cost function:

 walking_time(distance, (elev1, elev2), (res1, res2)) 

Or it could be just resistance

 walking_time(distance, (res1, res2)) 

Or whatever you like:

 walking_time(distance, (res1, res2), (somethingelse1, somethingelse1)) 

It could also be done with NamedTuple, like

 walking_time(distance, (elevation=(elev1, elev2), resistance=(res1, res2))) 

And then you could pass in RasterStack with all the layers you need in it, and just index them by name in your algorithm.

[rafaqz]

Also to clarify - you may notice that I don't like how cost_distance works in terra, and this is quite different!

For most of my uses resistance surface isn't the right way to work with elevation, I would like to be able to do path distance with it, more like: https://desktop.arcgis.com/en/arcmap/latest/tools/spatial-analyst-toolbox/path-distance.htm.

So this formulation may be a little more complicated but much more flexible. We can add various prebaked cost functions over time to fill out the functionality - but people should be able to write their own cost equations easily.

That's also possible here because user defined functions will be as fast as C (unless you really mess it up) so you can do whatever you like, its not really doable in R.

[pierrethiriet]

Thanks for your explanation.
A solution where elevation is just a resistance component would be quite nice. It provides a general approach without too much downside.
But, it is pretty strong "but", now I read the description of the path-distance tool from ArcMap, I do realize that vertical cost is different from having a cost based on slope.
Ultimately, I like the idea of a NamedTuple, which clarifies how the function works.

[rafaqz]

Great! Thanks for riffing on that I think we have something better than my current design.

Yes there are a lot of ways to define the cost of slope. Some functions just use the slope in a pixel as a generic cost like you say, but this ignores that paths going on contour are flat, with no vertical cost at all even if tge slope is steep. And a lot of roads and tracks are like that.

[pierrethiriet]

Thanks for the fast development. I had a quick look at the last version (Stencil 0.2.0), and I have just a few questions/remarks:

• Would it be possible for the mapstencil to return a Raster instead of a Matrix if used with Raster ? In most of the cases (if all ?), the result should keep the resolution and the geolocation of the initial raster, no ? Or is there an easy way to geolocate the results based on the source raster?
• I am not sure from the documentation that it is possible to adapt the stencil to the individual raster/Array values. In another way, would it be possible to get StencilArray(Array, x-> Window(x)) where x is the value of the central cell? In an even more general case, would it be feasible to adapt the focal shape according to the values of the different arrays of the same size? We did for example, retrieve the population concerned by odors emission. So we did a moving window shaped by both wind direction and emission plume size (distance...) to sum the population concerned:
  Here was the (..badly...) coded functions:
  First the fun to create a window or circular mask based on a distance and angle value

# Circular moving window function
## Function creating binary matrix
## with angle for windows function
function circle_mat_angle(r::Int, angle)
  focal = zeros(Int, 2r + 1, 2r + 1)
  if 337.5 <= angle <= 360 ||  0 <= angle < 22.5 # North
    [focal[i+r+1,j+r+1] = (i^2+j^2) <= r^2 for i ∈ -r:r, j ∈ -r:r if i >= 0 && abs(j) <= abs(i)]
  elseif 22.5 <= angle < 67.5 # North East
    [focal[i+r+1,j+r+1] = (i^2+j^2) <= r^2 for i ∈ -r:r, j ∈ -r:r if i >= 0 && j <= 0]
  elseif 67.5 <= angle < 112.5 # East
    [focal[i+r+1,j+r+1] = (i^2+j^2) <= r^2 for i ∈ -r:r, j ∈ -r:r if j <= 0 && abs(i) <= abs(j)]
  elseif 112.5 <= angle < 157.5 # South East
    [focal[i+r+1,j+r+1] = (i^2+j^2) <= r^2 for i ∈ -r:r, j ∈ -r:r if i <= 0 && j <= 0]
  elseif 157.5 <= angle < 202.5 # South
    [focal[i+r+1,j+r+1] = (i^2+j^2) <= r^2 for i ∈ -r:r, j ∈ -r:r if i <= 0 && abs(j) <= abs(i)]
  elseif 202.5 <= angle < 247.5 # South West
    [focal[i+r+1,j+r+1] = (i^2+j^2) <= r^2 for i ∈ -r:r, j ∈ -r:r if i <= 0 &&  j >= 0]
  elseif 247.5 <= angle < 292.5 # West
    [focal[i+r+1,j+r+1] = (i^2+j^2) <= r^2 for i ∈ -r:r, j ∈ -r:r if j >= 0 && abs(i) <= abs(j)]
  elseif  292.5 <= angle < 337.5 # North West
    [focal[i+r+1,j+r+1] = (i^2+j^2) <= r^2 for i ∈ -r:r, j ∈ -r:r if i >= 0 && j >= 0]
  end
  return focal
end

And apply the window function to each cell with 1) the array (e.g. population raster), 2) the radius = the plume size and 3) the angle given by the wind direction.

# The moving window function
function moving_window_dist_angle(array, radius, angle, fun::Function)
  nrow, ncol = size(array)
  result = similar(array)  # Empty array
  for j in 1:ncol, i in 1:nrow
      if ismissing(array[i,j]) || ismissing(radius[i,j]) || ismissing(angle[i,j])
        result[i,j] = missing
      else
          r = round(Int, radius[i,j]/1000) # Convert radius distance to km
          mask = circle_mat_angle(r, angle[i,j])
          # Calculate the window matrix boundaries
          imin = i<=r ? 1 : i-r
          imax = i+r > nrow ? nrow : i+r
          jmin = j<=r ? 1 : j-r
          jmax = j+r > ncol ? ncol : j+r
          # Reshape windows size according to the boundaries
          window = array[imin:imax,jmin:jmax]
          mask_sub = mask[(imin-i+r+1):(imax-i+r+1), (jmin-j+r+1):(jmax-j+r+1)]
          window_mask = window .* mask_sub
          result[i,j] = fun(skipmissing(window_mask)) # check the function
      end
  end
  return result
end

The function would be applied as the following:

moving_window_dist_angle(population, emission_distance, wind_direction, sum);

where the population is a raster with estimated population size, emission distance gives the distance of a plume emission in each cell (depending on the local climate) and the local wind direction.
Such an approach would be feasible with Stencil.

• I guess you are using some slope calculation in the future cost_distance function. Would it be possible to expose it publically? But I guess it would fit better in a more formal topo/hydro function package...

Thanks in advance, and maybe I should move those questions to the Stencil package...

[rafaqz]

As for returning Raster, we just need to override a few methods so that Raster rewraps AbstractStencilArray and forwards its methods. We can do this in the RastersStencilsExt extension. So it will plot as a raster and keep spatial information but also work with Stencils.jl.

I've thought about localized window shape a bit with DynamicGrids.jl but never actually implemented it, specifically for this wind direction problem too.

To do it, you can use specific stencils outside of mapwindow by just doing stencil(localstencil, array, I) where I is a cartesian index inside your own broadcast or KernelAbstractions.jl kernel.

So in each of your branches above you would just specify a different window shape. You will need to be careful to make all of this type stable and compile away - actually a big if else block like you use above is the way to do that.

But it could be much faster than your approach above by ensuring that r is known at compile time, then the kernels are unrolled instructions rather than small loops of unknown length. You may need to force this by passing r as Val{r}() and using the type paramter ::Val{R}) where R to build window shapes, like window = Moore{R,2}() .

Building kernels with a runtime radius will be very slow unless you specify the radius in a branch.

Maybe we can build a grouped stencil (like Layered) to do this kind of operation and write out the if/else block for all the options in a generated function. You would put your 8 values or value ranges in the type of the stencil, like Positional puts the offsets in the type. This would be matched with your 8 stencils in if/else block in the generated function.

For slopes, there is a big set of slope algorithms in the PR with cost_distance. See:
https://github.com/rafaqz/Rasters.jl/blob/43a34df63ecda59d510c881017047ed5c14c836f/ext/RastersStencilsExt/slope_aspect.jl

(I'm going to refactor it all to use VonNeumann windows but I think it may make no difference as the compiler can see that it doesn't have to read the corner values)

The cost distance model doesn't actually calculate slope directly, it just passes the cost function the start and end elevations that it can use to calculate cost however it needs to.

[pierrethiriet]

Thanks for all this information. I think I get your points regarding the localized focal approach. I'll check when I can (I would like to recode a previous work with your new tools to see what I can get). But I realize that my coding skills are still lacking...

As for the slope, I indeed saw a paper on the slopes algorithm. But in a broader perspective, would it be relevant to also add all the "standard" DEM calculations (aspect, ruggedness, etc.) here or in another package ?

[rafaqz]

No worries. aspect is already in that PR. And ruggedness is just some statistic of slope? We can totally add that too, feel free to PR.

[rafaqz]

As for the coding skills here I don't think you are far off from the examples above - the key thing is realizing why StaticArrays.jl is the fast way to do this. Having a fixed length and offsets known at compile time means the compiler knows everything about the stencil and is forced to compile specific code for all variants. In your code above the compiler does not know the value of r so has to do everything at runtime. r could be 1000 for all it knows - and in that case your method will be faster! But when r is small, static arrays are faster because the compiler can unroll all the instructions, skip instructions that are not required, and generally optimise things with more information (my knowledge pretty much ends there ;))

Otherwise Stencils.jl isn't doing anything fancy - there are no optimizations at all except generating a StaticArray for everything and running functions of them in a generic KernelAbstractions.jl kernel (thats just so it runs paralel on GPU/CPU without me writing any custom code for it).

[rafaqz]

This functionality should all go into Geomorphometry.jl