Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

This is a crosspost from stackoverflow. I know very little about GIS software and am doing all my mapping in R. Apologies in advance if this is too basic question. Let's say I have two shapefiles from different sources but with different attributes. Say, one is for adminstrative boundaries of Texas (boundaries.shp) and the other is for Texas rivers (rivers.shp). I also have a third file towns.csv showing the locations of towns within the state. After reading in the files, I can overplot town locations on the adimistrative boundaries in the maptools package:

plot(boundaries); points(towns$lon, towns$lat)

But how can I overlay all three? Surely there is an easy way of doing this?

share|improve this question
up vote 7 down vote accepted

PBSMapping should fit your needs. There's a tutorial at NCEAS. The code below is adapted from that tutorial. I'm making assumptions about your data btw. Please edit as appropriate for your situation.


#prepare towns
pts <- read.csv("towns.csv")
towns <- points(towns$lon, towns$lat)
# read in shapefiles 
rivers <- importShapefile("rivers.shp")
boundaries <- importShapefile("boundaries.shp")

# note that importShapefile reads the .prj file if it exists, but it
# does not adopt the proj4 format used by the above approaches
proj.abbr <- attr(boundaries, "projection") # abbreviated projection info
proj.full <- attr(boundaries, "prj") # full projection info
# [1] "LL"

# generate map using PBSmapping plotting functions
plotPolys(boundaries, projection=proj.abbr, border="gray",
    xlab="Longitude", ylab="Latitude")
addPoints(towns, pch=20, cex=0.8)
addLines(rivers, col="blue", lwd=2.0)
share|improve this answer
Thanks, R.K. I'm still trying to get the hang of it, but this was an excellent pointer. – user3671 Feb 10 '12 at 9:22
You're welcome. Have fun mapping :) – R.K. Feb 10 '12 at 23:12

The simplest way to overlay two plots might be using the add = TRUE option in plot. Here is an example with artificial data

# Load sp package for creating artificial data

# Create sample town points
towns <- data.frame(lon = sample(100), lat = sample(100))
towns <- SpatialPoints(towns)

# Create sample polygon grid
grd <- GridTopology(c(1,1), c(10,10), c(10,10))
polys <- as.SpatialPolygons.GridTopology(grd)

# Plot polygons

# Add towns (in red colour)
plot(towns, add = TRUE, col = 'red')
share|improve this answer
Thanks for chipping in, yellowcap. But it says "add is not a graphical parameter." – user3671 Feb 10 '12 at 9:26
The example works fine when I run it on my computer, but it seems that "add" does not work always and depends on the class of the input data, see this post. So my suggestion might not be the best way to go... – yellowcap Feb 10 '12 at 9:58
Not an issue. Thanks for speaking up. – user3671 Feb 10 '12 at 10:16

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.