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