为什么raster :: crop会更改我的RasterLayer的值？

Question

我对raster :: crop函数有一个奇怪的问题。当我裁剪包含唯一ID值（等于单元格编号）的大型栅格时，生成的较小栅格包含重复值，其中每个单元格应该是唯一的。

我有一个大型栅格（r）裁剪到非洲大陆的SpatialPolygonsDataFrame。我想为每个网格单元格生成一个唯一的ID，然后将较大的RasterLayer裁剪为较小的多边形（一个国家/地区） - 这样，我就可以根据其唯一的单元格将较小的栅格计算中的数据合并到较大的栅格中标识。

但是，当我裁剪包含唯一ID的较大栅格时，新栅格图层不包含唯一ID - 而是存在大量重复值。据我所知，似乎某些ID值被相邻的值替换，从而产生类似于例如的模式。 301,301,303,303,306,306,306而不是连续增加1（当然在同一行内）。较大的r仅包含唯一值，因此只有在将栅格裁剪为较小的栅格后才会出现问题。

多边形和栅格中的投影是相同的，我使用的是最新版本的Raster包。

问题似乎只出现在高分辨率栅格上。使用raster :: mask代替用NA替换多边形外部的单元格在某些情况下也会产生重复。

这个问题对我来说完全是个谜 - 我无法找到任何可能的原因。有谁知道裁剪功能是否存在问题以及它如何处理价值？即使我想出另一种方法来做到这一点，我想知道其他栅格是否会出现这个问题，因此，你的数据在crop功能中是如何被破坏的。我希望有人可以帮我弄清问题是什么。

我在下面创建了一个小例子来重现问题。如果您需要更多信息，请告诉我。

pack <- c("sp", "raster", "rgdal", "dplyr")
lapply(pack, require, character.only = TRUE); rm(pack)

r <- raster(africa, resolution = 1/60/2) ## Create empty raster layer based on extent of Africa polygons and resolution 30 arc seconds. 
values(r) <- 1:ncell(r) ## Generate unique cell ID (equal to cell number)
poly <- africa[3,] ## Subset one country

r_id <- crop(r, poly) ## Crop r to poly
## This is the function that seems to be responsible for the unexpected result. It should return a smaller raster containing the same values in the same cells as the larger raster. Therefore each grid cell value in r_id should be unique. 

as.data.frame(r_id) %>% ## This just to show that the resulting raster contains duplicate values where none should exist
  group_by(layer) %>% 
  summarise(n = n()) %>% 
  arrange(desc(n))

# A tibble: 137,270 x 2
      layer     n
      <dbl> <int>
 1 24774556     3
 2 24774560     3
 3 24774564     3
 4 24774568     3
 5 24774572     3
 6 24774576     3
 7 24774580     3
 8 24774584     3
 9 24774588     3
10 24774592     3
# ... with 137,260 more rows. 

Answer 1

这可能是由于使用4字节浮点数将值写入磁盘时，较大数字的数值不精确。

我怀疑这种情况发生在你的情况下，因为r由（临时）文件备份。设置单元格编号后，请查看数据源

values(r) <- 1:ncell(r)
r

我在这里说明一下

library(raster)
a <- 24774556:24774565
r <- raster(ncol=2, nrow=5)
values(r) <- a
x <- writeRaster(r, "test1.grd", overwrite=TRUE)
v <- values(x)
v
#[1] 24774556 24774556 24774558 24774560 24774560 24774560 24774562 24774564 24774564 24774564
table(v)
#v
#24774556 24774558 24774560 24774562 24774564 
#       2        1        3        1        3 

但是，如果将数据类型设置为“FLT8S”或“INT4U”：

y <- writeRaster(r, "test2.grd", datatype="FLT8S", overwrite=TRUE)
values(y)
#[1] 24774556 24774557 24774558 24774559 24774560 24774561 24774562 24774563 24774564 24774565


z <- writeRaster(r, "test3.grd", datatype="INT4U", overwrite=TRUE)
values(z)
# [1] 24774556 24774557 24774558 24774559 24774560 24774561 24774562 24774563 24774564 24774565

在你的情况下，你可以考虑，而不是values(r) <- 1:ncell(r)，做

 r <- init(r, "cell", datatype="FLT8S", filename="africa_id.grd")

或者，一起跳过这一点。您还可以裁剪零件，加工，然后使用merge将裁剪的部件重新组合在一起。