I am trying to calculate the slope from a large number of DEM rasters. It seems like Geomorphometry.jl has an implementation of the Horn (1981) method, but are there other library options?
I have some custom code using Stencils.jl, which will be much faster than Geomorphometry. But I never got around to testing it or putting it into Geomorphometry.jl, but we have talked about it, and about using Stencils.jl
using Rasters
using Stencils
using Stencils: Window
abstract type SlopeFilter end
abstract type SlopeConvolution <: SlopeFilter end
# Not all slope algorithms can provide aspect
abstract type SlopeAspectConvolution <: SlopeConvolution end
struct FD2 <: SlopeAspectConvolution end
struct FD3Reciprocal <: SlopeAspectConvolution end
struct FD3ReciprocalSquared <: SlopeAspectConvolution end
struct FD3Linear <: SlopeAspectConvolution end
struct FDFrame <: SlopeAspectConvolution end
struct SimpleDifference <: SlopeConvolution end
@inline function aspect_filter(method::SlopeConvolution, n::Window, d)
fx, fy = _slope_conv(method, n, d)
return _aspect(fx, fy)
end
@inline function slope_filter(method::SlopeConvolution, n::Window, d)
fx, fy = _slope_conv(method, n, d)
return _slope(fx, fy)
end
@inline function slopeaspect_filter(method::SlopeAspectConvolution, n::Window, d)
fx, fy = _slope_conv(method, n, d)
return _slope(fx, fy), _aspect(fx, fy)
end
_slope(fx, fy) = atan(√(fx^2 + fy^2))
function _aspect(fx, fy)
(ismissing(fx) || ismissing(fy)) && return missing
# Rotate - we want high Y (north) as the origin
# TODO: pass through the Order for X/Y dims
# So the result always has zero at North
-atan(fx, fy)
end
@inline function _slope_conv(::FD2, n::Window, d)
fx = (n[6] - n[4]) / 2d
fy = (n[8] - n[2]) / 2d
return fx, fy
end
@inline function _slope_conv(::FD3Reciprocal, n::Window, d)
fx = (n[3] -n[1] + √(2(n[6] - n[4])) + n[9] - n[7]) / (4 + 2 * √(2)) * d
fy = (n[7] -n[1] + √(2(n[8] - n[2])) + n[9] - n[3]) / (4 + 2 * √(2)) * d
return fx, fy
end
@inline function _slope_conv(::FD3Linear, n::Window, d)
fx = (n[3] - n[1] + n[6] - n[4] + n[9] - n[7]) / 6d
fy = (n[7] - n[1] + n[8] - n[2] + n[9] - n[3]) / 6d
return fx, fy
end
@inline function _slope_conv(::FD3ReciprocalSquared, n::Window, d)
fx = (n[3] - n[1] + 2(n[6] - n[4]) + n[9] - n[7]) / 8d
fy = (n[7] - n[1] + 2(n[8] - n[2]) + n[9] - n[3]) / 8d
return fx, fy
end
@inline function _slope_conv(::FDFrame, n::Window, d)
fx = (n[3] - n[1] + n[9] - n[7]) / 4d
fy = (n[7] - n[1] + n[9] - n[3]) / 4d
return fx, fy
end
@inline function _slope_conv(::SimpleDifference, n::Window, d)
fy = (n[5] - n[2]) / d
return fy, fy # No aspect
end
# These are the functions for users to call:
# slope, aspect and slopeaspect with the form `slope(raster, method; cellsize)
for (f, filt) in (:slope => :slope_filter, :aspect => :aspect_filter, :_slopeaspect => :slopeaspect_filter)
@eval begin
function $(f)(elevation::Raster, method=FD2(); cellsize=1)
hood = Window{1}();
padval = missingval(elevation)
newdata = Stencils.mapstencil(hood, parent(elevation); boundary=Stencils.Remove(padval)) do w
$(filt)(method, w, cellsize)
end
rebuild(elevation; data=newdata, name=$(QuoteNode(f)))
end
function $(f)(elevation::AbstractArray, method=FD2(); cellsize=1)
hood = Window{1}()
Stencils.mapstencil(hood, elevation, parent(elevation)) do w, e
$(filt)(method, w, cellsize)
end
end
end
end
function slopeaspect(elevation, method=FD2(); cellsize=1)
sa = _slopeaspect(elevation, method; cellsize)
slope = first.(sa)
aspect = last.(sa)
nt = (; slope, aspect)
return elevation isa Raster ? RasterStack(nt) : nt
end
I could also rewrite this using a NamedStencil to get rid of all the magic numbers
1 Like
There’s also a ZevenbergenThorne implementation in Geomorphometry (I’m one of the authors), but I’m not sure what you’re looking for? Is it performance or more/better algorithms, or anything else?
In terms of speed, the current implementations should be fast enough for most purposes (running in 0.05s for a raster with 1M cells), and there are indeed plans to use Stencils to make it even faster.
There’s a gdaldem executable in GDAL_jll (but not with different algorithms, and certainly not faster). I’m not aware of other native Julia libraries.
One can access gdaldem via GMT.jl or use the GMT’s grdgradient module to compute the slope.
Wow. First, I am blown away by how helpful the Julia community is. (I’m new). Thank you all for such rapid and helpful replies.
While speed is important, I was thinking more about sink/peak/line pre-processing, interpolation of missing data, and memory management for large datasets. It looks like I will need to process 1000 tiles (~2TB) of data, so I am curious about memory management methods for addressing the need to mosaic tiles in order to avoid edge effects.
Additional slope calculation methods are a bonus, but are not a critical concern for me at the moment.
To my knowledge there is currently no well-tested tooling for formulating and end-to-end custom overlapping moving window pipeline on an out-of-core mosaic in the Julia Geo ecosystem.
Just out of curiosity and maybe to make it easier to help, is there any library from another programming language that provides the functionality you are looking for?
I was wondering if this is something DiskArrayEngine could automate in future… with the halo padding required for slope filters?
Would be amazing if we could automate using Stencils.jl with DiskArrayEngine.jl, I guess a bit like ParallelStencil.jl
Yes, this is exactly what I had in mind as well. Seeing your code I think it would be quite easy to wrap into a moving window DiskArrayEngine function so that DiskArrayEngine can take care of memory management and IO and Stencils.jl does the fast moving window loop.
However, to me DiskArrayEngine is currently too experimental yet which kept me from suggesting it as a solution here. But as I said I would be very interested in solutions from other languages to get inspired by there APIs, compare performance etc…
I am not aware of a library in R that does all of these things. I am aware that the R raster/terra library has memory management settings that enable memory-based chunking depending on the specific raster size, I suspect that this makes the boundary problem worse, but I haven’t tried it.
I guess this would be falling into Dask territory: Introduction to Geospatial Raster and Vector Data with Python: Parallel raster computations using Dask