Hi,
I would like to count the number of points per pixel in a grid(there are
thousands of points and thousands of pixels). I have a list of occurrence
points and a grid with a given cell size. Ideally I'd like them saved to
one file where there is a column listing the pixel number and one listing
number of points for that pixel. I've tried overlay (sp package) but it
just spits out a list of all the points in my list with NA value attached
to each point and it doesn't provide a count:
this is my code so far:
grid=readGDAL("mask50")
proj4string(grid)=CRS("+proj=lcc")
x1<-read.csv(file.choose(),header=TRUE,sep=",", na.strings="");
proj4string(x1)=CRS("+proj=lcc")
o<-overlay(grid, x1)
o
NA.33354 (49.71, -97.91) NA
NA.33355 (50.7, -98.03) NA
NA.33356 (56.35, -94.66) NA
NA.33357 (50.83, -100) NA
NA.33358 (49.36, -96.11) NA
I'm wondering whether I need to set up a loop for each pixel.
I've also had no luck with count.points.id (adehabitat package). When I
use this function, I get a summary of my grid. Any thoughts would be much
appreciated!
Thanks,
Heather Kharouba
Counting number of points per pixel in a grid
3 messages · Heather Kharouba, Roger Bivand, Paul Hiemstra
On Tue, 20 Oct 2009, Heather Kharouba wrote:
Hi,
I would like to count the number of points per pixel in a grid(there are
thousands of points and thousands of pixels). I have a list of occurrence
points and a grid with a given cell size. Ideally I'd like them saved to
one file where there is a column listing the pixel number and one listing
number of points for that pixel. I've tried overlay (sp package) but it
just spits out a list of all the points in my list with NA value attached
to each point and it doesn't provide a count:
this is my code so far:
grid=readGDAL("mask50")
proj4string(grid)=CRS("+proj=lcc")
x1<-read.csv(file.choose(),header=TRUE,sep=",", na.strings="");
# here you seem to be missing a coordinates(x1) <- whatever
proj4string(x1)=CRS("+proj=lcc")
o<-overlay(grid, x1)
o
This is a SpatialPointsDataFrame, with the input coordinates of x1, the row names of grid, and the values of grid (all visible are NA, which maybe is reasonable for a mask).
NA.33354 (49.71, -97.91) NA NA.33355 (50.7, -98.03) NA NA.33356 (56.35, -94.66) NA NA.33357 (50.83, -100) NA NA.33358 (49.36, -96.11) NA I'm wondering whether I need to set up a loop for each pixel. I've also had no luck with count.points.id (adehabitat package). When I use this function, I get a summary of my grid. Any thoughts would be much appreciated!
No, saying summary(o) would show that you get a summary of the grid values at the points in x1. If you just want the grid cell indices for each point, coerce first: o <- overlay(as(grid, "SpatialGrid"), x1)
From that you can work out how many points there are in each grid cell,
most likely table(o) gets most of the way there, with the table being the counts and as.integer(names(table(o))) being the indices. Then: tab <- table(o) grid$counts <- as.integer(NA) grid$counts[as.integer(names(tab))] <- tab summary(grid) Hope this helps, Roger PS. Your point coordinate values look geographic (in Canada?), not projected, and probably (lat, long), not the required (long, lat); "+proj=lcc" values would be in metres.
Thanks, Heather Kharouba
_______________________________________________ R-sig-Geo mailing list R-sig-Geo at stat.math.ethz.ch https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Roger Bivand Economic Geography Section, Department of Economics, Norwegian School of Economics and Business Administration, Helleveien 30, N-5045 Bergen, Norway. voice: +47 55 95 93 55; fax +47 55 95 95 43 e-mail: Roger.Bivand at nhh.no
Hi Heather, Be sure to take a look a the raster package, in particular the pointsToRaster function. Does exactly what you want and is very fast. Going from sp-grids to raster-objects is very easy. cheers, Paul Heather Kharouba schreef:
Hi,
I would like to count the number of points per pixel in a grid(there are
thousands of points and thousands of pixels). I have a list of occurrence
points and a grid with a given cell size. Ideally I'd like them saved to
one file where there is a column listing the pixel number and one listing
number of points for that pixel. I've tried overlay (sp package) but it
just spits out a list of all the points in my list with NA value attached
to each point and it doesn't provide a count:
this is my code so far:
grid=readGDAL("mask50")
proj4string(grid)=CRS("+proj=lcc")
x1<-read.csv(file.choose(),header=TRUE,sep=",", na.strings="");
proj4string(x1)=CRS("+proj=lcc")
o<-overlay(grid, x1)
o
NA.33354 (49.71, -97.91) NA
NA.33355 (50.7, -98.03) NA
NA.33356 (56.35, -94.66) NA
NA.33357 (50.83, -100) NA
NA.33358 (49.36, -96.11) NA
I'm wondering whether I need to set up a loop for each pixel.
I've also had no luck with count.points.id (adehabitat package). When I
use this function, I get a summary of my grid. Any thoughts would be much
appreciated!
Thanks,
Heather Kharouba
_______________________________________________ R-sig-Geo mailing list R-sig-Geo at stat.math.ethz.ch https://stat.ethz.ch/mailman/listinfo/r-sig-geo