如何從點數據計算區域的覆蓋範圍？

Question

我有幾個月的數據文件，每個數據文件包含24個小時的魚x，y，z座標，記錄在兩個25x25x20m魚場圍欄中，用於40個標記的魚，每個標記每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.1米分辨率，但即使是1米分辨率也需要4個小時才能運行; 25×25米的網格陣列= 625個單元，因此365,000條魚觀測的座標文件必須與網格陣列625次交叉製表。如果網格分辨率爲0.1米，那麼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 
} 

這個計算的關鍵是計算每個觀測作爲入佛門的格箱位置(x.grid,y.grid)（對方的回答）沒有。然而，在這裏，這個計算是矢量化超過所有意見在dayfile和它的複雜性是獨立的分辨率的網格！訣竅是然後使用table來計算在每個組合(x.grid,y.grid)佔領的頻繁程度。這裏有兩個複雜的因素：

1. 的計算(xgrid,y.grid)位置可能是你的筆(xmin,xmax,ymin,ymax)之外。
2. 並非您所有的網格框都被佔用，因此可能會有整個表格的行和/或列缺少計數。

第二個問題是不恰當的，如果你只是在覆蓋率的百分比感興趣，但它是相關的，如果你真的關心哪個箱位置被佔用。上面的代碼通過同時處理：

1. 限制表格的範圍內，和tx.rangety.range，筆的。
2. 將表（可能帶有「孔」）映射回到筆的全網格grid.cov。這裏，grid.cov是與您的cov.grid變量對應的筆的矩陣。它的元素記錄了第i列第i列的箱子的職業數量，所以這實際上比occupied更多的信息，它只指定箱子是否已被佔用（至少一次）。爲了檢測箱子是否被佔用，我們評估grid.cv > 0。

在dayfile與365000模擬觀測上0.1米分辨率的網格上運行這個花了不到2秒，我的2千兆赫的MacBook：

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 

Answer 2

下面是一個解決方案，您可以創建一個表示網格的零點的矩陣，然後將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]])