R：如何使用特定缓冲区循环遍历空间点？

Question

所以我的问题很难描述所以我希望我能把我的问题尽可能清楚。
我使用 rLiDAR package 将 .las 文件加载到 R 中，然后使用 sp package 将其转换为 SpatialPointsDataFrame。所以我的 SpatialPointsDataFrame 非常密集。

现在我想定义一个0.5米的缓冲区，和他（缓冲区）一起循环（迭代）遍历这些点，总是选择缓冲区内Z值最高的点，作为下一个要跳转的点to.This 应该重复，直到缓冲区内没有任何点具有更高的 Z 值作为当前点。 这个"found"点的所有值（或者可能是X和Y值）应该被写入list/dataframe并且应该重复这个过程直到所有这样的最高值点被发现。那是我到目前为止得到的代码：

>library(rLiDAR)
>library(sp)
>rLAS<-readLAS("Test.las",short=FALSE)
>PointCloud<- data.frame(rLAS) 
>coordinates(PointCloud) <- c("X", "Y")

好吧，我在谷歌上进行了广泛的搜索，但我找不到任何关于如何进一步进行的线索...
我什至不知道哪些包可以提供帮助，我想也许 spatstat 因为我的问题可能会进入空间点模式分析。
有没有人知道如何在 R 中归档类似的东西？或者这样的事情是不可能的？（我是否必须跳到 python 才能完成这样的工作？）
很乐意提供帮助。

Answer 1

如果您想获得一组点，这些点是每个点周围 0.5 米半径圆内的局部最大值，这应该可行。它的要点是：

将 LAS 点转换为 SpatialPointsDataFrame
创建具有重叠多边形的缓冲多边形集
遍历所有缓冲的多边形并在缓冲区中找到所需的元素——在您的例子中，它是具有最大高度的元素。

代码如下：

library(rLiDAR)
library(sp)
library(rgeos)

rLAS <- readLAS("Test.las",short=FALSE)
PointCloud <- data.frame(rLAS) 
coordinates(PointCloud) <- c("X", "Y")

完成从 LAS 源创建 SpatialPointsDataFrame。我假设具有点高的字段是 PointCloud$value

pointCloudSpdf <- SpatialPointsDataFrame(data=PointCloud,xy)

使用 rgeos 库进行交集。 byid=TRUE 很重要，否则多边形将在它们相交的地方合并

bufferedPoints <- gBuffer(pointCloudSpdf,width=0.5,byid=TRUE)

# Save our local maxima state (this will be updated)
localMaxes <- rep(FALSE,nrow(PointCloud))
i=0
for (buff in 1:nrow(bufferedPoint@data)){
  i <- i+1
  bufPolygons <- bufferedPoints@polygons[[i]]
  bufSpPolygons <- SpatialPolygons(list(bufPolygons))
  bufSpPolygonDf <-patialPolygonsDataFrame(bufSpPolygons,bufferedPoints@data[i,])
  ptsInBuffer <- which(!is.na(over(pointCloudSpdf,spPolygonDf)))

  # I'm assuming `value` is the field name containing the point height
  localMax <- order(pointCloudSpdf@data$value[ptsInBuffer],decreasing=TRUE)[1]
  localMaxes[localMax] <- TRUE
}

localMaxPointCloudDf <- pointCloudSpdf@data[localMaxes,]

现在 localMaxPointCloudDf 应该包含来自原始点的数据（如果它们是局部最大值）。只是一个警告——如果你有很多点，这不会很快。如果这最终成为一个问题，您可能会更聪明地使用较小的网格和 raster 包中的 extract 预过滤您的点。

看起来像这样：

使像元大小足够小，以便每个 0.5m 的缓冲区将与至少 4 个栅格像元相交 -- 较小的错误，因为我们将圆形与正方形进行比较。

library(raster)
numRows <- extent(pointCloudSpdf)@ymax-extent(pointCloudSpdf)@ymin/0.2
numCols <- extent(pointCloudSpdf)@xmax-extent(pointCloudSpdf)@xmin/0.2
emptyRaster <- raster(nrow=numRows,ncol=numCols)

rasterize 将在一个单元格中创建一个具有给定字段最大值的网格。由于 square/circle 不匹配，这只是过滤掉明显非最大值的起点。在此之后，我们将得到一个栅格，其中所有局部最大值都由像元表示。但是，我们不知道在 0.5m 半径内哪些像元是最大值，也不知道它们来自原始特征层中的哪个点。

r <- rasterize(pointCloudSpdf,emptyRaster,"value",fun="max")

extract 将为我们提供每个点相交的栅格值（即每个像元的最大值）。从上面回想一下，所有局部最大值都将在此集合中，尽管有些值不会是 0.5m 半径的局部最大值。

rasterMaxes <- extract(r,pointCloudSpdf)

要使原始点与栅格最大值相匹配，只需从每个点的值中减去栅格值即可。如果该值为 0，则值相同，并且我们有一个具有潜在最大值的点。请注意，此时我们只是将这些点合并回栅格——我们将不得不丢弃其中的一些，因为它们是 "under" 一个 0.5m 的半径，具有更高的局部最大值，即使它们是最大的他们的 0.2m x 0.2m 单元格。

potentialMaxima <- which(pointCloudSpdf@data$value-rasterMaxes==0)

接下来，只需对原始 SpatialPointsDataFrame 进行子集化，我们将对这个点子集进行更详尽和更准确的迭代，因为我们应该丢弃一堆不可能是最大值的点。

potentialMaximaCoords <- coordinates(pointCloudSpdf@coords[potentialMaxima,])
# using the data.frame() constructor because my example has only one column
potentialMaximaDf <- data.frame(pointCloudSpdf@data[potentialMaxima,])
potentialMaximaSpdf <-SpatialPointsDataFrame(potentialMaximaCoords,potentialMaximaDf)

算法的其余部分是相同的，但我们正在缓冲较小的数据集并对其进行迭代：

bufferedPoints <- gBuffer(potentialMaximaSpdf, width=0.5, byid=TRUE)

# Save our local maxima state (this will be updated)
localMaxes <- rep(FALSE, nrow(PointCloud))
i=0
for (buff in 1:nrow(bufferedPoint@data)){
  i <- i+1
  bufPolygons <- bufferedPoints@polygons[[i]]
  bufSpPolygons <- SpatialPolygons(list(bufPolygons))
  bufSpPolygonDf <-patialPolygonsDataFrame(bufSpPolygons,bufferedPoints@data[i,])
  ptsInBuffer <- which(!is.na(over(pointCloudSpdf, spPolygonDf)))

  localMax <- order(pointCloudSpdf@data$value[ptsInBuffer], decreasing=TRUE)[1]
  localMaxes[localMax] <- TRUE
}

localMaxPointCloudDf <- pointCloudSpdf@data[localMaxes,]

R：如何使用特定缓冲区循环遍历空间点？

R: How do I loop through spatial points with a specific buffer?

r

loops

buffer

spatial

point