通过使用附加光栅（数字高程模型）对 RasterBrick 中的值进行重新分类
Reclassify values in a RasterBrick by the use of an additional Raster (Digital elevation model)
我有一个 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})
这是一些示例数据
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})
到 select 层超过 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
我有一个 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})
这是一些示例数据
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})
到 select 层超过 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