Search code examples
rrastershapefiler-maptools

Raster not cropping to shapefile


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.


Solution

  • 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