[R-sig-Geo] Rasterize polygons

[Robert J. Hijmans]

Tim,

To count the number of points (obsdeep) in each polygon, you can do:

table( over(SpatialPoints(obsdeep), poly) )

There is no role for rasterizing polygons here. But to answer your
question: you are getting a count of 1 because there is never more
than 1 polygon that matches a raster cell (no overlapping polygons).

Robert



On Tue, Jan 7, 2014 at 12:25 PM, Tim Sippel <tsippel at gmail.com> wrote:
> Hi-
> I need to count the number of observations in each of 8 polygons.  I'm
> using the raster package for this, but am having trouble getting what I
> need. Why am I getting a count of 1 in each of the polygons defined below?
>
> Sample data can be downloaded here:
> https://drive.google.com/file/d/0B0d3zfSSPFQsVS04aExhbl9hUU0/edit?usp=sharing
>
>
> library(raster)
>
> #polygons
> reg1<-Polygon(coords=cbind(c(200,220,220,200,200), c(0,0,10,10,0)));
> r1<-Polygons(list(reg1), 'reg1')
> reg2<-Polygon(coords=cbind(c(185,200,200,185,185), c(0,0,10,10,0)));
> r2<-Polygons(list(reg2), 'reg2')
> reg3<-Polygon(coords=cbind(c(200,225,225,200,200), c(10,10,20,20,10)));
> r3<-Polygons(list(reg3), 'reg3')
> reg4<-Polygon(coords=cbind(c(180,200,200,180,180), c(10,10,20,20,10)));
> r4<-Polygons(list(reg4), 'reg4')
> reg5<-Polygon(coords=cbind(c(200,225,225,200,200), c(20,20,30,30,20)));
> r5<-Polygons(list(reg5), 'reg5')
> reg6<-Polygon(coords=cbind(c(180,200,200,180,180), c(20,20,30,30,20)));
> r6<-Polygons(list(reg6), 'reg6')
> reg7<-Polygon(coords=cbind(c(200,235,235,200,200), c(30,30,45,45,30)));
> r7<-Polygons(list(reg7), 'reg7')
> reg8<-Polygon(coords=cbind(c(180,200,200,180,180), c(30,30,45,45,30)));
> r8<-Polygons(list(reg8), 'reg8')
> poly<-SpatialPolygons(Srl=list(r1,r2,r3,r4,r5,r6,r7,r8), pO=1:8)
>
> #raster setup
> xmn=170; xmx=240; ymn=-5; ymx=50
> grid<-raster(nrows=length(seq(ymn,ymx,by=5))-1,
> ncols=length(seq(xmn,xmx,by=5))-1, xmn=xmn, xmx=xmx, ymn=ymn, ymx=ymx); grid
>
> # rasterize & plot
> obsdeep<-rasterize(x=poly, y=grid, field=obs.deep$Longitude1, fun="count")
> plot(obsdeep); plot(poly, add=T)
>
> Many thanks,
>
> Tim
>
>         [[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo