使用 r 根据最近的像元值从特定点的栅格中提取值

Question

我正在尝试为某些站点数据分配值，这些数据恰好位于我拥有天气数据的区域之外。我正在尝试根据最近的单元格值提取，如果可能的话，提取 40 公里内的单元格值。

我的光栅 (r) 如下所示：

class(r)

class       : RasterBrick 
dimensions  : 201, 464, 93264, 23376  (nrow, ncol, ncell, nlayers)
resolution  : 0.25, 0.25  (x, y)
extent      : -40.5, 75.5, 25.25, 75.5  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 
data source : \\ueahome\eressci5\zuw13bqu\data\NTProfile\Desktop\EOBS European data\rr_0.25deg_reg_v10.0.nc 
names       : X1950.01.01, X1950.01.02, X1950.01.03, X1950.01.04, X1950.01.05, X1950.01.06, X1950.01.07, X1950.01.08, X1950.01.09, X1950.01.10, X1950.01.11, X1950.01.12, X1950.01.13, X1950.01.14, X1950.01.15, ... 
Date        : 1950-01-01, 2013-12-31 (min, max)
varname     : rr 

我使用以下代码根据纬度和经度数据提取

vals <- extract(r, matrix(c(issues[22,3], issues[22,2]), ncol = 2), buffer = 40000)

但不幸的是我得到以下输出：

*无法附上图片，因为我没有足够的声誉

X1950.01.01 X1950.01.02 X1950.01.03 X1950.01.04 X1950.01.05 X1950.01.06 X1950.01.07 X1950.01.08 X1950.01.09 X1950.01.10
1   0   4.8 4.6 0   0.2 0   0   0   0   0
2   0   4.7 4.5 0   1   0   0   0   0   0
3   0   4.7 4.5 0   1.1 0   0   0   0   0
4   0   4.6 4.3 0   1.2 0   0   0   0   0
5   NA  NA  NA  NA  NA  NA  NA  NA  NA  NA
6   NA  NA  NA  NA  NA  NA  NA  NA  NA  NA
7   NA  NA  NA  NA  NA  NA  NA  NA  NA  NA
8   0   4.1 3.9 0   0.7 0   0   0   0   0
9   0   4   3.7 0   0.9 0   0   0   0   0
10  0   4.1 3.8 0   1   0   0   0   0   0
11  0   4.1 3.8 0   1.1 0   0   0   0   0
12  NA  NA  NA  NA  NA  NA  NA  NA  NA  NA
13  NA  NA  NA  NA  NA  NA  NA  NA  NA  NA
14  NA  NA  NA  NA  NA  NA  NA  NA  NA  NA
15  NA  NA  NA  NA  NA  NA  NA  NA  NA  NA
16  NA  NA  NA  NA  NA  NA  NA  NA  NA  NA

*nb 我已经检查过这个网站，这些行实际上都不是最接近的。

如何选择最接近该点的单元格值而不减小缓冲区大小，直到出现一个单元格值（我有太多这样的站点，无法为每个站点执行此操作）？

提前致谢

Answer 1

r当找到到外部点的最小距离的栅格像元时r，我们不需要计算到 内部像元的距离r。

具体来说：

1. 对于具有x < xmin(r)和 的点y < ymin(r)，最接近的r是左下单元格；
2. 对于具有x > xmax(r)和 的点y < ymin(r)，最接近的r是右下单元格；
3. 对于具有x < xmin(r)和 的点y > ymax(r)，最接近的r是左上角的单元格；
4. 对于具有x > xmax(r)和 的点y > ymax(r)，最接近的r是右上单元格；
5. 对于具有xmin(r) < x < xmax(r)和 的点y < ymin(r)，最近的单元格r将位于底部边缘；
6. 对于具有xmin(r) < x < xmax(r)和 的点y > ymax(r)，最近的单元格r将位于顶部边缘；
7. 对于具有x < xmin(r)和 的点ymin(r) < y < ymax(r)，最近的单元格r将位于左边缘；
8. 对于具有xmin(r) < x < xmax(r)和 的点y > ymax(r)，最近的单元格r将位于顶部边缘；和
9. xmin(r) < x < xmax(r)和ymin(r) < y < ymax(r)已经在栅格内的点。

下面的函数可以计算出每个点属于这 9 个区域中的哪一个，识别它们对应的栅格列或行，并识别沿相关边缘/角的该列或行的像元。

该函数接受2 x n坐标矩阵（xthen y）或SpatialPoints对象。设置spatial=TRUE是否希望移动的点返回为SpatialPoints。

请注意，假设点位于平面坐标系中。

move.points <- function(r, pts, spatial=FALSE) {
  require(raster)
  require(sp)

  if (is(pts, 'SpatialPoints')) pts <- coordinates(pts)
  if (is(!r, 'Raster')) r <- raster(r)

  loc <- colSums(sapply(pts[, 1], '>', bbox(r)[1, ])) * 3 + 
    colSums(sapply(pts[, 2], '>', bbox(r)[2, ]))

  L <- split(as.data.frame(pts), loc)

  new.pts <- lapply(names(L), function(x) {
    switch(x, 
           '0' = xyFromCell(r, ncell(r) - ncol(r) + 1)[rep(1, nrow(L[[x]])), ],
           '1' = xyFromCell(r, cellFromXY(r, cbind(xmin(r), L[[x]][, 2]))),
           '2' = xyFromCell(r, 1)[rep(1, nrow(L[[x]])), ],
           '3' = xyFromCell(r, cellFromXY(r, cbind(L[[x]][, 1], ymin(r)))),
           '4' = {
             xy <- as.matrix(L[[x]])
             dimnames(xy) <- list(NULL, c('x', 'y'))
             xy
           },
           '5' = xyFromCell(r, cellFromXY(r, cbind(L[[x]][, 1], ymax(r)))),
           '6' = xyFromCell(r, ncell(r))[rep(1, nrow(L[[x]])), ],
           '7' = xyFromCell(r, cellFromXY(r, cbind(xmax(r), L[[x]][, 2]))),
           '8' = xyFromCell(r, ncol(r))[rep(1, nrow(L[[x]])), ]
    )
  })

  new.pts <- unsplit(mapply(function(x, y) {
    row.names(x) <- row.names(y)
    as.data.frame(x)
  }, new.pts, L, SIMPLIFY=FALSE), loc)

  colnames(new.pts) <- colnames(pts)
  if(isTRUE(spatial)) new.pts <- SpatialPoints(new.pts)  
  return(new.pts)
}

以下是随机放置 1000 个点的 1000 x 1000 像元栅格的示例（约 890 个像元落在栅格范围之外）：

r <- raster(matrix(runif(1e6), ncol=1e3))
pts <- cbind(lon=runif(1e3, -1, 2), lat=runif(1e3, -1, 2))
system.time(new.pts <- move.points(r, pts, spatial=FALSE))

#   user  system elapsed 
#   0.03    0.00    0.04 

对于 2000 x 2000 像元栅格，有 100 万个点（大约 890000 个点位于栅格范围之外）：

r <- raster(matrix(runif(4e6), ncol=2e3))
pts <- cbind(lon=runif(1e6, -1, 2), lat=runif(1e6, -1, 2))
system.time(new.pts <- move.points(r, pts, spatial=FALSE))

#   user  system elapsed 
#  12.71    0.23   13.12 

下面的图显示了上面第一个示例中的点是如何移动的。

plot(pts, asp=1, pch=20, panel.first=abline(h=0:1, v=0:1, lwd=2, col='gray'))
segments(pts[, 1], pts[, 2], new.pts[, 1], new.pts[, 2], col='#00000050')
plot(extent(r), add=TRUE, lwd=3)

完成此操作后，您就可以使用extract(r, new.pts).