ApproxNA 在大型栅格堆栈中产生不同且不正确的结果
ApproxNA producing different and incorrect results in stacks of large rasters
我正在尝试确认 R Raster 包的最新版本 (2.5-8) 中函数 approxNA 是否存在错误或至少无法解释的行为变化。 稍后注意:正确答案被视为下面答案字符串中的注释,它建议现已发布的软件包更新
我在下面有一个可重现的小例子。
最初的症状是,在长达一年的每日 eMODIS 程序块上,approxNA 命令的大部分输出似乎与输入相同——它有很多很多没有被插值的 NA 值……但是当我切换回旧版本的包 (2.3-12) 和旧版本的 R (3.3.1, Bug-in-your-hair) 它 returns 正确的插值结果。我已经尝试将栅格数据类型设置为 INT2S 并且没有设置特定的栅格数据类型。
我检查了文档,以防有某种理由认为结果会有所不同,但在文档中,可以正常工作的旧软件包版本和不能正常工作的新软件包版本之间没有任何变化。
在我的可重现示例中,错误较少,但仅从绘图视图(包含在示例中)中就可以很明显地看到它们。在这些图中,您可以看到几个图层的北部保留了本应插值的 NA 值。
我是 运行 64 位 Windows 7 和 64 位 R 3.3.2,带有光栅 2.5-8 和 R-Studio 接口,但我也在 R 上进行了测试-控制台具有相同的结果,在具有相同配置的不同计算机上,并且还在 32 位 R 中进行了测试,结果相同。
触发器似乎是第一层或最后一层中存在 NA 值,但错误发生在第一层或最后一层中没有 NA 值的单元格中:例如左上角的单元格仅缺失在原始堆栈的第 2 层中,因此对于从第 1 层到第 3 层的那些单元格应该是一个非常简单的插值。
几个月前我确实按照 R 的说明提交了错误报告,但我不确定该联系信息是否正确,或者它是否已进入垃圾邮件文件夹。
我也很乐意看到建议解决方法方向的答案,只要它不会太慢并且使用通用包或定期维护的包。
我似乎无法正确标记它们,但下面是
示例输入积木图
好输出砖的情节
坏输出砖图
感谢您的任何建议。
##### start Reproducible example
library(raster)
# make a large raster, it must be quite big to have a problem
r <- raster(ncols=8142, nrows=3285)
# make six rasters to stack, use different ranges to better see the interpolation results
r1 <- setValues(r, round(10 * runif(ncell(r)),0))
r2 <- setValues(r, NA)
r3 <- setValues(r, round(100 * runif(ncell(r)),0))
r4 <- setValues(r, round(1000 * runif(ncell(r)),0))
r5 <- setValues(r, round(50 * runif(ncell(r)),0))
r6 <- setValues(r, round(100 * runif(ncell(r)),0))
# insert some NA values
r1[100:600,] <- NA
r3[,1500:1950] <- NA
r5[,400:800] <- NA
r6[2500:2900,] <- NA
# insert some constant values to make it easier to see if interpolation is working
r4[300:500,] <- 750
r6[300:400,] <- 20
r6[400:500,] <- 100
# make a stack
s <- stack(r1,r2,r3,r4,r5,r6)
# see what the pre-interpolation stack looks like
#plot(s)
plot(s, col = colorRampPalette(c("royalblue","springgreen","yellow","red"))(255))
# interpolate, this takes a while
x2 <- approxNA(s, method = "linear", rule=2)
# see what the post interpolation looks like, there should be filled in values for all
# parts of all layers but instead it’s retaining many of the NA values
#plot(x2)
plot(x2, col = colorRampPalette(c("royalblue","springgreen","yellow","red"))(255))
## End reproducible example
我无法重现问题。我稍微简化了例子
library(raster)
r <- raster(ncols=81, nrows=32)
r1 <- setValues(r, 1)
r2 <- setValues(r, NA)
r3 <- setValues(r, 3)
r4 <- setValues(r, 4)
r5 <- setValues(r, 5)
r6 <- setValues(r, 6)
# insert some NA values
r1[1:6,] <- NA
r3[,15:19] <- NA
r5[,4:8] <- NA
r6[25:29,] <- NA
# make a stack
s <- stack(r1,r2,r3,r4,r5,r6)
# see what the pre-interpolation stack looks like
spplot(s)
# in memory
x1 <- approxNA(s, method = "linear", rule=2)
spplot(x1)
# to disk (that is, simulate a large raster
rasterOptions(todisk=TRUE)
x2 <- approxNA(s, method = "linear", rule=2)
rasterOptions(todisk=FALSE)
spplot(x2)
我看不出有什么不同。 x1
和 x2
看起来相同
x2 - x1
#class : RasterBrick
#dimensions : 32, 81, 2592, 6 (nrow, ncol, ncell, nlayers)
#resolution : 4.444444, 5.625 (x, y)
#extent : -180, 180, -90, 90 (xmin, xmax, ymin, ymax)
#coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
#data source : in memory01_223851_3440_81523.grd
#names : layer.1, layer.2, layer.3, layer.4, layer.5, layer.6
#min values : 0, 0, 0, 0, 0, 0
#max values : 0, 0, 0, 0, 0, 0
我正在尝试确认 R Raster 包的最新版本 (2.5-8) 中函数 approxNA 是否存在错误或至少无法解释的行为变化。 稍后注意:正确答案被视为下面答案字符串中的注释,它建议现已发布的软件包更新
我在下面有一个可重现的小例子。 最初的症状是,在长达一年的每日 eMODIS 程序块上,approxNA 命令的大部分输出似乎与输入相同——它有很多很多没有被插值的 NA 值……但是当我切换回旧版本的包 (2.3-12) 和旧版本的 R (3.3.1, Bug-in-your-hair) 它 returns 正确的插值结果。我已经尝试将栅格数据类型设置为 INT2S 并且没有设置特定的栅格数据类型。
我检查了文档,以防有某种理由认为结果会有所不同,但在文档中,可以正常工作的旧软件包版本和不能正常工作的新软件包版本之间没有任何变化。
在我的可重现示例中,错误较少,但仅从绘图视图(包含在示例中)中就可以很明显地看到它们。在这些图中,您可以看到几个图层的北部保留了本应插值的 NA 值。
我是 运行 64 位 Windows 7 和 64 位 R 3.3.2,带有光栅 2.5-8 和 R-Studio 接口,但我也在 R 上进行了测试-控制台具有相同的结果,在具有相同配置的不同计算机上,并且还在 32 位 R 中进行了测试,结果相同。
触发器似乎是第一层或最后一层中存在 NA 值,但错误发生在第一层或最后一层中没有 NA 值的单元格中:例如左上角的单元格仅缺失在原始堆栈的第 2 层中,因此对于从第 1 层到第 3 层的那些单元格应该是一个非常简单的插值。
几个月前我确实按照 R 的说明提交了错误报告,但我不确定该联系信息是否正确,或者它是否已进入垃圾邮件文件夹。
我也很乐意看到建议解决方法方向的答案,只要它不会太慢并且使用通用包或定期维护的包。
我似乎无法正确标记它们,但下面是
示例输入积木图
好输出砖的情节
坏输出砖图
感谢您的任何建议。
##### start Reproducible example
library(raster)
# make a large raster, it must be quite big to have a problem
r <- raster(ncols=8142, nrows=3285)
# make six rasters to stack, use different ranges to better see the interpolation results
r1 <- setValues(r, round(10 * runif(ncell(r)),0))
r2 <- setValues(r, NA)
r3 <- setValues(r, round(100 * runif(ncell(r)),0))
r4 <- setValues(r, round(1000 * runif(ncell(r)),0))
r5 <- setValues(r, round(50 * runif(ncell(r)),0))
r6 <- setValues(r, round(100 * runif(ncell(r)),0))
# insert some NA values
r1[100:600,] <- NA
r3[,1500:1950] <- NA
r5[,400:800] <- NA
r6[2500:2900,] <- NA
# insert some constant values to make it easier to see if interpolation is working
r4[300:500,] <- 750
r6[300:400,] <- 20
r6[400:500,] <- 100
# make a stack
s <- stack(r1,r2,r3,r4,r5,r6)
# see what the pre-interpolation stack looks like
#plot(s)
plot(s, col = colorRampPalette(c("royalblue","springgreen","yellow","red"))(255))
# interpolate, this takes a while
x2 <- approxNA(s, method = "linear", rule=2)
# see what the post interpolation looks like, there should be filled in values for all
# parts of all layers but instead it’s retaining many of the NA values
#plot(x2)
plot(x2, col = colorRampPalette(c("royalblue","springgreen","yellow","red"))(255))
## End reproducible example
我无法重现问题。我稍微简化了例子
library(raster)
r <- raster(ncols=81, nrows=32)
r1 <- setValues(r, 1)
r2 <- setValues(r, NA)
r3 <- setValues(r, 3)
r4 <- setValues(r, 4)
r5 <- setValues(r, 5)
r6 <- setValues(r, 6)
# insert some NA values
r1[1:6,] <- NA
r3[,15:19] <- NA
r5[,4:8] <- NA
r6[25:29,] <- NA
# make a stack
s <- stack(r1,r2,r3,r4,r5,r6)
# see what the pre-interpolation stack looks like
spplot(s)
# in memory
x1 <- approxNA(s, method = "linear", rule=2)
spplot(x1)
# to disk (that is, simulate a large raster
rasterOptions(todisk=TRUE)
x2 <- approxNA(s, method = "linear", rule=2)
rasterOptions(todisk=FALSE)
spplot(x2)
我看不出有什么不同。 x1
和 x2
看起来相同
x2 - x1
#class : RasterBrick
#dimensions : 32, 81, 2592, 6 (nrow, ncol, ncell, nlayers)
#resolution : 4.444444, 5.625 (x, y)
#extent : -180, 180, -90, 90 (xmin, xmax, ymin, ymax)
#coord. ref. : +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0
#data source : in memory01_223851_3440_81523.grd
#names : layer.1, layer.2, layer.3, layer.4, layer.5, layer.6
#min values : 0, 0, 0, 0, 0, 0
#max values : 0, 0, 0, 0, 0, 0