I am building a map for the northeastern U.S.A. The map background needs to be either an altitude map or a mean annual temperature map. I have two rasters from Worldclim.org which give me these variables but I need to clip them to the extent of the states I am interested in. Any suggestions on how to do this. This is what I have so far:

#load libraries
library (sp)
library (rgdal)
library (raster)
library (maps)
library (mapproj)

#load data
state<- data (stateMapEnv)
meantemp<-raster ("bio_1.asc")

#build the raw map
nestates<- c("maine", "vermont", "massachusetts", "new hampshire" ,"connecticut",
  "rhode island","new york","pennsylvania", "new jersey",
  "maryland", "delaware", "virginia", "west virginia")

map(database="state", regions = nestates, interior=T,  lwd=2)

#add site localities
sites<-read.csv("sites.csv", header=T)

map(database="state", regions = nestates, interior=T, lwd=2)
points (x=lon, y=lat, pch=17, cex=1.5, col="black")
library(maps)                                                                  #Add axes to  main map
map.scale(x=-73,y=38, relwidth=0.15, metric=T,  ratio=F)

#create an inset map

 # Next, we create a new graphics space in the lower-right hand corner.  The numbers are proportional distances within the graphics window (xmin,xmax,ymin,ymax) on a scale of 0 to 1.
  # "plt" is the key parameter to adjust
    par(plt = c(0.1, 0.53, 0.57, 0.90), new = TRUE)

  # I think this is the key command from http://www.stat.auckland.ac.nz/~paul/RGraphics/examples-map.R
    plot.window(xlim=c(-127, -66),ylim=c(23,53))

  # fill the box with white

  # draw the map
    map(database="state", interior=T, add=TRUE, fill=FALSE)
    map(database="state", regions=nestates, interior=TRUE, add=TRUE, fill=TRUE, col="grey")

The elevation and meantemp objects are the ones that need to be clipped to the area extent of the nestates object. Any input would help

  • 2
    Any chance you could make this reproducible by others, by perhaps creating rasters from random data with the same extent and resolution? – Spacedman May 20 '13 at 11:06
  • why not skip the crop piece and go to mask straight? – Diego Mar 13 at 20:12

I would drop using the maps package and find a state shapefile. Then load that into R using rgdal, and then do some polygon overlay work.

# use state bounds from gadm website:
# us = shapefile("USA_adm1.shp")
us <- getData("GADM", country="USA", level=1)
# extract states (need to uppercase everything)
nestates <- c("Maine", "Vermont", "Massachusetts", "New Hampshire" ,"Connecticut",
         "Rhode Island","New York","Pennsylvania", "New Jersey",
         "Maryland", "Delaware", "Virginia", "West Virginia")

ne = us[match(toupper(nestates),toupper(us$NAME_1)),]

# create a random raster over the space:        
r = raster(xmn=-85,xmx=-65,ymn=36,ymx=48,nrow=100,ncol=100)

# plot it with the boundaries we want to clip against:

# now use the mask function
rr <- mask(r, ne)

# plot, and overlay:

How's that? The gadm shapefile is quite detailed, you might want to find a more generalised one instead.

|improve this answer|||||
  • Cheers Robert, nice edit. I think I'd forgotten about mask. – Spacedman May 26 '13 at 11:58

Here is an approach using extract() from the raster package. I tested it with altitude and mean temperature data from the WorldClim website (I limit this example to altitude, temperature works similar), and an appropriate shapefile of the US containing state borders is to be found here. Just download the .zip data and decompress it to your working directory.

You need to load rgdal and raster libraries in order to proceed.


Let's import the US shapefile now using readOGR(). After setting the shapefile's CRS, I create a subset containing the desired states. Pay attention to the use of capital and small initial letters!

state <- readOGR(dsn = path.data, layer = "usa_state_shapefile")
projection(state) <- CRS("+proj=longlat +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +no_defs")

# Subset US shapefile by desired states
nestates <- c("Maine", "Vermont", "Massachusetts", "New Hampshire" ,"Connecticut",
             "Rhode Island","New York","Pennsylvania", "New Jersey",
             "Maryland", "Delaware", "Virginia", "West Virginia")

state.sub <- state[as.character(state@data$STATE_NAME) %in% nestates, ]

Next, import the raster data using raster() and crop it with the extent of the previously generated states subset.

elevation <- raster("/path/to/data/alt.bil")

# Crop elevation data by extent of state subset
elevation.sub <- crop(elevation, extent(state.sub))

As a final step, you need to identify those pixels of your elevation raster that lie within the borders of the given state polygons. Use the 'mask' function for that.

elevation.sub <- mask(elevation.sub, state.sub)

Here's a very simple plot of the results:

plot(state.sub, add = TRUE)

DEM of northeast US states


|improve this answer|||||
  • Where did you get the state shapefile? – I Del Toro May 23 '13 at 0:40
  • @IDelToro, I got it from Geocommons. – fdetsch May 23 '13 at 9:47
  • Why does this take so long (>>15 min, maybe hours) when working with a ~11mb rasterlayer and a single-polygon shape file? Is there a more efficient method? – ecologist1234 Oct 27 '17 at 17:21
  • @ecologist1234, can you provide an example? – fdetsch Aug 24 '18 at 7:39

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.