By I Del Toro


2013-05-20 03:07:18 8 Comments

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

2 comments

@fdetsch 2013-05-22 15:02:23

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([email protected]$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

@I Del Toro 2013-05-23 00:40:51

Where did you get the state shapefile?

@fdetsch 2013-05-23 09:47:39

@IDelToro, I got it from Geocommons.

@ecologist1234 2017-10-27 17:21:57

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?

@fdetsch 2018-08-24 07:39:27

@ecologist1234, can you provide an example?

@Spacedman 2013-05-20 11:31: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.

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.

@Spacedman 2013-05-26 11:58:43

Cheers Robert, nice edit. I think I'd forgotten about mask.

Related Questions

Sponsored Content

1 Answered Questions

Clipping raster brick object in R

  • 2016-04-11 07:47:48
  • Bhoj Raj
  • 1759 View
  • 6 Score
  • 1 Answer
  • Tags:   raster r clip

1 Answered Questions

[SOLVED] QGIS Clipping Raster Layer

  • 2017-02-14 20:35:56
  • joe d
  • 2167 View
  • 3 Score
  • 1 Answer
  • Tags:   qgis raster clip

1 Answered Questions

[SOLVED] QGIS Error clipping raster?

2 Answered Questions

[SOLVED] R GIS layers not lining up (spplot/plot)

  • 2016-07-11 16:14:52
  • Adam
  • 580 View
  • 1 Score
  • 2 Answer
  • Tags:   r mapping plot

1 Answered Questions

Sponsored Content