Combining polygons and calculating their area (i.e. number of cells) in R

Question

I have a simple raster (created with R-package: raster). Using the function "rasterToPolygons" I get polygons of all raster cells that contain the value "1":

library(raster) 
dat = list()
dat$x = seq(1.5, by = 10, len = 10)
dat$y = seq(3.5, by = 10, len = 15)
dat$z = matrix(sample(c(0,1), size = 10*15, replace = T), 10, 15)

r=raster(dat);plot(r)

r_poly = rasterToPolygons(r, fun = function(r) {r == 1}, dissolve = F)
plot(r_poly, add = T)

I do not use "dissolve = T" to avoid that all polygons are merged into one big polygon. Instead, I wish to obtain a new SpatialPolygonsDataFrame in which all polygons that share an edge or a point are combined. Polygons that are clearly separated should be identifiable as individual polygones. Based on the new SpatialPolygonsDataFrame I would like to analyze the size of the combined polygones as follows:

b = extract(r,r_poly_new) # "r_poly_new" contains the combined polygons
str(b)                    # list of clearly separated polygons
tab = lapply(b,table)      
tab

My question is twofold: 1) How to combine polygones that share an edge or point? 2) How to get this information into a format which allows analyzing the areas of the combined polygones? Thanks very much for your feedback.

Check out `rgeos::gUnaryUnion` and `rgeos::gUnionCascaded` - I think one of these is what you are after. — Simon O'Hanlon, Dec 18 '13 at 13:19

Josh O'Brien · Accepted Answer · 2014-04-18T02:10:00.127

10

You could first use raster::clump() to identify clusters of connected raster cells and then apply rasterToPolygons() to "polygonize" those cells. (Do note, though, that each clump's area can be computed directly from the RasterLayer without converting it to a SpatialPolygonsDataFrame, as shown below):

library(rgeos) ## For the function gArea

## Clump and polygonize
Rclus <- clump(r)
SPclus <- rasterToPolygons(Rclus, dissolve=TRUE)

## Check that this works
plot(SPclus, col = seq_along(SPclus))

## Get cluster areas from RasterLayer object
transform(data.frame(freq(Rclus)), 
          area = count*prod(res(Rclus)))

## Get cluster areas from SpatialPolygons object
transform(data.frame(SPclus), 
          area = gArea(SPclus, byid=TRUE))

enter image description here

edited Apr 18 '14 at 02:10

answered Dec 18 '13 at 16:52

Josh O'Brien

159,210
26
366
455

@SimonO101 -- The take-home message may be that **raster** and **rgeos** are *both* exceptionally comprehensive and smoothly functioning R packages. – Josh O'Brien Dec 18 '13 at 20:21
Absolutely agree. Add `sp` to that list. The triumvirate that rival ArcGIS for features and smash it into the ground for usability. – Simon O'Hanlon Dec 18 '13 at 21:10
@JoshO'Brien -- Thanks a lot! This works great and is exactly what I was looking for! However, I discovered that the last part of the R-code (i.e Get cluster areas from SpatialPolygons object) was not always working properly. Interestingly, (sometimes) the area of the largest cluster from the SpatialPolygons object was overestimated, whereas the cluster areas from RasterLayer object were always estimated correctly (I counted several times ;) ). Again, many thanks for the help! Alex – Alex Z Dec 18 '13 at 23:16
I am continuously surprised at how powerful `raster` is. – natario Jun 06 '16 at 11:57

score 3 · Answer 2 · answered Dec 18 '13 at 13:24

3

The rgeos package has many polygon manipulation tools. gUnion will union together touching polygons:

require(rgeos)
uni <- gUnion( r_poly , r_poly )
plot( uni , col = 2 )

enter image description here

answered Dec 18 '13 at 13:24

Simon O'Hanlon

58,647
14
142
184

Thanks a lot for this quick reply. However, the problem is that I still get one multipart polygon. What I need are singlepart polygons (in your example: 5 individual polygons). Any suggestions? Many thanks in advance! – Alex Z Dec 18 '13 at 13:38
1

@AlexZ yes there is a way to do this. It currently escapes me but I am looking. `rgeos` is very powerful and contains many polygon manipulation tools. Also asking this question on the `r-sig-geo` mailing list will undoubtedly raise a response from those with more know-how than me. – Simon O'Hanlon Dec 18 '13 at 13:55
Dear Simon, Thanks for your reply. I will ask the question on the r-sig-geo mailing list. Best regards, Alex. – Alex Z Dec 18 '13 at 15:03

score 3 · Answer 3 · edited Apr 18 '14 at 02:19

3

rasterToPolygons() is a computationally very expensive operation so, assuming the CRS is planar, I would go for:

m <- clump(r)
f <- freq(m)
f[,2] <- f[,2] * xres(r) * yres(r)

For lon/lat, I would use:

a <- area(r)
zonal(a, m, 'sum')

edited Apr 18 '14 at 02:19

Josh O'Brien

159,210
26
366
455

answered Dec 18 '13 at 17:41

Robert Hijmans

40,301
4
55
63

Combining polygons and calculating their area (i.e. number of cells) in R

3 Answers3

Linked