terra shade function not stacking different angles and directions with my own data (as per example)
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:
library(terra)
terra 1.7.73
alt = rast('~/Downloads/elevation_cropped.tif')
alt = rast('~/Downloads/elevation_cropped.tif')
alt_slope = terrain(alt, 'slope', unit='radians', filename = 'alt_slope.tif')
alt_aspect = terrain(alt, 'aspect', unit='radians', filename = 'alt_aspect.tif')
alt_shade = shade(alt_slope, alt_aspect, angle = c(45,45,45,80), direction = c(225,270,315,135), filename = 'alt_shade.tif')
alt_shade
class : SpatRaster
dimensions : 1864, 1593, 4 (nrow, ncol, nlyr)
resolution : 1000, 1000 (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
source : alt_shade.tif
varname : elevation_cropped
names : hs_45_225, hs_45_270, hs_45_315, hs_80_135
min values : 0.1426585, 0.1348058, 0.1915113, 0.6898556
max values : 0.9902307, 0.9877117, 0.9785777, 0.9999999
disagg
alt2 = disagg(alt, 10, method = 'bilinear')
alt2_slope = terrain(alt2, 'slope', unit='radians')
alt2_aspect = terrain(alt2, 'aspect', unit='radians')
hill2 = terra::shade(alt2_slope, alt2_aspect, angle = c(45,45,45,80), direction = c(225,270,315,135))
hill2
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__2.tif
memory
memory
memory
varnames : elevation_cropped
elevation_cropped
elevation_cropped
...
names : hs_45_225, hs_45_270, hs_45_315, hs_80_135
min values : 0.04666001, ? , ? , ?
max values : 0.99999017, ? , ? , ?
So, ? 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
f <- system.file("ex/elev.tif", package="terra")
r <- rast(f)
dim(r)
[1] 90 95 1
and this is very low resolution. Unfortunately, for the unwary (which would be both of us), the disagg process 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 further disagg were 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 wopts list allows processing to go
alt2_shade = shade(alt2_slope, alt2_aspect, angle = c(45,45,45,80), direction = c(225,270,315,135), filename = 'alt2_shade2.tif', overwrite = TRUE, wopt = list(memmax = 0.6, verbose = TRUE))
filename : /tmp/Rtmpt0OHpP/spat__2.tif
compute stats : 1, GDAL: 0, minmax: 0, approx: 1
driver : GTiff
disk available: 782.5 GB
disk needed : 1.1 GB
memory avail. : 0.6 GB
memory allow. : 0.36 GB
memory needed : 8.849 GB (4 copies)
in memory : false
block size : 454 rows
n blocks : 42
pb : 3
filename : /tmp/Rtmpt0OHpP/spat__2.tif
compute stats : 1, GDAL: 0, minmax: 0, approx: 1
driver : GTiff
disk available: 782.5 GB
disk needed : 1.1 GB
memory avail. : 0.6 GB
memory allow. : 0.36 GB
memory needed : 8.849 GB (4 copies)
in memory : false
block size : 454 rows
n blocks : 42
pb : 3
progress bar
filename : /tmp/Rtmpt0OHpP/spat__2.tif
compute stats : 1, GDAL: 0, minmax: 0, approx: 1
driver : GTiff
disk available: 782.5 GB
disk needed : 1.1 GB
memory avail. : 0.6 GB
memory allow. : 0.36 GB
memory needed : 8.849 GB (4 copies)
in memory : false
block size : 454 rows
n blocks : 42
pb : 3
filename : /tmp/Rtmpt0OHpP/spat__2.tif
compute stats : 1, GDAL: 0, minmax: 0, approx: 1
driver : GTiff
disk available: 782.5 GB
disk needed : 1.1 GB
memory avail. : 0.6 GB
memory allow. : 0.36 GB
memory needed : 8.849 GB (4 copies)
in memory : false
block size : 454 rows
n blocks : 42
pb : 3
progress bar
filename : alt2_shade2.tif
compute stats : 1, GDAL: 0, minmax: 0, approx: 1
driver : GTiff
disk needed : 4.4 GB
memory avail. : 0.6 GB
memory allow. : 0.36 GB
memory needed : 17.699 GB (2 copies)
in memory : false
block size : 227 rows
n blocks : 83
pb : 3
alt2_shade
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
source : alt2_shade2.tif
varname : elevation_cropped
names : hs_45_225, hs_45_270, hs_45_315, hs_80_135
min values : 0.634741, 0.634741, 0.634741, 0.634741
max values : 1.000000, 1.000000, 1.000000, 1.000000
But can it usefully be plotted? First Reduce(
alt2_shade_reduce = Reduce(mean, alt2_shade)
alt2_shade_reduce
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
source : alt2_shade2.tif
varname : elevation_cropped
name : hs_45_225
min value : 0.634741
max value : 1.000000
writeRaster(alt2_shade_reduce, 'alt2_shade_reduce_.tif')
Plot looks fine, and better.
ls -lah stackoverflow_shade
-rw-rw-r-- 1 chris chris 899M Mar 10 16:54 alt2.tif
-rw-rw-r-- 1 chris chris 942M Mar 11 10:04 alt2_slope.tif
-rw-rw-r-- 1 chris chris 931M Mar 11 10:05 alt2_aspect.tif
-rw-rw-r-- 1 chris chris 766M Mar 11 11:51 alt2_shade_reduce_.tif
- R package gtsave cropping image, despite setting vwidth
- Reason for unexpected output in subsetting data frame - R
- Fitting a custom model in R
- Difference between 3D plots using fitted and predicted values
- R lookup interpolate between values or closest match
- How to remove duplicate legend in R ggplot
- How can I change the number of breaks of x axis for each graph individually for each facet in a facet_grid in R?
- Get the specific index of occurence when working with period.apply
- Identify single outliers using `diff()`
- Why does `ifelse()` behave differently than `if (cond) {...} else {...}`?
- Getting infinite values when expecting integers after mutating
- FUN.VALUE missing when using vapply
- how to add tables to a forest plot
- Adding leading zeros in the character and numeric mixed column names in R
- I need to update a package
- R Dplyr: Summarizing a column, if it is present
- Cross-validation with nb method
- Reverse order in R leaflet continuous legend
- How to assign different colors to points based on values in a forest plot in R?
- How to visualize interaction using emmeans
- Adding 1 to the last and second last numbers in my data
- Horizontal f7Radio in shinyMobile
- How to use a template in vim
- Regex pattern to include string while excuding other patterns
- R >4.1 syntax: Error: function 'function' not supported in RHS call of a pipe
- Mutating to return the average of a column in each row
- R View() function does not show data frame
- Data manipulation in R in LINQ style
- Unexpected output of dplyr::top_n
- Get position index of second highest value by row for selected columns