[R-sig-Geo] Counting number of points per pixel in a grid

[Roger Bivand]

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