Dissolve polygons in R

Dissolving polygons is another fairly elementary GIS task that I need to perform regularly. With R this is can be a bit involved, but once done is fully reproducible and the code can be re-used. This post is essentially a companion piece to Clipping polygons in R; I wrote both because I often forget how to complete these tasks in R.

Getting started

Let’s gather together everything we need to complete this example. We need a shapefile of small geographies to ‘dissolve’, a lookup table to tell us which polygons dissolve into which, and we need a couple of R spatial packages to run everything. Let’s get started.

# The default behaviour of this script is to create a folder called 'dissolve-example'
# and download and run everything from here.
# You can change this by modifying the next two lines
dir.create("dissolve-example")
setwd("dissolve-example")

# Load a few packages. dplyr makes merges easier
require("rgdal")
require("rgeos")
require("dplyr")

# Set up shapefile to dissolve. I'm using English regions
download.file("http://census.edina.ac.uk/ukborders/easy_download/prebuilt/shape/England_gor_2011.zip",
              destfile = "lad-region-lookup.zip")
unzip("lad-region-lookup.zip", exdir = ".")
region <- readOGR(".", "England_gor_2011")

# Check the shapefile has loaded correctly
plot(region)
Regions

Regions

# We're going to dissolve all regions in to one country (England!)
# For this we'll create a lookup table and merge with the spatial data
# Hopefully for your 'real' data you have a lookup table of all polygons and
# their larger geography already!
lu <- data.frame()
lu <- rbind(lu, region@data)
lu$CODE <- as.character(lu$CODE)
lu$NAME <- as.character(lu$NAME)
lu$country <- NA
lu$country <- "England"

Merge

We now need to merge the lookup table into our spatial object data frame. We should end up with one row per zone to dissolve, each with a reference for the relevant larger geography. I think the trick is to make sure the row names match exactly, and if you can match the polygon IDs as well with spChFIDs().

# Merge lu (LookUp) into polygons,
region@data$CODE <- as.character(region@data$CODE)
region@data <- full_join(region@data, lu, by = "CODE")

# Tidy merged data
region@data <- select(region@data, -NAME.x)
colnames(region@data) <- c("code", "name", "country")

# Ensure shapefile row.names and polygon IDs are sensible
row.names(region) <- row.names(region@data)
region <- spChFIDs(region, row.names(region))

Dissolve

Now we can do the dissolve. I use gUnaryUnion (and indeed I think unionSpatialPolygons in the maptools package uses this by default).

# Now the dissolve
region <- gUnaryUnion(region, id = region@data$country)

If you want to just plot this using base plot you can stop there. If you want to do anything with the data or plot using ggplot you need to recreate the data frame.

# If you want to recreate an object with a data frame
# make sure row names match
row.names(region) <- as.character(1:length(region))

# Extract the data you want (the larger geography)
lu <- unique(lu$country)
lu <- as.data.frame(lu)
colnames(lu) <- "country"  # your data will probably have more than 1 row!

# And add the data back in
region <- SpatialPolygonsDataFrame(region, lu)

# Check it's all worked
plot(region)

Your plot should look like this:

Dissolved polygons

Example dissolved polygons

And your data frame should look like this:

region@data
##   country
## 1 England
Advertisements

6 thoughts on “Dissolve polygons in R

  1. Greg says:

    This tutorial rules, dude. Thanks! The R package gpclib doesn’t exist anymore and I was pulling my hair out trying to merge polygons.

      • Greg says:

        Hi again! Thanks again for this tutorial. I have an output problem which I hope you might be able to help me with: When I add writeOGR(region, “.”, “region”, driver=”ESRI Shapefile”) to your tutorial code the four output files are: region.dbf, .. .shp, .. .shx, and .. .prj. When I alter the code with my own data the output from the writeOGR command does not include the .prj file. Have you had this problem before? Any idea what could be going wrong?

      • Greg says:

        Nevermind! I figured out the problem. Somewhere in the last week I changed: readOGR() to readShapeSpatial() . I’m not sure why I did that but after changing it back I now get four output files.

      • Glad you got it sorted. As you discovered readShapeSpatial() from maptools doesn’t read in the projection (the .prj file) so if it’s not read in it can’t be written! 🙂

  2. Adam Smith says:

    The `aggregate` function in the `raster` package also handles SpatialPolygons objects:

    For example:
    library(raster)
    region <- readOGR(".", "england_gor_2011")
    region@data$country <- "England"
    england <- aggregate(region, by = "country")
    plot(england)

Comments are closed.