为啥 raster::crop 会改变我的 RasterLayer 的值？

Posted 2023-02-25

【中文标题】为啥 raster::crop 会改变我的 RasterLayer 的值？

【英文标题】：Why does raster::crop change values of my RasterLayer?为什么 raster::crop 会改变我的 RasterLayer 的值？

【发布时间】：2019-02-12 14:00:04

【问题描述】：

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

我有一个大栅格 (r) 被裁剪到非洲大陆的 SpatialPolygonsDataFrame。我想为每个网格单元生成一个唯一的 ID，然后将较大的 RasterLayer 裁剪为较小的多边形（一个国家）——这样，我将能够随后根据其唯一单元将来自较小栅格计算的数据合并到较大的栅格身份证。

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

多边形和光栅的投影是一样的，我使用的是最新版本的光栅包。

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

这个问题对我来说完全是个谜——我找不到任何可能的原因。有谁知道裁剪功能是否存在问题以及它如何处理值？即使我想出了另一种方法来做到这一点，我想知道其他栅格是否会出现此问题，因此，您的数据在裁剪功能中是否已损坏。我希望有人能帮我找出问题所在。

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

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. 

【问题讨论】：

抱歉，包现在包含在代码中。它们是“sp”、“rgdal”、“raster”和“dplyr”。

我收到Error in raster(africa, resolution = 1/60/2) : object 'africa' not found

抱歉，我为非洲使用了自定义的 SpatialPolygons（并且不知道如何附加文件）。尝试使用来自“tmap”的内置世界数据集，但无法重现该问题。无论如何，我目前也无法在原始数据/代码上重现问题，所以现在我希望这是我的数据/代码中的一个小故障，以某种方式修复（尽管我仍然不知道为什么问题发生在第一名）。

【参考方案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 将裁剪后的零件重新组合在一起。