1

我有一个RasterBrick，其中包含每日积雪数据，其值为 1、2 和 3（1= 下雪，2= 无雪，3= 被云遮挡）。

一天的积雪示例：

在此处输入图像描述

> snowcover
class       : Large RasterBrick 
dimensions  : 26, 26, 2938  (nrow, ncol, nlayers)
resolution  : 231, 232  (x, y)
extent      : 718990, 724996, 5154964, 5160996  (xmin, xmax, ymin, ymax)
crs         : +proj=utm +zone=32 +datum=WGS84 +units=m +no_defs +ellps=WGS84       
              +towgs84=0,0,0

现在我希望插入被云遮挡的像素（但仅限于单个 RasterLayer 中云量少于 90% 的情况，否则应为该图层保留原始值）。

对于空间插值，我想使用数字高程模型（相同的研究区域并且已经具有相同的分辨率）分别为RasterBrick 的每个图层提取上下雪线边界。上方的雪线表示所有无云像素都归类为雪的海拔高度。较低的雪线标识了所有无云像素也无雪的高度。

在此处输入图像描述

> dem
class       : RasterLayer 
resolution  : 231, 232  (x, y)
extent      : 718990.2, 724996.2, 5154964, 5160996 (xmin, xmax, ymin, ymax)
crs         : +proj=utm +zone=32 +datum=WGS84 +units=m +no_defs +ellps=WGS84       
              +towgs84=0,0,0 
values      : 1503, 2135  (min, max)

对于上部雪线，我需要雪覆盖像素的最小高度（值 = 1）。现在，RasterBrick 的 RasterLayer 中高于此最小高程的所有值为 3 的像素都应重新分类为值 1（假设被雪覆盖）。

另一方面，对于较低的雪线，我需要确定无雪像素的最大海拔（值 = 2）。现在，RasterBrick 的 RasterLayer 中高于此最大高程的所有值为 3 的像素都应重新分类为值 2（假设无雪）。

这可能使用R吗？

我试图利用覆盖功能，但我被困在那里。

# For the upper snowline:
overlay <- overlay(snowcover, dem, fun=function(x,y){ x[y>=minValue(y[x == 1])] <- 1; x})
4

1 回答 1

2

这是一些示例数据

library(raster)
dem <- raster(ncol=8, nrow=7, xmn=720145, xmx=721993, ymn=5158211, ymx=5159835, crs='+proj=utm +zone=32 +datum=WGS84')
values(dem) <- ncell(dem):1
snow <- setValues(dem, c(1, 1, rep(1:3, each=18)))
snow[,c(2,5)] <- NA
snow[3] <- 3


plot(snow)
lines(as(dem, 'SpatialPolygons'))
text(dem)

该图显示了雪等级（1、2、3），高程值位于顶部。 雪地图

我们可以使用掩码，但需要处理缺失值。

msnow <- reclassify(snow, cbind(NA, 0))
# mask to get only the snow elevations
x <- mask(dem, msnow, maskvalue=1, inverse=TRUE)

# minimum elevation of the snow-covered cells
minsnow <- minValue(x)
minsnow 
#[1] 37

# snow elevation = 1
snowy <- reclassify(dem, rbind(c(-Inf, minsnow, NA), c(minsnow, Inf, 1)))
newsnow <- cover(snow, snowy)

s <- stack(dem, snow, newsnow)
names(s) <- c("elevation", "old_snow", "new_snow")

在此处输入图像描述

你非常接近，你可以做到

 r <- overlay(dem, snow, fun=function(e, s){ s[e >= minsnow] <- 1; s})

但请注意，这也会覆盖没有雪的高单元格。

在此处输入图像描述

可以这样修复：

r <- overlay(dem, snow, fun=function(e, s){ s[e >= minsnow & is.na(s)] <- 1; s})

要选择具有超过 x% 单元格且值为 3 的层（这里我使用 34% 的阈值）：

threshold = .34
s <- stack(snow, snow+1, snow+2)
f <- freq(snow)
f 
#     value count
#[1,]     1    14
#[2,]     2    13
#[3,]     3    15
#[4,]    NA    14

nas <- f[is.na(f[,1]), 2]

ss <- subs(s, data.frame(from=3, to=1, subsWithNA=TRUE))
cs <- cellStats(ss, sum)
csf <- cs / (ncell(snow) - nas)
csf
#  layer.1   layer.2   layer.3 
#0.3571429 0.3095238 0.3333333 

i <- which(csf < threshold)
use <- s[[i]]
#use
class       : RasterStack 
dimensions  : 7, 8, 56, 2  (nrow, ncol, ncell, nlayers)
resolution  : 231, 232  (x, y)
extent      : 720145, 721993, 5158211, 5159835  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=utm +zone=32 +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0 
names       : layer.2, layer.3 
min values  :       2,       3 
max values  :       4,       5
于 2018-11-06T00:46:18.330 回答