2

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):

image1

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.

2

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")

enter image description here

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

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

enter image description here

  • Thank you so much! You just saved me hours of work! I'm very new to maps in R, so this was really starting to frustrate me. Again, thanks a million! – lassa8 Jun 13 '15 at 19:17

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.