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
- Reformat data to include 2 minute time intervals
- Modify pop-up message in leaflet map
- Is there a fast/clever way to return a logical vector if elements of a vector are in at least one interval?
- Conditional bslib tooltip
- Set CMU Sans Serif font on pdf quarto
- How to remove string before a certain delimiter in R?
- Leaflet map displays wrong location based on lat,long data
- Concatenate while dropping specific value strings
- R and ggplot with st_crop: map cropping does not work as wanted
- Generating a time series in 10 Hz
- Replace with 0 the first calculation in mutate and rollapply in R
- Change radius by site type geom_scatterpie
- How to stop R from executing when value is non-numeric
- How to solve pagedown_chrome_print failed to generate output due to fail to open favicon.ico
- Attempting to calculate ratios in multiple columns in dataframe via a for loop
- How to count the number of occurences in a column using dplyr
- Convert roman numerals to numbers in R
- Arranging character values of multiple columns that contain the same name with R
- Struggling to understand complete predictive model process in R
- Random change in obj_addr() output when including the objects into a list and vectorizing over them
- Pass vector of unquoted column names to a function and use first element
- retrieve list of urls from Amazon S3 bucket using R
- officer - save document - dont overwrite
- Apply function to list depending on values in another list
- How to return multiple columns from one function inside `mutate` and allow me to name the output columns?
- plotting 3-way interaction with ggeffects, theme changes not being applied to all facets
- Speed up R function that combines linestrings into a single linestring
- Find values in data frame that don't contain specific string and replace them in R
- Error downloading a csv in zip from website with get in R
- How to calculate row sums in R with dynamic columns