Convert SpatialGridDataFrame to polygons...?
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 %~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%~%