[R-sig-Geo] Convert SpatialGridDataFrame to polygons...?

[Gavin Simpson]

On Thu, 2010-07-15 at 21:06 +0200, Roger Bivand wrote:
> On Thu, 15 Jul 2010, Crowe, Andrew wrote:
> 
> > Gavin
> >
> > Gavin
> >
> > You can convert the raster grid to a polygon grid by coercing the 
> > SpatialGridDataFrame to a SpatialPixelsDataFrame then from that to a 
> > SpatialPolygonsDataFrame.
> >
> > You should then be able to recode the polygon data for each particular 
> > depth class and union them using unionSpatialPolygons.  This then needs 
> > joining with a data frame created from the depth classes to create a 
> > SpatialPolygonsDataFrame for output to a shapefile using writePolyShape

A little belatedly as I didn't get chance to work through the
suggestions from Andrew and Roger before I had to leave for an extended
fieldwork trip that I have just returned from...

Thanks Roger and Andrew for your suggestions. Roger asked me to let you
know how I got on, so here goes.

Roger's solution was what I had in mind when I was thinking how to do
things. And it worked nicely for some of the sites I needed to apply
this to. But for most sites it didn't work because many of the contour
lines produced from the grid were open lines, not closed lines. I view
this as a feature of the data I had to do this for (lake bathymetry data
where we don't have an accurate lake outline (0m depth) or indeed for
some sites even data over the entire lake).

I next tried Andrew's suggestion and this worked for all the sites I
needed to apply it to as it basically turns the grid into a series of
polygons and we then union the cell polygons that belong to same depth
interval. The only downside is that, as Roger mentioned, it produces
"jagged" sets of polygons, and these polygons are quite complex; some
depth classes are represented by numerous individual polygons
representing single pixels in the original grid. C'est la vie.

Finally, I've begun using the sp stack and related packages more and
more recently and I now have a reasonable working feel for what can be
done and how to achieve some fairly complex (for me, anyway) results. I
have been struck by how powerful and easy to use all this is. This is
due to the authors and maintainers of the various spatial R packages. As
it is easy to overlook this fact when deadlines loom etc; here's a BIG
thank you to you all, from me, for all the hard work!

G

> In addition, you can massage the output of contourLines() to make a less 
> jagged version:
> 
> library(maptools)
> data(volcano)
> x2 <- ContourLines2SLDF(contourLines(volcano))
> # volcano is a regular "image" style object, SpatialGridDataFrame
> # objects can be coerced with as.image.SpatialGridDataFrame()
> # we don't (yet) have contourLines() for sp objects
> x3 <- x2[(as.character(x2$level) > "140" & as.character(x2$level) <
>   "170"),]
> cx3 <- coordinates(x3)
> str(cx3)
> crds <- rbind(cx3[[1]][[1]], c(NA, NA), cx3[[2]][[1]], c(NA, NA),
>    cx3[[2]][[2]], c(NA, NA), cx3[[1]][[2]])
> # this could be automated in a loop
> colnames(crds) <- c("x", "y")
> pols <- map2SpatialPolygons(as.data.frame(crds), IDs=rep("1", 4))
> # number of IDs like input polygons
> pls <- slot(pols, "polygons")
> gpclibPermit()
> pls1 <- lapply(pls, checkPolygonsHoles)
> slot(pols, "polygons") <- pls1
> plot(pols, col="grey", pbg="yellow")
> 
> Both the use of checkPolygonsHoles() and unionSpatialPolygons require the 
> use of gpclib with its bad licence. We're at mid-term with the GSoC rgeos 
> project now, which provides similar facilities.
> 
> Using the contour lines may give a visually more pleasing, and maybe more 
> helpful rendering than the blocky polygonised raster cells (which anyway 
> will have been interpolated).
> 
> Please let us know how you get on,
> 
> Roger
> 
> 
> >
> > Hope that helps
> >
> > Andrew
> >
> > Dr Andrew Crowe
> >
> > Lancaster Environment Centre
> > Lancaster University
> > Lancaster    LA1 4YQ
> > UK
> >
> > Tel: +44 (0)1524 595879
> >
> > ________________________________
> >
> > From: r-sig-geo-bounces at stat.math.ethz.ch on behalf of Gavin Simpson
> > Sent: Thu 15/07/2010 12:28 PM
> > To: r-sig-geo at stat.math.ethz.ch
> > Subject: [R-sig-Geo] Convert SpatialGridDataFrame to polygons...?
> >
> >
> >
> > Dear list,
> >
> > I have been given a set of ASCII grids which I read into R using
> > readAsciiGrid(). The grids represent an interpolated water depth data
> > set for a set of lakes; each grid represents a separate lake. For some
> > reason I have been asked to convert this gridded data into a set of
> > polygons. For example; given the gridded data, I want to generate a
> > polygon that covers the area of the lake (from the grid) that is 1-2m
> > deep say. I would then need to write out this polygon as an ESRI shape
> > file.
> >
> > Are there tools in R to do this sort of conversion?
> >
> > It is a bit more complex that the simple example, as lakes can have
> > multiple basins (areas of deeper water). Say a lake has two deep areas
> > so a cross section through the lakes might look like a W, for such a
> > lake I would require an object that contained two polygons representing
> > the area of the lake between 4 and 5 meters say in *each* of the basins.
> > These would then be read out into a shape file. My GIS colleague
> > indicates that ArcGIS/ArcMAP has a function to polygonise a raster grid
> > in such a way. Is there something similar in the R spatial software
> > stack?
> >
> > An alternative approach would be to contour (at appropriate depth
> > intervals) the interpolated grid and convert the set of contours into
> > polygons for individual lakes. Is there a function to do this?
> >
> > As you can see, this is very much a look-see at the moment. We are
> > evaluating whether it would be quicker to work out how to do this in R
> > as opposed to doing it by hand in ArcGIS.
> >
> > Any suggestions would be most gratefully received indeed!
> >
> > All the best,
> >
> > Gavin
> > --
> > %~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%
> > Dr. Gavin Simpson             [t] +44 (0)20 7679 0522
> > ECRC, UCL Geography,          [f] +44 (0)20 7679 0565
> > Pearson Building,             [e] gavin.simpsonATNOSPAMucl.ac.uk
> > Gower Street, London          [w] http://www.ucl.ac.uk/~ucfagls/
> > UK. WC1E 6BT.                 [w] http://www.freshwaters.org.uk <http://www.freshwaters.org.uk/>
> > %~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%
> >
> > _______________________________________________
> > R-sig-Geo mailing list
> > R-sig-Geo at stat.math.ethz.ch
> > https://stat.ethz.ch/mailman/listinfo/r-sig-geo
> >
> > _______________________________________________
> > R-sig-Geo mailing list
> > R-sig-Geo at stat.math.ethz.ch
> > https://stat.ethz.ch/mailman/listinfo/r-sig-geo
> >
> 

-- 
%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%
 Dr. Gavin Simpson             [t] +44 (0)20 7679 0522
 ECRC, UCL Geography,          [f] +44 (0)20 7679 0565
 Pearson Building,             [e] gavin.simpsonATNOSPAMucl.ac.uk
 Gower Street, London          [w] http://www.ucl.ac.uk/~ucfagls/
 UK. WC1E 6BT.                 [w] http://www.freshwaters.org.uk
%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%