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:
- Read a set of polygons in ESRI Shape file format into an R SpatialPolygons data object
- Transform the polygon file using the desired map projection.
- Dissolve the polygon file's internal polygon boundaries into a set of larger polygons
- 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 consists of four steps:
- Import polygon data into R data objects; if desired, display them in a map context
- Transform all polygon and point files to be used into the same coordinate system and map projection.
- Aggregate the polygons using the maptools unionSpatialPolygon() function
- 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
Example: unionSpatialPolygons() usage:
InputPolygons = readShapePoly("NCCounties.shp",proj4string= CRS("projection and datum parameters"))
ProjectedPolygons = spTransform(InputShapefile,CRS(""projection and datum parameters"))
# get vector containing the label points (lat/long)
PolyLabelPts <- getSpPPolygonsLabptSlots(ProjectedPolygons)
# generate four assignment bins: quartiles of label pt longs
IDOneBin <- cut(lps[,1], range(lps[,1]), include.lowest=TRUE)
# dissolve the polygons into the four bins
DissolveResult <- unionSpatialPolygons(ProjectedPolygons ,IDOneBin)
# Convert to PolySet for PBSMapping area calculations
DissolveResultPS <- SpatialPolygons2PolySet(NcDissolve)
# one simple example: aggregating a polygon attribute
AreasOfResult = calcArea(DissolveResultPS)
Examples: Aggregating United States County Polygons
To download a .zip file containing the R script and input files for this example, click here.
Example 1: Aggregating North Carolina County Polygons
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 Part 1 R script...
#
# script demonstrares 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(gpclib) # General Polygon Clipping library
library(rgdal) # for map projection work; also loads sp
library(PBSmapping) # for GIS_like geospatial object manipulation / anslysis 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 boundries, 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 Shapefile...")
#
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 compatability 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)
plotPolys(NcDissolvePSFour, proj = TRUE,col="wheat1",xlab="longitude",ylab="latitude")
#
# projecton 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)
plotPolys(NcDissolvePSOne, proj = TRUE,col="wheat1",xlab="longitude",ylab="latitude")
#
# projecton 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 (1 region) area: %g sq km. Hit key to continue",NCDissolveArea))
#
browser()
#
# close all the graphics devices
#
Example 2: Northwestern United States County Polygons
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 Part 2 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 shapefile...")
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 (again)...")
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)
#
# calculate state area
#
attr(StateDissolvePS, "projection") <- "LL"
attr(StateDissolvePS, "zone") <- 11 # right in middle of the 5 state area
StateDissolveArea = calcArea(StateDissolvePS,rollup=1)
nextStateFIPS = StateFIPS[[2]][iCtr]
print(sprintf("State %s has %d counties. Area: %g km **2.",nextStateFIPS,nCountiesThisState,StateDissolveArea[,2]))
print("Press to see state outline.....")
browser()
#
# 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
}
Learning more:
Learning More:
The EPSG spatialreference.org Web SiteSearch for Map Projection information
UNC/Chapel Hill School of Information and Library Science course (Fall, 2003): GIS Analysis Functions
Point of Contact for this Use Case: reeves
nceas [dot] ucsb [dot] edu (Rick Reeves), NCEAS Scientific Programmer
This Use Case was updated June, 2010.
