R help - Mar 2006 - Getting polygons from contiguous grid cells having same value in sp

I have been using the sp package to deal with gridded data.  I would
like to know how to extract the vertices (x,y) of polygons that outline
areas of like-valued cells in a grid.  Here is a simple 3x3 grid:

2  2  1
1  2  1
1  1  3

x <- c(1,1,1,2,2,2,3,3,3) # define a 3 x 3 array of points
y <- c(1,2,3,1,2,3,1,2,3)
h <- SpatialPoints(cbind(x,y))  # make these an sp object of the points
class
h$z<- c(1,1,2,1,2,2,1,1,3) # add field values
gridded(h) <- T  # make it a grid

How could I get the vertices for the region of value=2 cells?  Here is
what they would be, in clockwise order from the upper left corner:
(0,3) (2,3) (2,1) (1,1) (1,2) (0,2)

Thanks,
Scott Waichler
Pacific Northwest National Laboratory
scott.waichler _at_ pnl.gov

On Tue, 21 Mar 2006, Waichler, Scott R wrote:
> 
> I have been using the sp package to deal with gridded data.  I would
> like to know how to extract the vertices (x,y) of polygons that outline
> areas of like-valued cells in a grid.  Here is a simple 3x3 grid:
> 
> 2  2  1
> 1  2  1
> 1  1  3
> 
> x <- c(1,1,1,2,2,2,3,3,3) # define a 3 x 3 array of points
> y <- c(1,2,3,1,2,3,1,2,3)
> h <- SpatialPoints(cbind(x,y))  # make these an sp object of the points
> class
> h$z<- c(1,1,2,1,2,2,1,1,3) # add field values
> gridded(h) <- T  # make it a grid
> 
> How could I get the vertices for the region of value=2 cells?  Here is
> what they would be, in clockwise order from the upper left corner:
> (0,3) (2,3) (2,1) (1,1) (1,2) (0,2)
Perhaps R-sig-geo would be a better place to ask (it would):

x <- c(1,1,1,2,2,2,3,3,3)
y <- c(1,2,3,1,2,3,1,2,3)
h <- SpatialPoints(cbind(x,y))
z <- c(1,1,2,1,2,2,1,1,3)
h1 <- SpatialPixelsDataFrame(h, data=list(z=z))
h2 <- as.SpatialPolygons.SpatialPixels(h1)
plot(h2)
h3 <- h2[h1$z == 2]
plot(h3, add=TRUE, col="grey")
library(spgpc) # install from the sourceforge repository, CRAN Spatial 
	       # task view link
h4 <- unionSpatialPolygons(h3, rep("z2", sum(z == 2)))
plot(h4, add=TRUE, border="red", lwd=2)

but note that you didn't allow for the coordinates being in the cell 
centres, so you are off by 0.5:
> slot(slot(slot(h4, "polygons")[[1]], "Polygons")[[1]],
"coords")
     [,1] [,2]
[1,]  2.5  2.5
[2,]  2.5  1.5
[3,]  1.5  1.5
[4,]  1.5  2.5
[5,]  0.5  2.5
[6,]  0.5  3.5
[7,]  2.5  3.5
[8,]  2.5  2.5

and because gpclib isn't checking to see if line segments are straight 
lines, there remain vertices at (2.5,2.5) and (1.5,3.5) that are strictly 
unneeded.

Next time R-sig-geo?

Roger

> 
> Thanks,
> Scott Waichler
> Pacific Northwest National Laboratory
> scott.waichler _at_ pnl.gov
> 
> ______________________________________________
> R-help at stat.math.ethz.ch mailing list
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide!
http://www.R-project.org/posting-guide.html
> 
-- 
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