如何在 R 中正确裁剪（）栅格数据范围

Question

我正在尝试裁剪一些 raster data and do some calculations (，特别是）。

但是，在进行计算之前比较裁剪栅格数据的范围时，我得到的结果与在裁剪结果数据之前进行计算的结果相同。

栅格数据的原始范围是-180, 180, -90, 90 (xmin, xmax, ymin, ymax)，我需要将其裁剪到由纬度和经度坐标定义的任何所需区域。

这是我用来测试的脚本：

library(raster) # Crop raster data
library(stringr)

# hadsstR functions ----------------------------------------
load_hadsst <- function(file = "./HadISST_sst.nc") {
    b <- brick(file)
    NAvalue(b) <- -32768 # Land
    return(b)
}

# Transform basin coordinates into numbers
morph_coords <- function(coords){
    coords[1] = ifelse(str_extract(coords[1], "[A-Z]") == "W", - as.numeric(str_extract(coords[1], "[^A-Z]+")),
                                         as.numeric(str_extract(coords[1], "[^A-Z]+")) )
    coords[2] = ifelse(str_extract(coords[2], "[A-Z]") == "W", - as.numeric(str_extract(coords[2], "[^A-Z]+")),
                                         as.numeric(str_extract(coords[2], "[^A-Z]+")) )
    coords[3] = ifelse(str_extract(coords[3], "[A-Z]") == "S", - as.numeric(str_extract(coords[3], "[^A-Z]+")),
                                         as.numeric(str_extract(coords[3], "[^A-Z]+")) )
    coords[4] = ifelse(str_extract(coords[4], "[A-Z]") == "S", - as.numeric(str_extract(coords[2], "[^A-Z]+")),
                                         as.numeric(str_extract(coords[4], "[^A-Z]+")) )
    return(coords)
}

# Comparison test ------------------------------------------
hadsst.raster <- load_hadsst(file = "~/Hadley/HadISST_sst.nc")

x <- hadsst.raster
nms <- names(x)
months <- c("01","02","03","04","05","06","07","08","09","10","11","12")
coords <- c("85E", "90E", "5N", "10N")
coords <- morph_coords(coords)
years = 1970:1974
range = 5:12

# Crop before calculating mean
x <- crop(x, extent(as.numeric(coords[1]), as.numeric(coords[2]),
                                        as.numeric(coords[3]), as.numeric(coords[4])))
xMeans <- vector(length = length(years)-1,mode='list')
for (ix in seq_along(years[1:length(years)])){
    xMeans[[ix]] <- mean(x[[c(sapply(range,function(x) grep(paste0(years[ix],'.',months[x]),nms)))]], na.rm = T)
}
mean.brick1 <- do.call(brick,xMeans)

# Calculate mean before cropping
x <- hadsst.raster
xMeans <- vector(length = length(years)-1,mode='list')
for (ix in seq_along(years[1:length(years)])){
    xMeans[[ix]] <- mean(x[[c(sapply(range,function(x) grep(paste0(years[ix],'.',months[x]),nms)))]], na.rm = T)
}
mean.brick2 <- do.call(brick,xMeans)
mean.brick2 <- crop(mean.brick2, extent(as.numeric(coords[1]), as.numeric(coords[2]),
                                                                                as.numeric(coords[3]), as.numeric(coords[4])))

# Compare the two rasters
mean.brick1 - mean.brick2

这是mean.brick1 - mean.brick2的输出：

class       : RasterBrick 
dimensions  : 5, 5, 25, 5  (nrow, ncol, ncell, nlayers)
resolution  : 1, 1  (x, y)
extent      : 85, 90, 5, 10  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 
data source : in memory
names       : layer.1, layer.2, layer.3, layer.4, layer.5 
min values  :       0,       0,       0,       0,       0 
max values  :       0,       0,       0,       0,       0

如您所见，两个 RasterBricks 完全相同，对于任何任意选择的坐标来说这应该是不可能的，如下面的小矩阵示例：

我做错了什么吗？ 在对数据进行计算之前裁剪数据应该会明确地给出不同的结果。

Answer 1

好的，我将从您上一个问题中的继续：

我们从完整的 hadsst.raster 积木开始（为了获得可重现的示例，可以使用我之前回答中的解决方案的第一部分进行伪造）。

因此此数据集的维度为 180, 360, 516，即 180 行、360 列和 516 个时间层。

从技术上讲，栅格是一个矩阵，它可能是这样的：

只是一堆矩阵层（准确地说是 516 个），其中每个像素都精确对齐。这里我只有三个示例图层，其余的用三个点表示。

因此，如果我们进行时间平均，我们基本上会提取单个像素的所有值并取它们的平均值（或任何其他平均操作）。此处用红色方块表示。

这也说明了为什么裁剪不会影响时间平均：

如果我们说橙色方块是我们感兴趣的范围并且我们在平均之前执行裁剪操作，我们基本上会丢弃这个方块周围的所有值。之后，我们再次获取所有层上每个像素的所有值并计算平均值。

现在应该清楚了，为什么丢弃橙色方块周围的像素并不重要。您还可以计算它们的平均值，然后丢弃这些值，只剩下橙色方块的值。如果您已经确定不需要它们进行进一步计算，那么它就没有任何实际意义。无论如何，方块内的值不会受到影响。

当我们谈论空间平均时，它通常意味着对单个图层内的像素进行平均，在这种情况下可能是对橙色矩形内的值进行平均。

两个常见的操作是

焦点平均（也称为邻域平均）
聚合

焦点平均将对每个像素取定义数量的相邻像素的所有值的平均值（最常见的是 3x3 正方形，其中要定义的像素是中心像素） .

聚合实际上是将多个像素组合成一个更大的像素。这意味着不仅该像素的值将被平均，而且生成的栅格将具有更少的单个像素和更粗糙的分辨率。

好的，为您提供实际解决方案：

我假设您有一个由范围定义的感兴趣区域 aoi:

aoi <- extent(xmin,xmax,ymin,ymax)

您要做的第一件事是裁剪初始砖块以减少计算负担：

hadsst.raster_crp <- crop(hadsst.raster,aoi)

下一步是时间平均，我们使用我在其他 :

的解决方案中定义的函数

hadsst.raster_crp_avg <- hadSSTmean(hadsst.raster_crp, 1969:2011, first.range = 11:12, second.range = 1:4)

好的，现在您有了您感兴趣区域的时间平均值。下一步取决于您的最终目标。据我了解，对于您感兴趣的区域，您只需要每个时间平均值的单个平均值。

如果是这样，可能是时候离开实际的栅格域并继续使用基本 R:

res <- lapply(1:nlayers(hadsst.raster_crp_avg),function(ix) mean(as.matrix(hadsst.raster_crp_avg[[ix]])))

这将为您提供一个包含与您的积木 hadsst.raster_crp_avg 一样多的元素的列表。

使用 lapply，我们遍历层，将每一层转换为矩阵，然后计算所有元素的平均值，为整个感兴趣区域留下每个平均时间步长的单个值。

更进一步，您可以使用 unlist 将其转换为矢量并将其添加到 data.frame 或执行您喜欢的任何其他操作。

希望这很清楚，这就是您要找的东西。

最佳

如何在 R 中正确裁剪（）栅格数据范围

How to properly crop() raster data extent in R

r

r-raster