地形阴影功能不使用我自己的数据堆叠不同的角度和方向（如示例）

Question

我尝试使用自己的数据从

shade

中的

terra

函数复制示例的最后一部分，以使用多个方向和角度来改进山体阴影，但生成的对象没有值。因此，

Reduce

函数不会产生均值。具有简单阴影参数（具有一个角度和方向，如示例）的对象效果很好。

没有任何错误消息，只是最终输出栅格中没有值，如下所示。不确定这是否是因为

terra

的大小限制（也许？）。

我的代码：

library(terra)
elevation <- terra::rast(file.choose())


alt<- terra::mask(elevation, bcr)


alt <- terra::disagg(alt, 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)

无论是在

Reduce

之前还是之后，h都有这个输出（并且它似乎在不到一秒的时间内运行，而对于只有1个角度和方向的运行，山则不是这种情况）：

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 

我尝试使用的海拔可以在这里找到： https://drive.google.com/file/d/12yezqqNutfD19eyxJmbL9WqvwTdcjT3e/view?usp=sharing

如有任何意见，我们将不胜感激！

干杯和感谢！！

Answer 1

使用您的数据我得到以下信息：

alt = rast('~/Downloads/elevation_cropped.tif')

alt2 = disagg(alt, 10, method = 'bilinear')
alt2_slope = terrain(alt2, 'slope', unit='radians')
alt2_aspect = terrain(alt2, 'aspect', unit='radians')
hill1 = terra::shade(alt2_slope, alt2_aspect, 40, 270)

hill1                                   
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      : spat_12cdd925ddc96_1232345.tif 
varname     : elevation_cropped 
name        :  hillshade 
min value   : 0.01107622 
max value   : 0.98851299 

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,        ? ,        ? ,        ?  

嗯，到目前为止一切顺利；但是，我们应该如何处理

hill2

中的那些？（s）标记呢？显然，就价值观而言，

hill1

和

hill2

之间还发生了其他事情。我们可以检查一下吗？

hill2[,,2]
   [99999,]       NaN
 [ reached getOption("max.print") -- omitted 296835201 rows ]
Warning message:
[readValues] raster has no values

有点良性 - 我们可以写入文件吗？

hill2 = terra::shade(alt2_slope, alt2_aspect, angle = c(45,45,45,80), direction = c(225,270,315,135), filename = 'hill2_multi_shade.tif')
Error: [shade] there are no cell values

不。我们可以分配 hill2[,,2] 吗？

hill2_2 = hill2[,,2]
hill2_2
[99999,]       NaN
 [ reached getOption("max.print") -- omitted 296835201 rows ]

所以，

Reduce

在hill2的这个过程中非常实用，因为更好的值已经写入第1层。除此之外，我会等待其他人发表评论。否则，如果目标是“查看”角度/方向影响，则这些将作为

hill1

完成，改变角度/方向，因为

hill2

方法将值合并到第一层。没有Reduce调用，并尝试执行以下操作：

plot(hill2[,,1])
Killed
chris@jacie2:~/stackoverflow_r$ 

这不是情节一所期望的良性结果。也许其他人会解释。