1

I'm attempting to merge some states from a shapefile, and produce a raster that I can use downstream. I have gotten the states merged, however when I am creating an empty raster to rasterize with the crop function seems to fail. I'm pretty new to the GIS features in R and really appreciate the help. Shapefile is from http://www.arcgis.com/home/item.html?id=f7f805eb65eb4ab787a0a3e1116ca7e5

library(maptools)
library(shapefiles)
library(raster)
usa.states <- readOGR(dsn = "states_21basic/", layer = "states")
head(usa.states)
Co=usa.states[usa.states@data$STATE_NAME== "Colorado",]
Nm=usa.states[usa.states@data$STATE_NAME== "New Mexico",]
Az=usa.states[usa.states@data$STATE_NAME== "Arizona",]
Ut=usa.states[usa.states@data$STATE_NAME== "Utah",]
Corners= spRbind(spRbind(spRbind(Co,Ut),Nm),Az)
CRS="+proj=longlat +datum=WGS84"
Corners=spTransform(Corners, CRS(CRS))
> extent(Corners)
class       : Extent 
xmin        : -114.8218 
xmax        : -102.0372 
ymin        : 31.33563 
ymax        : 42.0023 
cor.ext=extent(Corners)
r<-raster(ncol=ncol(Corners), nrow=nrow(Corners), crs=CRS)
Corners.crop= crop(r,cor.ext, snap="out")

When I then call the extent of the 'Corners.crop' however I receive:

> extent(Corners.crop)
class       : Extent 
xmin        : -180 
xmax        : -36 
ymin        : 0 
ymax        : 45

I'm confused to what I'm missing to get this to work. I am also looking to have a 1Km resolution and am curious if it would be better to change the resolution on the empty raster or after I rasterize shape.

Reed
  • 308
  • 2
  • 14
  • `readOGR` is in the `rgdal` package, `crop` in the `raster` package. No need to reproject the data, they already are in `+proj=longlat +datum=NAD83 +no_defs +ellps=GRS80 +towgs84=0,0,0`. And did you check the output of `ncol(Corners)` and `nrow(Corners)`? –  Feb 26 '15 at 06:18
  • @Pascal did leave off some of the packages I had installed, I'm sorry. I know I didn't really need to reproject the data, but it seemed like an easy place to do it for the final product. As far as the 'ncol' and 'nrow' go that segment will be replaced with something that gets me better resolution. – Reed Feb 26 '15 at 14:45

1 Answers1

4
library(rgdal)
library(raster)
library(rgeos)

usa.states <- readOGR("states.shp", layer = "states")

# Here we subset once
Corners <- usa.states[usa.states$STATE_NAME %in% c("Colorado", "New Mexico","Arizona","Utah"),]

enter image description here

# Dissolve polygons into one
Corners <- gUnaryUnion(Corners)

enter image description here

# Create a 20x20 raster using the extent of Corners
# The number of rows and columns can be change to increase/reduce the resolution
r <- raster(extent(Corners), ncol=20, nrow=20, crs=CRS(proj4string(Corners)))

# Rasterize
Corners.crop <- rasterize(Corners, r)

enter image description here

  • Great answer for the most part, however I really need the 1Km finial resolution. Is there some way to accurately set that (that take into effect the unequal spacing of lines of Latitude) or is it just easier to approximate and change later? – Reed Feb 26 '15 at 14:34
  • 1
    As written, "The number of rows and columns can be change to increase/reduce the resolution". With this projection, it will be around 1km. I let you calculate the number of rows and columns, according to the range in latitude and longitude. –  Feb 26 '15 at 14:38