Wednesday, 21 February 2018

Clipping a raster in R


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



Answer



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.



No comments:

Post a Comment

arcpy - Changing output name when exporting data driven pages to JPG?

Is there a way to save the output JPG, changing the output file name to the page name, instead of page number? I mean changing the script fo...