[R-sig-Geo] overlay points and polygons and show maps

[Roger Bivand]

On Sat, 19 Jan 2008, Takatsugu Kobayashi wrote:

(Please note that the PBSmapping maintainer and authors are not subscribed 
to this list, so several of your questions about the use of PBSmapping 
functions have not been answered here; I have added a CC to the 
maintainer).

> Hi,
>
> I am trying to show grid cells with colors. Colors depends on the number
> of points falling within grid cells. Here is what I did;
>
>
> rm(list=ls(all=T))
> library(sp)
> library(maptools)
> library(PBSmapping)
>
> ### Define parameters for a grid
> start.point <- -5
> cco <- c(start.point,start.point)
> csize <- c(1,1)
> n <- abs(cco[1]*2)+1
> cnum <- c(n,n)
> grd <- GridTopology(cco, csize, cnum)
> r.grd <- as(grd, "SpatialPolygons")
> r.grd.ps <- SpatialPolygons2PolySet(r.grd)
>
> ### Plot random points
> x <- rnorm(10000)
> y <- rnorm(10000)
> events <- as.EventData(data.frame(EID = 1:10000, X=x, Y=y))
>
> ### Count # of points within cells
> n_in_cells <- findPolys(events, r.grd.ps)
> n_cells <- aggregate(n_in_cells, list(n_in_cells[,2]), sum)[,c(1,4)]
>
> Then I would like to add n_cells to r.grd.ps, and show maps with colored
> cells using plotPolys(). Any suggestions?

If an sp solution is of interest:

library(sp)
start.point <- -5
cco <- c(start.point,start.point)
csize <- c(1,1)
n <- abs(cco[1]*2)+1
cnum <- c(n,n)
grd <- GridTopology(cco, csize, cnum)
r.grd <- as(grd, "SpatialPolygons")

set.seed(1)
x <- rnorm(10000)
y <- rnorm(10000)
xy <- SpatialPoints(cbind(x, y))

o1 <- overlay(xy, r.grd)
res <- numeric(prod(cnum))
tres <- table(o1)
ntres <- as.integer(names(tres))
res[ntres] <- tres
IDs <- sapply(slot(r.grd, "polygons"), function(x) slot(x, "ID"))
resdf <- data.frame(cnt=res, row.names=IDs)
r.grd.df <- SpatialPolygonsDataFrame(r.grd, resdf)
spts <- list("sp.points", xy, pch=".", col="black")
spplot(r.grd.df, "cnt", sp.layout=spts)

or

g.grd <- SpatialGrid(grd)
o2 <- overlay(g.grd, xy)
res <- numeric(prod(cnum))
tres <- table(o1)
ntres <- as.integer(names(tres))
res[ntres] <- tres
g.grd.df <- SpatialGridDataFrame(grd, data=data.frame(cnt=res))
spts <- list("sp.points", xy, pch=".", col="black")
spplot(g.grd.df, "cnt", sp.layout=spts)

avoiding casting to polygons, and side-stepping the ID key needed to 
construct a SpatialPolygonsDataFrame.

Hope this helps,

Roger


>
> Thank you very much
>
> tk
>
> _______________________________________________
> 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