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’

[chris-english]

Can be achieved increasing memory in wopt=list(memmax=.06, or .07. See SO where worked out and likely close. Sorry for my bad advice, but workable now, knowing more.

[francisvolh]

Thanks @chris-english !

So modifying the wopt does help but only if filename is given. I attempted to run it without writing the file (without filename) and the code for shade did run and produced an output, but with no layers at all.

I think terra does give out warning or error messages for other process where it runs out of memory (it has happened to me already), so maybe they can add this to the package in a future release? Plus a warning on the documentation that dissag is needed only for low resolution rasters, and using it in larger rasters will produce unexpected results.

I will leave it to the developer team to close the issue if they see fit or whenever they address some of these suggestions!

cheers!

[rhijmans]

Sorry for the long response time. I have updated the manual (warning about disagg) and I have made the method more memory safe.