I have tried to replicate the las t section of the example from the shade function in terra with my own data to improve the hillshade with the use of several directions and angles, but the resulting object has no values. And therefore, no mean is produced with the Reduce function. The object with the simple shade parameter (with one angle and direction, as per the example) works fine.
There is no error message whatsoever, just no values in the final output raster as show below. not sure if this is because of a size limitation of terra (maybe?).
My code:
library(terra)
elevation <- terra::rast(file.choose())
alt <- terra::disagg(elevation , 10, method="bilinear")
slope <- terra::terrain(alt, "slope", unit="radians")
aspect <- terra::terrain(alt, "aspect", unit="radians")
#basic hillshade
hill <- terra::shade(slope, aspect, 40, 270)
terra::plot(hill)
terra::plot(hill, col=grey(0:100/100), legend=FALSE, mar=c(2,2,1,4))
#get a better shade with different angles and directions
h <- terra::shade(slope, aspect, angle = c(45, 45, 45, 80), direction = c(225, 270, 315, 135))
h <- Reduce(mean, h)
terra::plot(h)
Either before the Reduce or after, h has this output (and its seems to run in less than a second, which is not the case for run with only 1 angle and direction, hill):
h
class : SpatRaster
dimensions : 18640, 15930, 1 (nrow, ncol, nlyr)
resolution : 100, 100 (x, y)
extent : -2791000, -1198000, 1522000, 3386000 (xmin, xmax, ymin, ymax)
coord. ref. : +proj=laea +lat_0=45 +lon_0=-100 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs
The elevation I am trying to use could be found here: https://drive.google.com/file/d/12yezqqNutfD19eyxJmbL9WqvwTdcjT3e/view?usp=sharing
Any input would be appreciated!
Cheers and Thanks!!
EDIT: After Chris pointed out that it may be a memory issue of the function attempting to run and stack 4 shade at a time to improve the shading I attempted something that may be a "fix", or just shows that the function is breaking internally with large rasters when trying to produce several shades.
I just tried this:
h1 <- terra::shade(slope, aspect, angle = 45, direction = 225)
h2 <- terra::shade(slope, aspect, angle = 45, direction = 270)
h3 <- terra::shade(slope, aspect, angle = 45, direction = 315)
h4 <- terra::shade(slope, aspect, angle = 80, direction = 135)
stack1 <- c(h1, h2, h3, h4)
#check the stacked shades
stack1
class : SpatRaster
dimensions : 18640, 15930, 4 (nrow, ncol, nlyr)
resolution : 100, 100 (x, y)
extent : -2791000, -1198000, 1522000, 3386000 (xmin, xmax, ymin, ymax)
coord. ref. : +proj=laea +lat_0=45 +lon_0=-100 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs
sources : spat_5bc87444484_23496.tif
spat_5bc8105f55ed_23496.tif
spat_5bc82d144698_23496.tif
spat_5bc819d7aac_23496.tif
varnames : elevation_cropped (1)
elevation_cropped (1)
elevation_cropped (1)
...
names : hillshade, hillshade, hillshade, hillshade
min values : 0.04666001, 0.09771556, 0.04562845, 0.634741
max values : 0.99999017, 0.99748391, 0.99444449, 1.000000
meanh <- terra::app(stack1, "mean")
par(mfrow= c(1,2))
terra::plot(hill, col=grey(0:100/100), legend=FALSE, mar=c(2,2,1,4))
terra::plot(meanh, col=grey(0:100/100), legend=FALSE, mar=c(2,2,1,4))
The mean shade is produced, and its surely different than the original shade. Memory didn't explode. So this is maybe a bug or issue within the function.
Edit2: I compared this 1 by 1 approach, and the "group shading" from the example in the documentation, just to check if the raster results were the same, and there is a slight difference. The max and min value don't match. But they are visually similar. So there is something else the function is doing internally. I haven't been able to look at source code. Not sure I want to dig into that.
Example from the documentation
Original h object from the documentation example (grouped shading)
h
class : SpatRaster
dimensions : 900, 950, 1 (nrow, ncol, nlyr)
resolution : 0.0008333333, 0.0008333333 (x, y)
extent : 5.741667, 6.533333, 49.44167, 50.19167 (xmin, xmax, ymin, ymax)
coord. ref. : lon/lat WGS 84 (EPSG:4326)
source(s) : memory
varname : elev
name : hs_45_225
min value : 0.5087930
max value : 0.8495533
My object ("individual shading and mean after")
meanh
class : SpatRaster
dimensions : 900, 950, 1 (nrow, ncol, nlyr)
resolution : 0.0008333333, 0.0008333333 (x, y)
extent : 5.741667, 6.533333, 49.44167, 50.19167 (xmin, xmax, ymin, ymax)
coord. ref. : lon/lat WGS 84 (EPSG:4326)
source(s) : memory
name : mean
min value : 0.6748505
max value : 0.8538219
Using your data I get the following, one just using the resolution of the cropped raster, the next disagg-d:
disaggSo, ? marks or no layers at all, or memory not mapped crash. Multidirectional shade points out that it provides 'better' results for low resolution rasters, and the example file used was
and this is very low resolution. Unfortunately, for the unwary (which would be both of us), the
disaggprocess made it into the ?shade example as a 'better' final product, and we didn't pause to think that perhaps 'elevation_cropped' was already at a good resolution. The ensuing results of furtherdisaggwere certainly unexpected... I still think this is worthy of mention at terra issues at least as to help.As it turns out, it can be done. Using wopt and insights from mem_info, memmax that highlights the memory required but writes to disk, so giving more memory 'elbow room' via
woptslist allows processing to goprogress bar
progress bar
But can it usefully be plotted? First Reduce(
Plot looks fine, and better.