从栅格中沿自定义圆形条带（恒定径向距离）提取值

Question

我有一个从 .tif 图像文件导入的大型光栅对象。我想在距中心给定径向距离处分析该栅格的像素，以识别我可以在图像中注意到的某些轴对称现象（下图）。为此，我想从图像中心（如下图所示）提取与给定半径（和宽度）的圆形条带相交的像素值。

我已经探索了几个选项来执行此操作，包括提取功能和成像器包。在 imager 包中，您可以轻松地沿行或列提取值，但我找不到像我需要的那样提取具有自定义形状的值的函数。

我可能会使用带有 SpatialPolygons 参数的 extract 函数，但是为此我必须输入多边形对象中所有点的位置，这将需要非常高密度的点位置（因为我认为这些点，由线元素连接）。此外，我想改变我提取像素值（以及后来的平均值）的条带的数量密度，因此这种方法既乏味又不灵活。因此，我想知道你们是否有任何建议来解决这个问题。

> str(imRaster)
Formal class 'RasterLayer' [package "raster"] with 12 slots
  ..@ file    :Formal class '.RasterFile' [package "raster"] with 13 slots
  .. .. ..@ name        : chr "C:\\Users\\Nandu\\input_images\\All\\101_1A_1000ms.tif"
  .. .. ..@ datanotation: chr "INT2U"
  .. .. ..@ byteorder   : chr "little"
  .. .. ..@ nodatavalue : num -Inf
  .. .. ..@ NAchanged   : logi FALSE
  .. .. ..@ nbands      : int 1
  .. .. ..@ blockrows   : int 1
  .. .. ..@ blockcols   : int 1392
  .. .. ..@ driver      : chr "gdal"
  .. .. ..@ open        : logi FALSE
  ..@ data    :Formal class '.SingleLayerData' [package "raster"] with 13 slots
       .
       .
       .

或者，如果有其他方法（除了将 .tif 图像导入光栅之外）可以允许这样的操作，那也可能有用。请告诉我。

提前致谢！

编辑：为了使问题可重现，我在此处进一步添加了指向我要使用的图像的链接。这是一个tiff图像，可以从这里下载

Link to the tiff image that I am trying to process in R

library(tiff)
library(raster)
imRaster = raster(file_path to the image)
plot(imRaster)
xy = cbind(684.4228, 599.0458) # Gives a rough location of the center in the image coordinates

可以运行此代码以获得上述栅格并将其可视化。希望对您有所帮助！

Answer 1

让我们首先加载栅格并为其分配一个假坐标参考系统（CRS）。我们使用度量单位，因此您图像中的 1 个像素指的是 pseudo-地理栅格中的一个“米”。

library(raster)

imRaster <- raster("~/Downloads/46-ECN300448-216.tif")
imRaster@crs <-  CRS("+proj=merc +lon_0=0 +k=1 +x_0=0 +y_0=0 +a=6378137 +b=6378137 +units=m +no_defs")

plot(imRaster, col = grey.colors(255))

之后，正如您提到的，我们使用您的中心点 xy 并将其转换为 SpatialPoint

xy <- cbind(684.4228, 599.0458)
xySp <- SpatialPoints(coords = xy)
plot(xySp, add = T)

然后我们使用rgeos 包在该点周围创建两个径向缓冲区。在此示例中，我们使用 200 和 300 px/m。

library(rgeos)

outer <- 300
inner <- 200

bufOut <- gBuffer(xySp, width = outer)
bufIn <- gBuffer(xySp, width = inner)

strip <- bufOut - bufIn
plot(strip, add = T, col = "#FF000050")

最后，我们可以使用strip 屏蔽光栅图像并使用getValues() 计算统计数据（甚至使用原始值）。

m <- mask(imRaster, mask = strip)

# plot(m) # plot the mask if you want to see what it looks like

mean(getValues(m), na.rm = T)
# [1] 17004.24

重要提示：蒙版后仅保留中心位于多边形内的单元格。您可以通过使用缓冲区来解决这个问题。