r - 如何根据点数据计算区域的覆盖范围？

Question

我有几个月的数据文件，每个文件包含两个 25x25x20m 养鱼围栏中 40 条标记鱼的鱼 x、y、z 坐标的 24 小时记录，每个标记每 6-9 秒定位一次。每个文件包含大约 365,000 个观察值。

我想计算鱼每天覆盖围栏的比例。我已经编写了一些 R 代码来完成这项工作，但由于文件很大，运行大约需要 4 个小时。这是我的代码:

xmin <- 8
ymin <- 11.5
xmax <- 33
ymax <- 36.5
boxsize <- 1

# define coverage grid
cov.grid <- matrix(c(xmin,ymin), nrow = 1, ncol = 2, byrow = FALSE)
colnames(cov.grid) <- c('x','y')
x <- xmin
y <- ymin
while(x < xmax)
  {
  while(y < ymax)
    {
    y <- y+boxsize
    cov.grid <- rbind(cov.grid, c(x,y))  
    }
  x <- x+boxsize
  y <- ymin
  cov.grid <- rbind(cov.grid, c(x,y))  
}
cov.grid <- as.data.frame(cov.grid)


# count grid cells occupied by fish
day.row <- 1
grid.row <- 1
bin <- 0
cov.grid$occupied <- NA

for(grid.row in 1:nrow(cov.grid)){
x1 <- cov.grid[grid.row,1]
y1 <- cov.grid[grid.row,2]
x2 <- x1+boxsize
y2 <- cov.grid[grid.row+1,2] 
repeat
  {
  if(dayfile[day.row,'PosX'] > x1 & dayfile[day.row,'PosX'] < x2 &         dayfile[day.row,'PosY'] > y1 & dayfile[day.row,'PosY'] < y2) {bin <- 1} else    {bin <- 0}
  day.row <- day.row+1
  if(bin == 1 | day.row == nrow(dayfile)){break}
  }
cov.grid[grid.row,'occupied'] <- bin
day.row <- 1
}

# return coverage summary

coverage <- matrix(c(length(which(cov.grid$occupied == 1)), nrow(cov.grid),     length(which(cov.grid$occupied == 1))/nrow(cov.grid)), ncol = 3)
colnames(coverage) <- c('occupied', 'total', 'proportion')
coverage

代码逻辑如下:

1. 创建笔区的矩阵网格。
2. 对于每个网格单元格，查看鱼坐标文件以检查是否有鱼占据该单元格；如果是，则为 1，如果否，则为 0。
3. 在网格矩阵中添加一个新列来记录每个单元格是否被一条鱼占据。
4. 统计占用的单元格数，计算笔的覆盖比例。

理想情况下，我希望网格分辨率为 0.1m 分辨率，但即使是 1m 分辨率也需要 4 小时才能运行； 25x25m 网格阵列 = 625 个单元格，因此 365,000 条鱼类观测值的坐标文件必须与网格阵列交叉制表 625 次。如果网格分辨率为 0.1m，则 365,000 个观测值需要交叉制表 625,000 次，这可能需要数周时间!

我相信一定有更有效的方法来做到这一点。但是，我现在才学习 R 几个月，所以我不确定如何改进代码。

如有任何帮助或建议，我们将不胜感激!

Answer 1

您根本不需要使用循环。执行以下操作:

compute.coverage <- function(xmin, xmax, ymin, ymax, boxsize, dayfile) {
  x.grid <- floor((dayfile$PosX - xmin) / boxsize) + 1
  y.grid <- floor((dayfile$PosY - ymin) / boxsize) + 1
  x.grid.max <- floor((xmax - xmin) / boxsize) + 1
  y.grid.max <- floor((ymax - ymin) / boxsize) + 1
  t.x <- sort(unique(x.grid))
  t.y <- sort(unique(y.grid))
  tx.range <- c(min(which(t.x > 0)), max(which(t.x <= x.grid.max)))
  ty.range <- c(min(which(t.y > 0)), max(which(t.y <= y.grid.max)))
  t <- table(y.grid, x.grid)[ty.range[1]:ty.range[2],tx.range[1]:tx.range[2]]
  grid.cov <- matrix(0,nrow=y.grid.max,ncol=x.grid.max)
  t.x <- t.x[(t.x > 0) & (t.x <=x.grid.max)]
  t.y <- t.y[(t.y > 0) & (t.y <=y.grid.max)]
  eg <- expand.grid(t.y,t.x)
  grid.cov[cbind(eg$Var1,eg$Var2)] <- as.vector(t)  
  coverage <- matrix(c(length(which(grid.cov > 0)), length(grid.cov), length(which(grid.cov > 0))/length(grid.cov)), ncol = 3)
  colnames(coverage) <- c('occupied', 'total', 'proportion')
  coverage
}

此计算的关键是像 Rufo(另一个答案)那样为每个观察计算网格框位置 (x.grid,y.grid)。然而，这里的计算是向量化所有 dayfile 中的观察，其复杂度独立的分辨率网格!诀窍是然后使用 table 计算每个 (x.grid,y.grid) 组合的占用频率。这里有两个复杂的因素:

1. 计算出的 (xgrid,y.grid) 位置可能在您的笔外 (xmin,xmax,ymin,ymax)。
2. 并非所有的网格框都被占用，因此表格中可能有整行和/或整列的计数缺失。

如果你只对覆盖百分比感兴趣，第二个问题是不相关的，但如果你真正关心哪个框位置被占用，它是相关的。上面的代码通过以下方式处理:

1. 将表格限制在笔的范围内，tx.range 和 ty.range。
2. 将表格(可能带有“孔”)映射回笔的完整网格 grid.cov。此处，grid.cov 是对应于您的 cov.grid 变量的笔矩阵。它的元素记录了第 i 行和 j 列的盒子的占用数量，所以这实际上比你的 occupied，它只指定该框是否已被占用(至少一次)。为了检测一个框是否已被占用，我们评估 grid.cv > 0。

在我的 2 GHz Macbook 上用 365,000 次模拟观察在 dayfile 上以 0.1 米分辨率的网格运行此程序花费了不到 2 秒:

xmin <- 8
ymin <- 11.5
xmax <- 33
ymax <- 36.5
boxsize <- 0.1

## simulate dayfile
set.seed(123)
PosX <- runif(365000,xmin-2,xmax+2)
PosY <- runif(365000,ymin-2,ymax+2)
dayfile <- data.frame(PosX=PosX,PosY=PosY)

print(system.time(coverage <- compute.coverage(xmin,xmax,ymin,ymax,boxsize,dayfile)))
##   user  system elapsed 
##  1.096   0.052   1.193 

print(coverage)
##     occupied total proportion
##[1,]    62168 63001   0.986778