Take the 2-minute tour ×
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

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)
elevation<-raster("alt.bil")
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)
map.axes()

#add site localities
sites<-read.csv("sites.csv", header=T)
lat<-sites$Latitude
lon<-sites$Longitude

map(database="state", regions = nestates, interior=T, lwd=2)
points (x=lon, y=lat, pch=17, cex=1.5, col="black")
map.axes()
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
    polygon(c(0,360,360,0),c(0,0,90,90),col="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

share|improve this question
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
add comment

2 Answers

up vote 8 down vote accepted

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.

library(raster)
# 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)
r[]=runif(100*100)

# plot it with the boundaries we want to clip against:
plot(r)
plot(ne,add=TRUE)

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

# plot, and overlay:
plot(rr);plot(ne,add=TRUE)

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

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

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.

library(rgdal)
library(raster)

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(elevation.sub)
plot(state.sub, add = TRUE)

DEM of northeast US states

Cheers,
Florian

share|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. –  flowla May 23 '13 at 9:47
    
Thanks @RobertH for the improvement using mask :-) –  flowla May 27 '13 at 15:56
add comment

Your Answer

 
discard

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.