Trying to limit my shapemap to state boundary lines in R

Question

I have a feeling I'm either trying to do this incorrectly, or I'm just not searching for the correct terms, but I cannot figure out how to do what I want to do.

I am trying to plot the ranges of fish species (shapefile #1) in Maryland/Delaware/Virginia superimposed over the state maps, then add major rivers (shapefile #2) in R.

I have no trouble getting everything plotted, but since both shapefiles are for the entire United States, they plot polygons/lines on my maps outside the MD/DE/VA region that I'm interested in.

Here's what I've done so far in R:

> library(maps)
> library(maptools)
> library(scales)
> mudsunfish=readShapePoly("~/Desktop/Maryland/Acantharchus_pomotis.shp")
> map('state', regions=c('maryland', 'virginia', 'delaware'), lwd=3)
> plot(mudsunfish, add=TRUE, col=alpha("skyblue", 0.3), border=FALSE)

This gives me the below image (not worried about the rivers in this instance):

Notice the blue shading in the upper right, where New Jersey and Pennsylvania would be. I don't want that. I can remove the shading in Photoshop, but with >120 species and rivers to add, I'd like to not have to do that.

How can I limit my shapefiles to the borders of the states? Should I be going about this with different packages? I'm sure someone else has run into this before, but I'm at a loss for what phrase(s) to Google to find my answer.

Answer 1

Convert your map to a SpatialPolygons object and then use gIntersection from the rgeos package. This clips your blue polygons exactly to the state boundaries. Inspired by this answer.

library(rgeos)
library(maps)
library(maptools)

mmap <- map('state', regions=c('maryland', 'virginia', 'delaware'), fill=TRUE)
IDs <- sapply(strsplit(mmap$names, ":"), function(x) x[1])
mmap.sp <- map2SpatialPolygons(mmap, IDs=IDs, proj4string=CRS("+proj=longlat +datum=WGS84"))

# make up some polygons to clip
mpoly <- readWKT("MULTIPOLYGON(((-80 36, -80 38, -75 38, -75 36, -80 36)), ((-77 39, -77 41, -74 41, -74 39, -77 39)))")
proj4string(mpoly) <- CRS("+proj=longlat +datum=WGS84")

plot(mmap.sp)
plot(mpoly, add=T, col="blue")

# clip
mpoly.clipped <- gIntersection(mmap.sp, mpoly, byid=TRUE)

plot(mmap.sp)
plot(mpoly.clipped, add=T, col="blue")