NCEAS Scientific Computing: Solutions Center Use Case
Use Case: Creating Maps for Publication using R Graphics
The R software includes a rich set of graphics tools. This case demonstrates how to generate publication-quality maps using these tools.
Scientists often wish to produce maps to include in publications. In the past, they have had to resort to using GIS systems (and the expertise of other GIS specialists) in order to produce these maps.
The R software includes both a rich set of functions for manipulating geospatial data AND a sophisticated graphics toolset. These make it possible for scientists to produce publication-quality maps without resorting to expensive and complicated GIS software. We will review three examples in this Use Case.
In addition to the base R graphics functions, four add-on geospatial packages provide map development capability:
- 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 Function interface to the Geospatial Data Abstraction Library (GDAL), which provides services for importing and exporting vector and raster format geospatial databases into R data objects. GDAL also provides a map projection library .
First Example: Basic Thematic Map
The first map is a basic thematic map displaying Western United States Dams: consisting of point, line, and polygon data layers along with standard map components:map feature annotation, geographic coordinates along the axes, and a legend.
Note: soon we will add two more complex map examples - A satellite image map with vector overlay, and a three-dimensional terrain surface map - to this Use Case.
Here are the typical steps required to produce a map using R:
- 1. Acquire and Ingest the geospatial data layers into the R programming environment.
- 3. Convert the data to the correct R data type.
- 2. If necessary, transform all of the data layers into a common coordinate system and/or map projection.
- 3. Display the spatial data layers in the appropriate order and using complimentary color and pattern attributes.
The sample map is produced by an R function; we call the functions from a 'driver' function as follows:
TestDriver <- function() { MakeWesternUsaDamMap() print("done") }Note: In the near future, additional sample maps will be added to this use case.
Example Map: Western United States Dams
This example produces a simple thematic map showing the location of major dams in the Western United States. The data layers are available from the USGS,and are stored in the industry-standard Shape File format. an industry standard for storing vector geospatial data.
Click here to download a .zip file archive containing all of the R scripts and sample geospatial data files used to produce these maps!
Important Note for UBUNTU users: When running this script, if you receive the error: "Error in X11() : could not find any X11 fonts", you might need to change the default fonts from UTF-8. One way to do so: start your R session from the Linux command line as follows: $ LANG=C R.
Here is the map:
Here is the R script that produced this map:
MakeWesternUsaDamMap<- function() { # library(maptools) # Loads sp package also library(PBSmapping) # Loads sp package also, if not already loaded # # attach to the folder containing map data. NOTE: You may need to modify this path name. # setwd("F:\\Projects\\GeoSpatSemUseCases\\MapProdWithRGraphics\\MapData\\WesternUSA") # # Lets use the PBSMapping package tools to assemble a basic map with polygons, points, lines # # read in, transform, plot the Western State polygons # # Always good to do: save the default graphics/plotting parameters, # as we will be changing a few of them to make these maps.. # DefGraphicsParams <- par() # # Read in the 'base map' polygons: Western United States outlines, # and major Western states reservoirs. # StateOutlines = readShapePoly("WesternStates") # LakeOutlines = readShapePoly("FifteenWesternLakes") # # We would like our map to have labeled X and Y axes and a 'neat line' with tick marks. # One good way to get these is to use the PBSmapping plotPolys routine, which is # enhanced to add these features to the basic polygon plot. # # To specify other graphics parameters, we will use a combination of the inputs to # the PBSmapping functions (e.g., col=), and the par() function, which we use to # globally set the R graphics environment parameters. # # Note: R provides several ways to select graphics colors. Here, we use explicit color names. # To get a (large) list of R graphics color names, run the colors() command. # StatePolysPS = SpatialPolygons2PolySet(StateOutlines) par(cex.lab = 1.1) plotPolys(StatePolysPS,col="wheat1",xlab="Deg W Longitude",ylab="Deg N Latitude") # # Add two sets of drainage: one for California, one for the rest of the country # WestDrainageLines = readShapeLines("WesternUSADrainage") WestPlotLinesPS = SpatialLines2PolySet(WestDrainageLines) addLines(WestPlotLinesPS,col="darkolivegreen3") # CalifDrainageLines = readShapeLines("CentralCaliforniaDrainage") CalifPlotLinesPS = SpatialLines2PolySet(CalifDrainageLines) addLines(CalifPlotLinesPS,col="dodgerblue3") # # Overlay the Reservoir polygons on top of the rivers, using dark green color. # The 'fg' parameter determines the color used to plot polygons. # par(fg="darkolivegreen") plot(LakeOutlines,add=TRUE) # # Add the set of Dam locations. Within the PBSMapping realm, we have to translate # the point locations to an Event Set data object # WestDams = readShapePoints("ThirteenWesternDams") nDams = length(WestDams@data[,1]) DamEventSet = as.EventData(data.frame(EID=1:nDams,X=WestDams@coords[,1],Y=WestDams@coords[,2])) # # Create a 'custom' three-layer symbol for the dams. # Technique: plot the bottom layer (a black 'X') first in bold font (font = 2). # then use the base graphics plot() function to overprint a (larger) orange circle, # and then a (smaller) red circle. graphics parameter 'cex' determines size # of the printing character. # par(font=2) addPoints(DamEventSet,pch=4,cex=2.3,col="black") plot(WestDams,pch=19,cex=1.6,col="darkorange",add=TRUE) plot(WestDams,pch=19,cex=1.0,col="red",add=TRUE) # # Now, label the points with the name of each dam. Get these # from the WestDams SpatialPointsDataFrame object. # damLabels = as.character(WestDams$DAM_NAME) # # use par() to set text size (a bit smaller) and string rotation (5 degrees up) parameters. # par(cex=0.7) par(srt=0.0) # experiment with values from 0.0 - 10.0 to find best 'tradeoff' rotation # text(WestDams,labels=damLabels,pos=4,col="darkslateblue") # # Finally, add the map title and legend. # # First, reset key graphics parameters: Enlarge and de-rotate text, # par(cex = 1.0) par(fg="black") title("Major Dams of the Western United States") # # Set up the (sometimes complex) legend() function parameters # Brief descriptions of the input parameters are supplied in comments. # for more details, consult the help file for legend() >?legend. # par(cex = 0.9) # Shrink text size slightly leg.txt = c("West US Rivers","CA Rivers","Reservoirs","Dams") # These are the legend text strings legend("bottomleft", # placement of legend box inset=0.0125, # offset from map 'neat line' bg="whitesmoke", # legend background color title="Legend", # title legend=leg.txt, # legend element labels lty = c(1,1,-1,-1), # line types for legend. '-1' value indicates that no line is printed pt.cex=2.0, # expand size of legend printed characters pch=c(-1,-1,16,13), # printed character types for legend. Note relationship of '-1' values to those in lty! col = c("darkolivegreen3", # ordered list of colors used to print the legend symbols. "dodgerblue3", "darkolivegreen", "darkred")) # # Restore the graphics environment to the state existant at start of program. # par(DefGraphicsParams) # }Second Example: Raster Satellite Image Base Map with Point and Polygon vector overlays
The second sample map displays a satellite image of the Puget Sound region of Washington state along with the outlines of the corresponding counties and the centroid point for each county. This example demonstrates how to construct a map consisting of multiple, spatially coinciding data layers, in which the base layer is a grid or raster image and the remaining layers are point and polygon vector layers. Note that all of the map layers share the same map projection.
Here are the three base layers:
Raster Grid: Satellite Image Polygons: Counties Points: County Centroids
Here is the final map:
Here is the R code that produced the map:
# # R code will go here # MakeRasterVectorLayerMap <- function() { library(sp) library(rgdal) library(maptools) # # Import three geospatial data layers: # 1) Raster satellite image (GeoTIFF format), used as the base map, # 2) Vector county boundary file (Shape File), the first overlay layer, # 3) Vector county centroid file (Shape File, the second overlay layer. # psImg <- readGDAL("PugetSoundSub1.img") Centroids <- readShapePoints("PSCentroidPointShape") Counties <- readShapePoly("PugetSoundCountiesClp2.shp") # # Create data layers required by the spplot function # Note: in order for polygons to be plotted in the correct sequence # (AFTER the base layer), we need to convert them to SpatialLines data type. # Counties_lines <- as(Counties, "SpatialLines") points <- list("sp.points", Centroids, pch = 21,col="green") polys <- list("sp.lines", Counties_lines, col="white") # # create color scale for the base layer # greys <- grey(0:256 / 256) # # To annotate the county centroids, get the county names # from the polygons, split off the 'County' string from # each name (last word in each name),creating a new list. # Nam = Centroids$COUNTY ns = length(Nam) ss = 1:ns ss = as.character(" ") for (i in 1: ns ) { ss[i] = as.character("") if (!(is.na(Nam[i]))) { subs = strsplit(as.character(Nam[i])," ") nsubs = length(subs[[1]]) - 1 # # Concatenate the name string components. the sub() call removes the leading # the leading " " from the completed string. # for (r in 1: nsubs) { ss[i] = sub(" ","",(paste(ss[i],subs[[1]][r],sep = " "))) } } } len = length(Centroids@coords) LatLongs = Centroids@coords dim(LatLongs) = c(len/2,2) # # arrange for point labels # ptLabels = list("panel.text",LatLongs[,1],LatLongs[,2],labels=ss,font="bold",col="lightblue2",pos=2) grob2 = spplot(psImg, "band1",col.regions = greys, sp.layout = list(points, ptLabels, polys), cuts = length(greys), colorkey = FALSE, scales = list(draw = TRUE), main = "Puget Sound Study Area", sub="Point locations: county centroids", legend = list(right = list(fun = mapLegendGrob(layout.north.arrow())))) plot(grob2) # # TBA: create an actual map legend, add to display list. # print("done") }Learning More:
The CRAN (R) Task View summarizes the R Analytical Environment's spatial data analysis and management tools..
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 edition of the R Newsletter contains an in-depth article describing the R spatial data classes; parts of the article appear in this web page.
This ESRI White Paper describes the Shape File format.
This Web Site describes the standard .dbf format which is a component of the ESRI Shape File.
