unexpected results with shade when stacking different angles and directions with large file

[francisvolh]

I am having issues with the example code for the shade function in terra when using a larger raster and attempting the stacked example of angles and directions. . It seems the function does run but the produced object has no layers (as expected from the stacked shade angles).

We have been trying to figure out the issue in this stackoverflow questions: https://stackoverflow.com/questions/78116098/terra-shade-function-not-stacking-different-angles-and-directions-with-my-own-da?noredirect=1#comment137741126_78116098

Coming to the conclusion that it is a memory issue that breaks the function internally and no error is reported. Also, when attempting to produce 4 individual shades, and averaging them, it works (with my data).

This is a reproducible example, with a different map from the stackoverflow (so I don't have to paste the link here) but its producing the same results in my computer:

library(giscoR)
library(elevatr)
country_sf <- giscoR::gisco_get_countries(
  country = "PE",
  resolution = "1"
)

elev <- elevatr::get_elev_raster(
  locations = country_sf, #warn about locations
  z = 7, clip = "locations"
)




alt<- terra::mask(terra::rast(elev), terra::vect(country_sf))


alt <- terra::disagg(alt, 10, method="bilinear")

slope <- terra::terrain(alt, "slope", unit="radians")

aspect <- terra::terrain(alt, "aspect", unit="radians")


hill <- terra::shade(slope, aspect, 40, 270)

terra::plot(hill, col=grey(0:100/100), legend=FALSE)

The previous output works fine. But the following just runs h but produces either two things:

1. no layers or no error message
2. An object as shown below (when more memory is available maybe? but still not enough) where there is one layer and the other 3 are "question marked"/empty, but the layers "exist".

#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))

Output with bug:

h
class       : SpatRaster 
dimensions  : 33820, 23310, 4  (nrow, ncol, nlyr)
resolution  : 0.0005423858, 0.0005423858  (x, y)
extent      : -81.32385, -68.68084, -18.35433, -0.01084772  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 +no_defs 
sources     : spat__2.tif  
              memory  
              memory  
              memory  
names       :  hs_45_225, hs_45_270, hs_45_315, hs_80_135 
min values  : -0.3026005,        ? ,        ? ,        ?  
max values  :  0.9999985,        ? ,        ? ,        ?  

A raster with some more resolution, may produce zero layer but an SpatRaster with dimensions, resolution, extent, but no min or max values.

I tried just producing the individual shades, and getting a mean via terra::app() which seemed to work.

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)

meanh <- terra::app(stack1, "mean")

A secondary issue, was that when looking But when I looked into the example in the documentation, using my "optional method", the average products seem different (just checking the min and max values of the mean output rasters), which Chris mentioned in the SO question that should not be the case.

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 version of the object with the documentation example data ("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

packageVersion("terra")
[1] ‘1.7.73’