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.

Answer 1

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))

Answer 2

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 )

Answer 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')