R：如何计算网格的每个单元格中有多少点？答案

【问题标题】：R: How can I count how many points are in each cell of my grid?R：如何计算网格的每个单元格中有多少点？
【发布时间】：2015-12-29 14:56:56
【问题描述】：

我根据带项圈动物的 GPS 位置制作了一个参考网格，50x50 米的单元格。我想做相当于 ArcGIS 中的空间连接，并计算每个单元格中的点数。

我已经制作了一个参考网格，使用 SpatialPointsDataFrame 对象（数据框已经投影，使用 UTM 坐标系）

RESO <- 50 # grid resolution (m)
BUFF <- 500 # grid extent (m) (buffer around location extremes) 
XMIN <- RESO*(round(((min(dat.spdf$Longitude)-BUFF)/RESO),0))
YMIN <- RESO*(round(((min(dat.spdf$Latitude)-BUFF)/RESO),0))
XMAX <- XMIN+RESO*(round(((max(dat.spdf$Longitude)+BUFF-XMIN)/RESO),0))
YMAX <- YMIN+RESO*(round(((max(dat.spdf$Latitude)+BUFF-YMIN)/RESO),0))
NRW <- ((YMAX-YMIN)/RESO)
NCL <- ((XMAX-XMIN)/RESO)
refgrid<-raster(nrows=NRW, ncols=NCL, xmn=XMIN, xmx=XMAX, ymn=YMIN, ymx=YMAX) 
refgrid<-as(refgrid,"SpatialPixels")

为了确保我的网格与 SpatialPoints 处于同一投影中：

proj4string(refgrid)=proj4string(dat.sp.utm) #makes the grid the same CRS as point

adehabitatMA 中的count.point 函数似乎可以解决问题

cp<- count.points(dat.spdf, refgrid)

但我收到此错误：

Error in w[[1]] : no [[ method for object without attributes

这不是实现我的目标的正确途径吗？或者我该如何解决这个错误？还是over函数（sp包）更合适？

SpatialPointsDataFrame 的输出 (dat.spdf)

dput(head(dat.spdf, 20))
structure(list(Latitude = c(5.4118432, 5.4118815, 5.4115713, 
5.4111541, 5.4087853, 5.4083702, 5.4082527, 5.4078161, 5.4075528, 
5.407321, 5.4070598, 5.4064237, 5.4070621, 5.4070555, 5.4065127, 
5.4065134, 5.4064872, 5.4056724, 5.4038751, 5.4024223), Longitude = c(118.0225467, 
118.0222841, 118.0211875, 118.0208637, 118.0205413, 118.0206064, 
118.0204101, 118.0209272, 118.0213827, 118.0214189, 118.0217748, 
118.0223343, 118.0227079, 118.0226916, 118.0220733, 118.02218, 
118.0221843, 118.0223316, 118.0198153, 118.0196021), DayNo = c(1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 
2L, 2L, 2L)), .Names = c("Latitude", "Longitude", "DayNo"), row.names = c(1L, 
2L, 3L, 4L, 5L, 6L, 7L, 8L, 9L, 10L, 11L, 12L, 13L, 15L, 16L, 
17L, 18L, 19L, 20L, 21L), class = "data.frame")

总结：

summary(dat.spdf)
Object of class SpatialPointsDataFrame
Coordinates:
           min      max
Longitude 610361.0 613575.5
Latitude  596583.5 599385.2
Is projected: TRUE 
proj4string : [+proj=utm +zone=50 +ellps=WGS84]
Number of points: 5078
Data attributes:
Latitude       Longitude       DayNo      
Min.   :5.396   Min.   :118   Min.   :  1.0  
1st Qu.:5.404   1st Qu.:118   1st Qu.: 92.0  
Median :5.406   Median :118   Median :183.0  
Mean   :5.407   Mean   :118   Mean   :182.6  
3rd Qu.:5.408   3rd Qu.:118   3rd Qu.:273.0  
Max.   :5.422   Max.   :118   Max.   :364.0

【问题讨论】：

您是否确认refgrid 属于SpatialPixels 类（或继承类）？
是的，似乎是（对格式感到抱歉...）summary(refgrid) Object of class SpatialPixels Coordinates: min max x 609850 614100 y 596100 599900 Is projected: TRUE proj4string : [+proj=utm +zone=50 +ellps=WGS84] Number of points: 6460 Grid attributes: cellcentre.offset cellsize cells.dim s1 609875 50 85 s2 596125 50 76
如果您将refgrid 保留为栅格，则可以在table(cellFromXY(refgrid, dat.spdf)) 上构建。
谢谢@jbaums，我的表格看起来像每个单元格内的点数。我将如何进行空间显示？
我现在无法测试，但也许 r <- raster(refgrid); r[] <- 0; tab <- table(cellFromXY(refgrid, dat.spdf)); r[as.numeric(names(tab))] <- tab

标签： r spatial sp

【解决方案1】：

光栅化功能可以为您做到这一点：

library(raster)
r <- raster(xmn=0, ymn=0, xmx=10, ymx=10, res=1)
xy <- spsample(as(extent(r), 'SpatialPolygons'), 100, 'random')

x <- rasterize(xy, r, fun='count')

【讨论】：

【解决方案2】：

这是一种方法，首先将由点表示的单元格编号的频率制成表格，然后将这些频率分配给单元格的值，最后提取单元格的坐标和值。

library(raster)
r <- raster(xmn=0, ymn=0, xmx=10, ymx=10, res=1)
r[] <- 0
xy <- spsample(as(extent(r), 'SpatialPolygons'), 100, 'random')
tab <- table(cellFromXY(r, xy))
r[as.numeric(names(tab))] <- tab

现在我们有这样的东西：

plot(r)
points(xy, pch=20)

我们可以用coordinates()提取单元格的坐标，用values(r)或简单的r[]提取它们的值：

d <- data.frame(coordinates(r), count=r[])

head(d)
##     x   y count
## 1 0.5 9.5     0
## 2 1.5 9.5     1
## 3 2.5 9.5     1
## 4 3.5 9.5     3
## 5 4.5 9.5     2
## 6 5.5 9.5     3

【讨论】：