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?

2 Answers 2


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')
  • Thanks for chipping in, yellowcap. But it says "add is not a graphical parameter."
    – user3671
    Feb 10, 2012 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, 2012 at 9:58

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)
  • Thanks, R.K. I'm still trying to get the hang of it, but this was an excellent pointer.
    – user3671
    Feb 10, 2012 at 9:22
  • You're welcome. Have fun mapping :)
    – R.K.
    Feb 10, 2012 at 23:12

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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