NCEAS Scientific Computing: Solutions Center Use Case

Generate Convex Hull polygons from ESRI Shape file Point Data using R

This case presents an R script that: 1) Reads and displays in an outline map ESRI point and polygon Shape Files for a three-county Central California region; 2) demonstrates two different R programming methods for generating the Convex Hulls bounding the polygons 'belonging' to each county, 3) calculates the areas of each county and Convex Hull polygon.

Getting Started:

A Convex Hull polygon is defined as the smallest convex polygon bounding all of the members of a point set. Ecologists often calculate convex hull polygons from point data measurements (such as species foraging data) as estimates of feeding ranges and other population attributes.

As with many R programming tasks, there is more than one way (and associated set of R packages) to calculate Convex Hulls. In this use case, we will demonstrate two methods: from the sp and maptools packages to read and convert the incoming shape files into plot-friendly Polygon data objects. Routines from the splancs package determine the polygons containing each point (the well-defined 'point-in-polygon' problem).

This example uses the following R packages:

• The base grDevices package provides base graphics services to other R packages, including point, line, and polygon drawing and color and pattern control.
  
  
• PBSMapping provides two-dimensional spatial analysis, transformation, and plotting functions typically found in Geographic Information Systems. This package was developed by scientists conducting statistical fisheries research at the Pacific Biological Station in Nanaimo, British Columbia.


• splancs provides functions for spatial and space-time point pattern analysis of polygons.


• sp is a critical base package for spatial operations in R. It provides definitions for basic spatial classes (points, lines, polygons, pixels, and grids) in an attempt to unify the way R packages represent and manage these sorts of data. It also includes several core functions for creating and manipulating these data structures. The sp package is automatically loaded when the splancs package is used. 

If you have not already added these packages to your local R installation, you can do so using the Packages/Install Packages option in the R GUI toolbar (Ms Windows XP users), or using the following statement at the R command line:

install.packages(c("maptools","splancs","PBSmapping")) # sp loaded automatically with maptools or splancs

For tips on installing R on the Ubuntu Linux platform widely used at NCEAS visit the link:

http://help.nceas.ucsb.edu/index.php/R

To run the sample code discussed here, download and expand this .zip archive: (ConvexHullRDemoFiles.zip-20 kb). This file includes the R script and input files. After starting an R session from the same directory in which you have expanded this archive, you can run the example R code directly over the internet: Simply enter the following statement at the R command line; the imported data will then be stored as objects in your R environment:

source("http://www.nceas.ucsb.edu/scicomp/GISSeminar/UseCases/PointPolyRegionsOfInterest/PointInPolyProc.r")

You may also enter the statements one at a time by typing or pasting them into the R command line. To view the online help page for any R method, be sure the appropriate library has been loaded and then enter the command name proceeded by a question mark:

?coordinates
?inout

Working through the example program: CalculateConvexHull.r

The ConvexHullDemo() function:

1) Reads and converts an ESRI polygon Shape File and an ASCII .csv format point location file into graphics-compatible formats,

2) Assigns the point locations to their bounding polygons using a Point-in-Polygon algorithm

3) Generates convex hull polygons around several groups of points

4) Displays the input and calculated spatial features (e.g., convex hulls) using a simple map graphic.

ConvexHullDemo() demonstrates two different R programming approaches for geospatial data processing and display: routines from the maptools library (combined with methods from the base grDisplay graphics library , and routines from the PBSmapping package.

First, configure the R libraries, read the spatial data files into R geospatial data objects, and plot them using base R graphics tools

ConvexHullDemo <- function()
{
#
# R Programming convex hull demonstration script
#
# Load the required libraries:
# 
   library(maptools)   
   library(PBSmapping) 
   library(splancs)
#
# When called, debug function enables single-step through the ConvexHullDemo
# function using 'next' command
#
#  debug(ConvexHullDemo)
#   browser()
#
# read in two spatial data sets: Three California county polygons
#
   TriCounty <-read.shape("TriCountyPolyDD.shp")    
   sites <- read.csv("ThreeCountyPoints.csv") 
# 
# Create a 1-dim matrix of unique county keys
# FIPS code is the unique key for counties 
# Note for later: The first key,"079", corresponds to SLO County
#
   indexes = as.matrix(TriCounty$att.data["COUNTYFIPS"])
#
# get number of points in the 'sites' data set
# 
   len = length(sites$LatDD) 
# 
# promote the point data frame to an SpatialPointsDataFrame object 
# so that we can add new attributes to the data frame.
# 
   sites2 = sites
   coordinates(sites2) = c("LonDD","LatDD") 
#
# get the lat/lon coordinates for the points
#
   xy=coordinates(sites2)
#
# append an attribute column to the point SPDF to hold the covering polygon
#   
   ar = array(0,len)
   sites2$inpoly = as.vector(ar)  
   TcPolys = Map2poly(TriCounty)
# 
# Lets have a look at the data: plot all the polygons and then all the points 
# 
   plot(TcPolys)
   points(xy,pch="*") 

Here are the baseline counties and their point locations

Next, compute Point-in-Polygon, convex hull, and convex hull area for San Luis Obispo County:

#
# Next, perform Point-in-Polygon and convex hull calculations: 
#
# Loop through all of the individual polygons in the Counties shape file. 
#

# First, we demonstrate the chull() method 
# using the single San Luis Obispo County polygon (and the points that it bounds).
# The San Luis Obispo County points have the greatest spatial extent. SLO county
# is the first polygon in TcPolys, so the remaining 'subsites' object contains the 
# correct points. 
#
   ThisPoly = TcPolys[[1]] 

   inside = inout(xy,ThisPoly)
#
# assign the FIPS code for current polygon to the attribute column of all points falling within.
# 
   sites2$inpoly[inside] = indexes[1]  # SLO County key = "079" 
#
# create a two-column matrix containing coordinates of sites within current polygon
#
  subsites = as.matrix(cbind(xy[inside==TRUE,1],xy[inside==TRUE,2]))
#
   hullpts = chull(subsites)
# 
# chull() returns a vector of integer coordinates. So, for plotting purposes,
# close the 'polygon' by appending the first to the last.
#
   hullpts = c(hullpts,hullpts[1])
#
# plot the hull by referencing the 'subset of the subset' points and connecting the dots.
#
   lines(subsites[hullpts, ],col="green")
#
   matChull = as.matrix(subsites[hullpts, ])
#

 Here is the result of point-in-polygon and convex hull calculations for San Luis Obispo County, in which the point data exhibit the greatest spatial dispersion.


Then, compute the convex hull area:

#
# Next, compute the area of the convex hull 
# This requires some preprocessing as we need to construct
# a polygon whose vertices are geographic coordinates and
# which is projected into UTM Coordinates. 
# To do so, we use methods and data objects from the 
# PBSMapping package. 
#
# First,create a matrix containing the geographic coordinates
# of the point coverage elements that make up the convex hull.
#
# Create PBSMapping data types: Polygon and Polygons (list of Polygon)
#
   poly = Polygon(matChull)
   ID = c("P1") # character vector label required by Polygons data type
# 
# we need to create a PolySet object in order to use the calcArea method.. 
# Remember that 'mat'is an R 2-D matrix with the closed coordinates
# (last point = first point) for a single polygon
#
   ps = appendPolys(NULL,matChull,1,1,FALSE)
#
# to compute area in units of kilometers ** 2, set the projection
# attribute of the PolySet to UTM and Zone 9.
#
   attr(ps,"projection") = "LL"
   attr(ps,"zone") = 9
   psUTM = convUL(ps)
#
   polygonArea = calcArea(psUTM,rollup=1)
#
# polygonArea is a compound data type, see them with the command attributes(polygonArea)
# we will print just the $area attribute.
#
   print(sprintf("SLO County convex hull area (using chull()) is %f km**2.",polygonArea$area))

Here is the Convex Hull for SLO County:

Finally, perform the same set of operations using methods from the PBSmapping package:

  print("hit key to see results of PBSmapping package-based solution......")
   browser()
#
# Now, calculate the Convex Hull for the same county using only
# the spatial data processing routines in PBSmapping. 
# Note that we used some PBSmapping routines in the first example. 
#
# remember, 'sites' is the data frame containing the three-county point data
#
# To compute the convex hull, convert the subsites data frame (with SLO county 'inpoly'
# field marking the points of interest) to a PBSmapping-format 
# EventData object with longitude and latitude columns named X and Y, respectively.
#
   SloSites = sites2[sites2$inpoly == indexes[1],]
#
# Lets now use the  the PBSMapping methods for reading and displaying ESRI Shapefiles
#
# Often, it is necessary to use a new data object/format when using a different R geospatial package. 
# In this case, the PBSMapping package uses the EventData and PolySet data objects to 
# handle point and polygon data, respectively
#
# First, store the point data in EventSet objects.
#
# The calcConvexHull() requires an EventData (data frame) objectleng
# We create one using the @coords slot from the SloSites SpatialPointsDataFrame object.
# To see the underlying structure of SloSites, use '> Str(SloSites)'
#
   len = length(sites2@data[,1])
   AllCountyEventSet <- as.EventData(data.frame(EID=1:len,X=sites2@coords[,1],Y=sites2@coords[,2]))
   len = length(SloSites@data[,1])
   SloEventSet <- as.EventData(data.frame(EID=1:len,X=SloSites@coords[,1],Y=SloSites@coords[,2]),projection="LL",zone=9)
#
# compute the convex hull for SLO county points; results stored in a PolySet data object.
#
   convexHullPoly = calcConvexHull(SloEventSet)
#
# set the projection attribues to: Latitude/Longitude, UTM Zone 9
#
   attr(convexHullPoly ,"projection") = "LL"
   attr(convexHullPoly ,"zone") = 9
#
#  Read the Polygon shape file into a SpatialPolygonsDataFrame, and convert this 
#  to a PolySet object compatible with the PBSMapping plotPolys function. 
#  Set projection attributes to match the convex hull polygon.
#
   CountyPolysSP <- readShapePoly("TriCountyPolyDD", proj4string=CRS("+proj=longlat +zone=9"))
   CountyPolysPS <- SpatialPolygons2PolySet(CountyPolysSP)
#
# Open a second graphics device: dev.set(1) opens the first available new device
#
   dev.set(1)
#
#  Plot the County polygons. Note: the data to be displayed has been assigned a projection
#  (rather than just plotted as x/y coordinates), so the resulting map has a slightly
#  different appearence. 
#
   plotPolys(CountyPolysPS,projection=TRUE)
#
# plot only the SLO points in the eventSet, using addPoints()
# 
   addPoints(SloEventSet)
#
# Next, draw the latest version of the convex hull on the new map. 
#
   addPolys(convexHullPoly,border = "red")
#
# Finally, compute the area of the convex hull. 
# First, convert the coordinates to UTM so that area calculations are in units of km ** 2.
#
   convexHullPolyUTM = convUL(convexHullPoly)
   polygonArea = calcArea(convexHullPolyUTM)
#
  print(sprintf("SLO County convex hull area (using calcConvexHull()) is %f km**2.",polygonArea$area))
  print("end of example")
#
# Done
#
		

Here are the results of the different Convex Hull calculations: First, as generated using grDevices/chull()

SLO County convex hull area (using chull()) is 6596.811490 km**2.

 

Next, the convex hull as generated by PBSmapping/CalcConvexHull()

SLO County convex hull area (using CalcConvexHull()) is 6596.811490 km**2.

Learning More:

The R Programming Language / CRAN Task View: Analysis of Spatial Data

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

Site Home | NCEAS Home | KNB Home