Compute Polygon Centroids

Download data and R Code for this example

Project Requirement:

The centroid of a polygon is the geometric center of that polygon, calculated using a mathematical algorithm. Polygon centroids are widely used to estimate the 'center of gravity' of individual polygons in a collection. In addition,  centroids are used when plotting maps as the default location for displaying polygon labels.

The R language provides several methods for computing polygon centroids. This use case demonstrates the methods in the PBSmapping package, which includes a collection of easy-to-use geospatial analysis and display functions. Routines from the maptools package provide an interface between PBSmapping data objects and standard R data objects and file formats.

Input Data / Format:

Polygon File: USGS Hydrologic Unit Code (HUC4) Polygons -
ESRI Polygon Shape File (MSA Centroid Points in RED)

Workflow:

1) Read input file (polygons) into an R spatial data object.

2) Calculate polygon centroids for each polygon.

3) Write output file: Point Shape File and CSV file containing HUC4 polygon centroids.

Discussion:

To review and run this example:

1) Download the Zip file archive, unpack into a folder.

2) Be sure that your R (version 2.9 or later) environment includes sp, maptools, rgdal packages.

3) Start R, set working directory to the folder containing the downloaded and unpacked example.

4) At the R command prompt, enter 'source("./ComputePolyCentroids.r")'

5) Enter the command 'ComputePolyCentroids("HucsLower48Simple","NewHuc4Centroids")'

6) When graphic/map is displayed, press any key to complete script execution.

Results:

R Script:


# rgdal package loads sp package

require(rgdal)
require(PBSmapping)
require(maptools)
     
# Read Hydrologic Unit Code polygons, convert to PBSMapping PolySet object.

huc4.polygons <- readOGR(".", "HucsLower48Simple")
huc4.ps <- SpatialPolygons2PolySet(huc4.polygons)

# Calculate the centroids 
 
huc4.centroids <- calcCentroid(huc4.ps, rollup=1)

# Combine the centroid coordinates with the HUC4 codes,
# and promote to a Spatial Points Data Frame

huc4.with.cent <- cbind(huc4.centroids, huc4.polygons)
coordinates(huc4.with.cent) <- c("X", "Y")

# Write the updated HUC cetroids to a CSV file.

write.csv(huc4.with.cent, file="Huc4WithCentroids.csv", row.names=FALSE)

# Alternative: Write an ESRI Point Shape File.

writeSpatialShape(huc4.with.cent, "Huc4WithCentroids.shp")

# Plot the new centroids against the background of the source polygons
# Use the PBSMapping package graphics routines.
 
plotPolys(huc4.ps,col="wheat1",xlab="longitude",ylab="latitude",projection = TRUE)
addPoints(huc4.centroids,col="blue",pch=19)
   
# In Depth: Computing centroids without using the R Spatial data objects

# Read HUC4 polygons directly into a PBSmapping PolySet object

huc4.ps2 <- importShapefile("HucsLower48Simple.shp")

# Calculate the centroids 

huc4.centroids2 <- calcCentroid(huc4.ps2, rollup=1)

# Combine the centroid coordinates with the HUC4 codes, 
# and promote to a Spatial Points Data Frame

huc4.with.cent2 <- merge(huc4.centroids2, attr(huc4.ps2, "PolyData"))
coordinates(huc4.with.cent2) <- c("X", "Y")

# Does this method result in centroids identical to the first method? NO.

all.equal(huc4.centroids, huc4.centroids2)

# Extract PIDs of centroids that differ between the two methods

which.diff <- setdiff(huc4.centroids$PID, merge(huc4.centroids, 
                            huc4.centroids2, sort=FALSE)$PID)

# Plot the discrepancies

plotPolys(subset(huc4.ps, PID %in% which.diff), projection="LL")
addPoints(subset(huc4.centroids, PID %in% which.diff), col="blue")
addPoints(subset(huc4.centroids2, PID %in% which.diff), col="red")

Learning More:

The GIS Primer: The Nature of Spatial Information

Overview: Polygon Centroid Calculations