NCEAS Scientific Computing: Solutions Center Use Case

Read and Write ESRI Shape Files Using R

The ESRI shape file is a widely used standard data format for storing vector-format geophysical data. R provides several different techniques to import ESRI shape file. This use case demonstrates use of five of the most widely-used packages using point, polygon, and line shapefiles.

The R analytical environment includes a rich set of tools for reading, processing, and displaying spatial data. Most of these tools are compatible with the ESRI Shapefile, which has become an industry standard for storing and exchanging point, line, and polygon geospatial data. The R geospatial tools are provided by an evolving set of R packages, developed to meet the (sometimes overlapping) needs of the scientific community. Within the past two years, the sp package has emerged as a well-defined base class providing a foundation of geospatial object data types for the other spatial classes. Here is a summary of five widely-used shape file-aware R packages:

  • sp contains base classes for the four major spatial data types: point, line, polygon, and grid. sp also includes methods for plotting, printing, segmenting, and summarizing spatial objects. The components of sp can be used individually, or they can and have been used to provide base spatial object services to domain-specific R classes such as PBSmapping (see below).

  • maptools contains tools for input/output, processing, and analysis of geophysical data, in particular, ESRI shape files. maptools also provides software interfaces ('wrappers') for exchanging spatial object with other R packages.

  • PBSMapping was developed by a team of fisheries scientists. It extends the R language to include two-dimensional plotting features similar to those commonly available in a Geographic Information System (GIS).

  • rgdal provides an R interface to the widely-used, open-source Geophysical Data Abstraction Library (GDAL). This package provides services for processing and analysis of vector and raster geophysical data, as well as performing map projections and datum transformations, using the external PROJ.4 library.

  • shape files. provides functions for reading and writing the components of the ESRI Shape File standard. This package is used by R programmers who need to parse and modify the individual shape file components, in particular, the data attribute (.dbf) layer.
Note: These packages' internal representation of the data read from ESRI shape files. may vary. This has implications for use of the data; different functions (as provided by different contributed packages) operate on different types of R objects. In addition, there are some subtle variations in how the data are imported. For example, are character variables in the attribute table converted to factors, or left as-is? Are blank cells in the attribute table represented as NA or as empty characters? Often these subtleties are unimportant or can be dealt with after you read in the data; furthermore, in some cases there exist functions that convert between the different types of R spatial formats. Nevertheless, it pays to be cognizant of these differences. Recent efforts to standardize the foundation-level geophysical data formats (in particular, emergence of sp as the de facto foundation class) will eventually resolve these conflicts and allow researchers to focus their attention on the content of the data rather than on the details of their internal representation.

Examples: Reading Shape files.

The following R code examples demonstrate how to read shape files. using each of the five packages discussed above. In order to run the examples, you may need to install the packages into your R environment using either the R graphical interface packages/install packages option (Microsoft Windows version), or the following R command line entry:

  • > install.packages(c("sp","maptools", "PBSmapping","rgdal", "shape files."))

Note to Linux and Macintosh users: The rgdal package depends on external libraries that must be installed on the local system. For Windows users, these libraries are included in the package binary, and thus no other action is required beyond installing rgdal itself. For Linux and Mac users, it may be necessary to obtain and install both the open-source GDAL and PROJ.4 libraries before installing rgdal.

For more information and instructions, visit http://cran.r-project.org/src/contrib/Descriptions/rgdal.html.

Click here to download a .zip file archive containing the R code and sample shape files. used here.

You can also run these examples by running the R code directly over the internet: Enter the following statement at the R command line; the imported data will be stored as objects in your R environment:

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

I. Using sp and maptools package methods

This section includes the preliminaries: loading the required R packages

Note: The maptools package provides data input method and SpatialPolysDataFrame class; the sp class provides data plotting method 

#
# load the appropriate R package libraries
#
   library(sp)
   library(maptools)
   library(rgdal)
   library(shape files.)
   library(PBSmapping)
#
# First, use methods in the base sp class to read and display
# polygon and point files.
#
   FiveStatePolyFrame = readShapePoly("FiveNWStateAlbersWArea")
# 
#  Review the class of our shape file. object. 
#  
   class(FiveStatePolyFrame)
#
# display the polygons
#
   plot(FiveStatePolyFrame)
#
   print("hit key to see point data in a separate window..")
   browser()
# 
# Lets look at a different data set: Projected polygons and points
# for counties in the five western states, with polygon areas,  
# along with the county polygon centroids.
#
   FiveStateCtyCentr = readShapePoints("FiveNWStateCentrdAlbers")
#
# display the points in a second graphics window.
#
   dev.set(1)
#
   plot(FiveStateCtyCentr)
#
# Next, use the rgdal library method readOGR() to 
#
   print("hit key to try 'rgdal' methods..")
   browser()
#

Here are the shape files.:

Five State County Polygons Five State County Polygons

II. Using rgdal methods

The rgdal package provides more encapsulated methods 

# 
# Next, use methods in the rgdal package to read the same file
# 
   FiveStatePolyOGR = readOGR("FiveNWStateAlbersWArea.shp","FiveNWStateAlbersWArea")
#
# compare the data objects produced by the maptools and rgdal methods
#
   class(FiveStatePolyFrame)
   class(FiveStatePolyOGR)
#
# display the polygons in a new graphics window.
#
   plot(FiveStatePolyOGR)

   print("hit key to view points in a separate window..")
   browser()
   dev.set(1)
#
# read and display the point file.
#
   FiveStateCentrOGR = readOGR("FiveNWStateCentrdAlbers.shp","FiveNWStateCentrdAlbers")
   plot(FiveStateCentrOGR)
#
# Finally, compare the class of the data objects produced by the sp and rgdal methods -
# both return a SpatialPointsDataFrame object.
#
   class(FiveStateCtyCentr)
   class(FiveStateCentrOGR)
#

III. Using shape files. package (shapefile component input)

The shape files. package is intended to read some or all of the shapefile component files into list objects so that the programmer can use their contents in some way. One example of such use would be to read the shapefile attribute table (the .dbf component) so that the attributed could be used in calclations.

Note that the shape files. package does not include data display methods, and that conversion of the lists generated by this package are not easily converted into
spatial data frame objects expected by R packages (e.g., maptools) used to display spatial objects. 

#
# Reading a an ESRI polygon shapefile into a shapefile list using shape files. package:
# 
     FiveStateShapePolys <- read.shapefile("FiveNWStateAlbersWArea")

IV. Using deprecated, but still useful, maptools classes (Map objects)

The first option below also produces a 'SpatialPolygonDataFrame' object. However, because maptools includes some handy functions that operate on data stored as an older 'Map' object type, a second import function exists that will produce a 'Map' object.

#
# Next: use the deprecated Map object classes
# from the maptools package to read the point
# and polygon shape files. and then
# display them with the (plot.Map object). 
#
# deprecated methods still work, but R will display 
# an error message warning user to switch to newer methods. 
#
   MapPolyObj = read.shape("FiveNWStateAlbersWArea")
   plot.Map(MapPolyObj)
#
# do the same thing with the point data set
# display the polygons in a third graphics window.
#
   print("hit key to see points in a new window..")
   browser()
#
   dev.set(1)
#
   MapPointObj = read.shape("FiveNWStateCentrdAlbers")
   plot.Map(MapPointObj)
#
# Finally, lets demonstrate the PBSmapping package.
# PBSmapping methods require spatial data frames be 
# converted to the package's PolySet and EventData 
# objects. Note that here, we use additional method 
# arguments to specify map projection information
# 
#  PolySet and EventData classes.
#
   print("hit key to try 'PBSmapping' methods..")
   browser()

Here are the shape files. Note the addition of X and Y Axes:

Five State County Polygons Five State County Polygons

V. Using PBSMapping package: comprehensive data frame-derived classes.

This approach imports shape files. into a PolySet data object derived from R base data frames, where it is available to the PBSMapping geophysical analysis methods 

#
# finally, display the data using the PBSmapping package
# PBSmapping methods 
#
#  Read the Polygon shape file into a SpatialPolygonsDataFrame, and convert 
#  to a PolySet object compatible with the PBSMapping plotPolys() function. 
#
# Note that we transform the points and polygons into the UTM coordinate system.
#
   CountyPolysSP <- readShapePoly("FiveNWStateAlbersWArea", proj4string=CRS("+proj=utm +zone=9"))
   CountyPolysPS <- SpatialPolygons2PolySet(CountyPolysSP)      
#
# The PBSmapping plotPolys() method adds map bounding box and better and plot axis formatting...
#
   plotPolys(CountyPolysPS,projection=TRUE)
#
# To plot the centroids, we need to transform the SpatialPointsDataFrame to an EventData object
# We create one using the @coords slot from the FiveStateCtyCentr SpatialPointsDataFrame object.
# 
   len = length(FiveStateCtyCentr @data[,1])
   CentroidEventSet <- as.EventData(data.frame(EID=1:len,X=FiveStateCtyCentr@coords[,1],
                                    Y=FiveStateCtyCentr@coords[,2]),projection="utm",zone=9)
#
# add the centroid points to the existing map display.Read
# 
   addPoints(CentroidEventSet)

Here is the plot produced by PBSmapping package:

Five State County Polygons

Learning More:

The CRAN (R) Task View summarizes the R Analytical Environment's spatial data analysis and management tools..

This edition of the R Newsletter contains an in-depth article describing the R spatial data classes; parts of the article appear here.

This ESRI White Paper describes the Shape File format.

Search for information (including PROJ4 string) about your desired Map Projection: EPSG spatialreference.org Web Site

Alternate PROJ4 information: remotesensing.org Geotiff projection Web Site

This Web Site describes the standard .dbf format which is a component of the ESRI Shape File.

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 July, 2007.

Site Home | NCEAS Home | KNB Home