NCEAS Scientific Computing: Polygon Dissolve Using R

Polygon Dissolve Operations using R

This Use Case demonstrates use of Polygon Dissolve operations to aggregate polygons (and their descriptive attributes) into larger polygons. and includes a working end-to-end example R language script.

In this example you will learn how to:

1) Read a set of polygons in ESRI Shape file format into an R SpatialPolygons data object

2) Transform the polygon file using the desired map projection.

3) Dissolve the polygon file's internal polygon boundaries into a set of larger polygons

4) Aggregate descriptive attributes assigned to the input polygons into averages and assign them to the new, larger polygon(s).

Getting started:

The Polygon Dissolve operation removes the interior polygons from a multi-polygon data set. The output is usually a single polygon that defines the outer boundary or 'hull' of the input data set.

The R maptools package provides a set of geospatial data processing and analysis methods, including the unionSpatialPolygons() method for performing polygon dissolve.

We demonstrate polygon dissolve using unionSpatialPolygons() and a map of North American biomes obtained from (web site)

Polygon analysis usually consists of four steps:

1) Import polygon data into R data objects; if desired, display them in a map context
2) Transform all polygon and point files to be used into the same coordinate system and map projection.
3) Aggregate the polygons using the maptools unionSpatialPolygon() function
4) compute the areas of the new, larger polygons, display the areas along with map plots.

We will demonstrate polygon dissolve operations with two data sets: We demonstrate the basic polygon dissolve / attribute aggregation using the counties in North Carolina. Then, we work with the counties of the five Northwestern states (Idaho,Montana, Oregon, Washington, Wyoming) to create outline polygons and polygon areas for the entire region and each state.

The Fundamentals : R Polygon Dissolve operations

Here is a brief example of unionSpatialPolygons() usage:

    InputPolygons = readShapePoly("NCCounties.shp",proj4string= CRS("projection and datum parameters"))
    ProjectedPolygons = spTransform(InputShapefile,CRS(""projection and datum parameters"))
#	
    PolyLabelPts <- getSpPPolygonsLabptSlots(ProjectedPolygons)         # get vector containing the label points (lat/long)
    IDOneBin <- cut(lps[,1], range(lps[,1]), include.lowest=TRUE)       # generate four assignment bins: quartiles of label pt longs
	DissolveResult <- unionSpatialPolygons(ProjectedPolygons ,IDOneBin) # dissolve the polygons into the four bins
    DissolveResultPS <- SpatialPolygons2PolySet(NcDissolve)             # Convert to PolySet for PBSMapping area calculations
    AreasOfResult = calcArea(DissolveResultPS)                          # one simple example: aggregating a polygon attribute

Examples: Aggregating United States County Polygons

Part One: North Carolina Counties
Part Two: Five-State Region

To download a .zip file containing the R script and input files for this example, click here.

Example 1: Aggregating North Carolina County Polygons

back to questions

In this first example, we aggregate the North Carolina county polygons into larger polygons, then compute the areas of the new polygons. We assign the input county polygons to the aggregated polygons using the longitudes of their label points (the approximate center of each polygon). Here are the polygons, followed by the R script that produced them:


North Carolina County Polygons	Aggregated Polygons: Four Regions

Aggregated Polygons: NC State Outline

Calculating North Carolina polygon areas..

North Carolina (100 county polygons) area: 127099 sq km. Hit key to continue

Dissolving North Carolina polygons...

North Carolina (four regions) area: 127099 sq km. Hit key to continue

Dissolving North Carolina polygons (again)...

North Carolina (one region) area: 127099 sq km. Hit key to continue

Here is the R script...

# Polygon Dissolve Demonstration using R geospatial libraries
# Author: Rick Reeves, Scientific Programmer, NCEAS.
# Copyright 2007 National Center for Ecological Analysis 
# and Synthesis
#
# Usage from R command line: 
#
#  1) > setwd("path to folder location of R script and sample files")
#  2) > source("PolyDissolveDemo.r")
#  3) > DemoDriver()
# 
DemoDriver <- function()
{
#   debug(PolyDissolveDemo) # to step through the code, un-comment this line....
   PolyDissolveDemo()
}
PolyDissolveDemo <- function()
{
#
# script demonstrates polygon dissolve using maptools
# and area calculation using PBSmapping package.
#
# first, load the R packages that we will need
#
   library(maptools)   # for geospatial services; also loads foreign and sp
   library(pickling)     # General Polygon Clipping library 
   library(rgdal)      # for map projection work; also loads sp
   library(PBSmapping) # for GIS_like geospatial object manipulation / analysis including poly
# 
# Part 1:
# First example, derived from the maptools unionSpatialPolygon() method documentation in the R 
# maptools user manual: We calculate and sum the area of each polygon in the North Carolina map. 
# Then, use unionSpatialPolygons() to dissolve the county polygon boundaries, and assign each
# polygon's area to one of four regions, based on longitude thresholds.
#
   print("start of demo. Hit key...")
   browser()
   print("reading and transforming North Carolina Shape file")
#
         NorthCaroBase <- readShapePoly("sids.shp",proj4string=CRS("+proj=aea +ellps=GRS80 +datum=WGS84"))
#
# Transform the polygons (which were read in as unprojected geographic coordinates)
# to an Albers Equal Area projection.
#
   NorthCaroProj = spTransform(NorthCaroBase,CRS("+proj=aea +ellps=GRS80 +datum=WGS84"))
#
# Convert to a PolygonSet for comparability with PBSmapping package routines.
#
   NorthCaroProjPS = SpatialPolygons2PolySet(NorthCaroProj)
#
   plotPolys(NorthCaroProjPS, proj = TRUE,col="wheat1",xlab="longitude",ylab="latitude")
#
# polygon area calculations
#
   print("Calculating North Carolina polygon areas...")
   attr(NorthCaroProjPS, "projection") <- "LL"
   NCPolyAreas = calcArea(NorthCaroProjPS,rollup=1)
#
# Compute the area of state by summing the area of the counties.
#
   numCountyPolys = length(NCPolyAreas[,1])
   NCArea = sum(NCPolyAreas[1:numCountyPolys,2])
   print(sprintf("North Carolina (%d county polygons) area: %g sq km. Hit key to continue",
         numCountyPolys,NCArea))
   browser()
#
# create a set of 'breakpoints' that unionSpatialPolygons() method will use to 
# place the dissolved polygons into one of four longitude 'zones' on the output map: 
# Each county/polygon's x/y label coordinates gives the location of that polygon's center.
# 
   print("Dissolving North Carolina polygons...")

   lps <- getSpPPolygonsLabptSlots(NorthCaroProj)
# 
# Assign each county to one of four longitudinal 'bins': 
#
   IDFourBins <- cut(lps[,1], quantile(lps[,1]), include.lowest=TRUE)
#
# Dissolve operations: result is a SpatialPolygons object.
# Convert to PolySet to display new polygons and calculate areas.0
#
  NcDissolve   <- unionSpatialPolygons(NorthCaroProj ,IDFourBins)
  NcDissolvePSFour <- SpatialPolygons2PolySet(NcDissolve)
#
# display in a new window....
#
  dev.set(1) 
  plotPolys(NcDissolvePSFour, proj = TRUE,col="wheat1",xlab="longitude",ylab="latitude")
#
# projection attribute must be "UTM" or "LL" for areas in km ** 2
#
  attr(NcDissolvePSFour, "projection") <- "LL"
  NCDissolvePolyAreas = calcArea(NcDissolvePSFour ,rollup=1)
#
# sum the areas of all polygons in the poly set - 
# the expression 'length(NCDissolvePolyAreas [,1]' returns the number of polygons
#
  NCDissolveArea = sum(NCDissolvePolyAreas [1:length(NCDissolvePolyAreas [,1]),2])
#
  print(sprintf("North Carolina (four regions) area: %g sq km. Hit key to continue",
        NCDissolveArea))
  browser()
#
# next, lets put all county polygons into a single 'bin' to get just a county outline: 
# replace the quantile() call with one to range(), that returns just the min/max of the lps[,1] vector.
#
  print("Dissolving North Carolina polygons (again)...")
  IDOneBin <- cut(lps[,1], range(lps[,1]), include.lowest=TRUE)
  NcDissolve   <- unionSpatialPolygons(NorthCaroProj ,IDOneBin)
  NcDissolvePSOne <- SpatialPolygons2PolySet(NcDissolve)
#
# display in a new window....
#
  dev.set(1) 
  plotPolys(NcDissolvePSOne, proj = TRUE,col="wheat1",xlab="longitude",ylab="latitude")
#
# projection attribute must be "UTM" or "LL" for areas in km ** 2
#
  print("Calculating North Carolina polygon area (again)...")
  attr(NcDissolvePSOne, "projection") <- "LL"
  NCDissolvePolyAreas = calcArea(NcDissolvePSOne ,rollup=1)
  print(sprintf("North Carolina (one region) area: %g sq km. Hit key to continue",NCDissolveArea))
#
  browser()
#
# close all the graphics devices
#
   print("Press any key to see Part 2..")  
   browser()
   while(dev.cur() > 1)
     dev.off()

Example 2: Northwestern United States County Polygons

back to questions

In the second example, we dissolve the county polygons for the five Northwest United States in two stages: First, we aggregate all of the polygons into one regional outline. Then, we generate boundary polygons for each of the five states. To do so, we select the counties for each state, using the state FIPS code attached to each county as a primary key, and then perform the polygon dissolve and area calculations for each state. Here are maps showing the results for the five-state region, and two sample states: Idaho and Montana. The R script that produced these results follows the pictures.


County Polygons: 5 Northwest States	Aggregated Polygons: 5 NW State Region

Idaho Counties	Idaho Boundary with Area

Montana Counties	Montana Boundary with Area

Map with Five State Outines

Reading and transforming Five State County shape file

Calculating Five State County polygon area...

Five State area: 1.28598e+06 sq km. Number of counties: 208. Hit key to continue

State Idaho has 44 counties. Area: 216563 km **2. Hit key to see outline....

State Montana has 56 counties. Area: 384037 km **2. Hit key to see outline....

State Oregon has 36 counties. Area: 251817 km **2. Hit key to see outline....

State Washington has 49 counties. Area: 176539 km **2. Hit key to see outline....

State Wyoming has 23 counties. Area: 257025 km **2. Hit key to see outline....

Here is the R script...

#
# Part 2: Using the county polygons for 5 NOrthwest United States, generate 
# the state boundaries and state areas using the polygons for each state
# the key: use the State FIPS code attribute to select the counties for each state.
#
#   State         FIPS CODE
#   Idaho             16
#   Montana           30
#   Oregon            41
#   Washington        53
#   Wyoming           56
#
# create a list with two columns: state name, FIPS code
#
    StateFIPS <- data.frame(name="StateFIPS",States = c("Idaho","Montana","Oregon","Washington","Wyoming"),FIPSCode = c(16,30,41,53,56))
# 
#  Accessing Data Frame elements: try these combinations:
#
#  StateFIPS[[1]]    = "StateFIPS"
#  StateFIPS[[1]][1] = "StateFIPS"
#  StateFIPS[[2]][1] = "Idaho"
#  StateFIPS[[3]][1] = 16
#  StateFIPS[[2]][2] = "Montana"
#  StateFIPS[[3]][2] = 30
#
   print("Reading and transforming Five State County shape file")
   FiveStateBase <- readShapePoly("FiveNWStatesCounties.shp",proj4string=CRS("+proj=aea +ellps=GRS80 +datum=WGS84"))
   FiveStateBasePS = SpatialPolygons2PolySet(FiveStateBase)
#
# Transform the unprojected input polygons to Albers Equal Area projection.
#
   FiveStateProj = spTransform(FiveStateBase,CRS("+proj=aea +ellps=GRS80 +datum=WGS84"))
#
# Convert to a Polygon set so that PBSmapping routines for area calculation can be used
#
   FiveStateProjPS = SpatialPolygons2PolySet(FiveStateProj)
#
   plotPolys(FiveStateProjPS, proj = TRUE,col="wheat1",xlab="longitude",ylab="latitude")
#
# polygon area calculations
#
   print("Calculating Five State County polygon area...")
   attr(FiveStateProjPS, "projection") <- "LL"
   attr(FiveStateProjPS, "zone") <- 11 # right in middle of the 5 state area
   FiveStatePolyAreas = calcArea(FiveStateProjPS,rollup=1)
#
# Compute the area of state by summing the area of the counties.
#
   numCountyPolys = length(FiveStatePolyAreas[,1])
   FiveStateArea = sum(FiveStatePolyAreas[1:numCountyPolys,2])
   print(sprintf("Five State area: %g sq km. Number of counties: %d. Hit key to continue",
         FiveStateArea,numCountyPolys))
   browser()
#
# create a set of 'breakpoints' that unionSpatialPolygons() method will use to 
# place the dissolved polygons into one of four longitude 'zones' on the output map: 
# Each county/polygon's x/y label coordinates gives the location of that polygon's center.
# 
   lps <- getSpPPolygonsLabptSlots(FiveStateProj)
#
# Generate the five-state boundary polygon.
#
# perform the dissolve, creating the new polygon
#
   print("Dissolving Five State polygons...")
   lps <- getSpPPolygonsLabptSlots(FiveStateProj)
   IDOneBin <- cut(lps[,1], range(lps[,1]), include.lowest=TRUE)
   FiveStateDissolve   <- unionSpatialPolygons(FiveStateProj,IDOneBin)
   FiveStateDissolvePS <- SpatialPolygons2PolySet(FiveStateDissolve)
   FiveStateDissolveAreas = calcArea(FiveStateDissolvePS,rollup=1)
   plot(FiveStateDissolve)
   print(sprintf("Five State outline.Hit key to continue"))
   browser()

#
# start a vector of SpatialPolygons. First element is the 'five state outline'.
# Below, we will append the outlines for each state.
#
   DissolvePolyVec = FiveStateDissolve
#
# Finally: generate the individual state outlines by selecting and then dissolving
# the county polygons for each state. Use FIPS codes to select polygons for each set
# 
# lets create a list of polygons: The first one is the outline of the entire 5-state region,
# followed by the outline of each state.

   print("Creating outline polygons for each state....")
   numberOfStates = length(StateFIPS[[2]]) 
   for (iCtr in 1: numberOfStates) # number of states
   {
      nextStateFIPS = StateFIPS[[3]][iCtr] # List element 3 (FIPS codes), member 'iCtr'  
#
# build logical vector with TRUE values for data frame elements (polygons) with this FIPS code
#
      select = FiveStateProj@data$STATEFIPS == nextStateFIPS
#
# make new DF with counties for this state, then do the dissolve operation.
# to check out structure, check out: str(nextStatePolys)
#
      nextStatePolys = FiveStateProj[select,]     
      nCountiesThisState = length(nextStatePolys@data[,1])
#
# perform the dissolve, create the next state outline polygon
#
      lps <- getSpPPolygonsLabptSlots(nextStatePolys)
      IDOneBin <- cut(lps[,1], range(lps[,1]), include.lowest=TRUE)
      StateDissolve   <- unionSpatialPolygons(nextStatePolys ,IDOneBin)
      StateDissolvePS <- SpatialPolygons2PolySet(StateDissolve)
      plot(nextStatePolys)
      browser()
#
#  calculate state area
#
      attr(StateDissolvePS, "projection") <- "LL"
      attr(StateDissolvePS, "zone") <- 11 # right in middle of the 5 state area
      StateDissolveArea = calcArea(StateDissolvePS,rollup=1)
      print(sprintf("State %s has %d counties. Area: %g km **2. Hit key to see outline...",
      nextStateFIPS = StateFIPS[[2]][iCtr],nCountiesThisState,StateDissolveArea))

#
# last thing: add the latest polygon to the end of the SpatialPolygon vector
#
      DissolvePolyVec = c(DissolvePolyVec,StateDissolve)
# 
# update the area and name attributes of the polygons slot in this new element.
# to determine the slots in this, enter:  > str(Vec[[1]]
#
      attr(DissolvePolyVec [[iCtr+1]]@polygons,"area") = StateDissolveArea
#
      idlabel = sprintf("%s (FIPS Code %d) Area= %g km ** 2",
                        StateFIPS[[2]][iCtr],StateFIPS[[3]][iCtr],StateDissolveArea[2])
      attr(DissolvePolyVec[[iCtr+1]]@polygons,"ID") = idlabel
      plot(DissolvePolyVec [[iCtr+1]])
      title(main = idlabel)
      print(sprintf("State %d of %d. Hit key to see next state...",iCtr,numberOfStates))
      browser()
   }
# 
# Finally, draw a map with the six 'dissolved' polygons in the data frame, using PBSmapping methods
# For this demo, plot the SpatialPolygons in 'DissolvePolyVec list separately.
#
   NextPolyPS = SpatialPolygons2PolySet(DissolvePolyVec[[1]])
   plotPolys(NextPolyPS, proj = TRUE,col="wheat1",xlab="longitude",ylab="latitude")
   print("five state outline (again). Hit key to plot the individual states...")
      browser()
#
# select a color for the first state..
#
   PlotColor = 39
   for (iCtr in 1:numberOfStates+1)
   {
        NextPolyPS = SpatialPolygons2PolySet(DissolvePolyVec[[iCtr]])
        addPolys(NextPolyPS,col=PlotColor)
        PlotColor = PlotColor + 5
   }
   print("End of demo...hit key...")
#
# erase graphics windows
#
   browser()
   while(dev.cur() > 1)
     dev.off()
}

Learning more:

Search for information (including PROJ4 string) about your desired Map Projection: EPSG spatialreference.org Web Site

Alternate PROJ4 information: remotesensing.org Geotiff projection Web Site

UNC/Chapel Hill School of Information and Library Science course (Fall, 2003): Summary of GIS Analysis Functions

Point of Contact for this Use Case: Rick Reeves, NCEAS Scientific Programmer reeves@nceas.ucsb.edu