I have made a reference grid, cells 50x50m, based on GPS locations of a collared animal. I want to do the equivalent to a spatial join in ArcGIS, and count the number of points in each cell.
I have made a reference grid, using a SpatialPointsDataFrame object (the data frame is already projected, using UTM coordinate system)
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")
To make sure my grid was in the same projection as the SpatialPoints:
proj4string(refgrid)=proj4string(dat.sp.utm) #makes the grid the same CRS as point
count.point
function in adehabitatMA seems like the function that will do the trick
cp<- count.points(dat.spdf, refgrid)
But I get this error:
Error in w[[1]] : no [[ method for object without attributes
Is this not the right route to take to achieve my goal? Or how can I resolve this error? Or would the over
function (sp package) be more suitable?
output from 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")
And summary:
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
The rasterize function can do that for you:
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')
Here's one way to do it, first tabulating the frequency of cell numbers represented by points, then assigning these frequencies to the cells' values, and finally extracting the cells' coordinates and values.
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
Now we have something like this:
plot(r)
points(xy, pch=20)
We can extract the cells' coordinates with coordinates()
and their values with values(r)
or simply 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
If you love us? You can donate to us via Paypal or buy me a coffee so we can maintain and grow! Thank you!
Donate Us With