Download data and R Code for this example
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)
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.
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.
# 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")