我该如何从点数据计算一个区域的覆盖范围？

Question

我有几个月的数据文件，每个文件都包含两个25x25x20米的鱼场笼内40条已标记鱼在24小时内每6-9秒记录一次的鱼的x、y、z坐标。每个文件包含约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.1米，但即使分辨率设置为1米，运行时间也需要4小时；25x25米的网格数组= 625个单元，因此需要将包含365,000条鱼观测数据的坐标文件与网格数组交叉表。如果分辨率为0.1米，则365,000条观测数据需要与网格数组交叉表625,000次，这可能需要几周时间！

我相信肯定有更有效的方法来完成这个任务。然而，我学习R语言才几个月，不知道如何改进代码。

非常感谢您的帮助和建议！

- Adamaki

---

你可以通过先找到最北、最南、最东和最西的点，并直接将其外部的任何内容分配为0，来削减搜索空间的边缘。此外，您可以从更粗的分辨率开始，然后仅针对具有1的区域进行分辨率的细化（如果大正方形中没有鱼，则没有理由检查其较小的正方形）。但更好的方法可能是翻转您的过程：识别每条鱼的坐标，然后在您的网格中将它们绘制为1。 - ddunn801

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行和第列的框的占用次数，因此实际上比你的occupied提供更多信息，后者仅指定该框是否被占用（至少一次）。要检测一个框是否被占用，我们评估grid.cv > 0。

在0.1米分辨率网格上运行这个代码，使用365,000个模拟观测值的dayfile在我的2 GHz Macbook上只需要不到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

- aichao

---

太棒了，非常感谢！我已经尝试了你们两个的解决方案，它们都很好用，但是aichao的最快，因为它不使用循环。我不完全理解你的代码，但我看到它依赖于“floor”和“table”命令，这些我以前没有使用过，所以我会花一些时间来研究你的代码，以便下次知道该怎么做。现在我已经批处理了代码，可以运行60天的数据，只需要几分钟就能在我的3.2GHz电脑上完成。再次感谢你们的努力，我非常感激！ - Adamaki

Answer 2

这里是一种解决方案，您可以创建一个矩阵，用零表示网格，然后将每条鱼所在的单元格加1。然后区分具有1条或多条鱼和没有鱼的单元格，最后进行比例计算。我没有检查效率，但我想它会工作得更好（没有比较，只有一个 for）。

我将解决方案封装在一个函数中（更加优雅，可更轻松地应用于多个场合）。

请告诉我这是否适用于您！

dayfile<-data.frame(PosX=c(30.5,25.5,28.5), PosY=c(30,24,20))

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

coveragefun<-function(xmin, xmax, ymin, ymax, boxsize, dayfile){

  ncols <- ceiling((xmax-xmin)/boxsize)
  nrows <- ceiling((ymax-ymin)/boxsize)

  matspace <- matrix(0,nrow=nrows, ncol=ncols)

  for(i in 1:(dim(dayfile)[1])){
    xpos <- 1 + (dayfile$PosX[i]-(xmin))/boxsize
    ypos <- 1 + (dayfile$PosY[i]-(ymin))/boxsize
    matspace[xpos,ypos]<-matspace[xpos,ypos]+1
  }

  matcount<-matspace>=1

  coverage <- c(sum(matcount), dim(matcount)[1]*dim(matcount)[2], sum(matcount)/(dim(matcount)[1]*dim(matcount)[2]))
  names(coverage) <- c('occupied', 'total', 'proportion')
  return(coverage)
}

coverageres <- coveragefun(xmin, xmax, ymin, ymax, boxsize, dayfile)
coverageres

您还可以从函数中恢复matspace对象，以便进行摘要并了解网格中有多少填充单元。为此，您可以将代码的最后几行更改如下。

  return(list(coverage, matspace))
}

coverageres <- coveragefun(xmin, xmax, ymin, ymax, boxsize, dayfile)
coverageres[[1]]
table(coverageres[[2]])

- Rufo