28

I have two SpatialPolygonsDataFrame files: dat1, dat2

extent(dat1)
class       : Extent 
xmin        : -180 
xmax        : 180 
ymin        : -90 
ymax        : 90 


extent(dat2)
class       : Extent 
xmin        : -120.0014 
xmax        : -109.9997 
ymin        : 48.99944 
ymax        : 60

I want to crop the file dat1 using the extent of dat2. I don't know how to do it. I just handle raster files using "crop" function before.

When I use this function for my current data, the following error occurs:

> r1 <- crop(BiomassCarbon.shp,alberta.shp)
Error in function (classes, fdef, mtable)  : 

 unable to find an inherited method for function ‘crop’ for signature"SpatialPolygonsDataFrame"’
agstudy
  • 119,832
  • 17
  • 199
  • 261
Jian Zhang
  • 1,173
  • 5
  • 15
  • 20
  • this is hopeless, work on the details for a question involving dat1 and dat2, or those other things – mdsumner Dec 21 '12 at 06:03

4 Answers4

65

Streamlined method added 2014-10-9:

raster::crop() can be used to crop Spatial* (as well as Raster*) objects.

For example, here's how you might use it to crop a SpatialPolygons* object:

## Load raster package and an example SpatialPolygonsDataFrame
library(raster) 
data("wrld_simpl", package="maptools")

## Crop to the desired extent, then plot
out <- crop(wrld_simpl, extent(130, 180, 40, 70))
plot(out, col="khaki", bg="azure2")

Original (and still functional) answer:

The rgeos function gIntersection() makes this pretty straightforward.

Using mnel's nifty example as a jumping off point:

library(maptools)
library(raster)   ## To convert an "Extent" object to a "SpatialPolygons" object.
library(rgeos)
data(wrld_simpl)

## Create the clipping polygon
CP <- as(extent(130, 180, 40, 70), "SpatialPolygons")
proj4string(CP) <- CRS(proj4string(wrld_simpl))

## Clip the map
out <- gIntersection(wrld_simpl, CP, byid=TRUE)

## Plot the output
plot(out, col="khaki", bg="azure2")

enter image description here

Josh O'Brien
  • 159,210
  • 26
  • 366
  • 455
  • 1
    Thanks for providing a code example. I just did not have time when I posted my response. +1 for your clever solution in using the raster package to create the bounding polygon. – Jeffrey Evans Dec 21 '12 at 17:04
  • 1
    @JeffreyEvans -- Yeah. I find that **raster** does a *much* better job of providing user friendly conversion funtions/methods for users than does **sp**. I appreciate that the **sp** authors may want to keep their package spare and lightweight (since its intended to undergird so many other packages) but I still hope they eventually add a few more utility functions. – Josh O'Brien Dec 21 '12 at 19:04
  • 1
    Neat: `as(extent(130, 180, 40, 70), "SpatialPolygons")` Thank you! – geotheory Jul 09 '14 at 21:55
  • crop is great - but is there an inverse crop in raster that does the same thing, but with the inverse of the 2nd shapefile? – jebyrnes Jan 12 '18 at 17:29
  • @jebyrnes Have a look at the example in `?raster::erase` for one way to do that. – Josh O'Brien Jan 12 '18 at 17:50
7

Here is an example of how to do this with rgeos using the world map as an example

This comes from Roger Bivand on R-sig-Geo mailing list. Roger is one of the authors of the sp package.

Using the world map as an example

library(maptools)

data(wrld_simpl)

# interested in the arealy bounded by the following rectangle
# rect(130, 40, 180, 70)

library(rgeos)
# create  a polygon that defines the boundary
bnds <- cbind(x=c(130, 130, 180, 180, 130), y=c(40, 70, 70, 40, 40))
# convert to a spatial polygons object with the same CRS
SP <- SpatialPolygons(list(Polygons(list(Polygon(bnds)), "1")),
proj4string=CRS(proj4string(wrld_simpl)))
# find the intersection with the original SPDF
gI <- gIntersects(wrld_simpl, SP, byid=TRUE)
# create the new spatial polygons object.
out <- vector(mode="list", length=length(which(gI)))
ii <- 1
for (i in seq(along=gI)) if (gI[i]) {
  out[[ii]] <- gIntersection(wrld_simpl[i,], SP)
  row.names(out[[ii]]) <- row.names(wrld_simpl)[i]; ii <- ii+1
}
# use rbind.SpatialPolygons method to combine into a new object.
out1 <- do.call("rbind", out)
# look here is Eastern Russia and a bit of Japan and China.
plot(out1, col = "khaki", bg = "azure2")

enter image description here

mnel
  • 113,303
  • 27
  • 265
  • 254
1

You cannot use crop on sp polygon objects. You will need to create a polygon representing the bbox coordinates of dat2 and then can use gIntersects in the rgeos library.

Edit: This comment was in relation to the version available in 2012 and this is no longer the case.

Jeffrey Evans
  • 2,325
  • 12
  • 18
0

see ?crop

corp(x, y, filename="", snap='near', datatype=NULL, ...)

x Raster* object

y Extent object, or any object from which an Extent object can be extracted (see Details

You need to rasterize the first SpatialPolygon, using rasterize function from the raster package

I create some data to show how to use rasterize:

n <- 1000
x <- runif(n) * 360 - 180
y <- runif(n) * 180 - 90
xy <- cbind(x, y)
vals <- 1:n
p1 <- data.frame(xy, name=vals)
p2 <- data.frame(xy, name=vals)
coordinates(p1) <- ~x+y
coordinates(p2) <- ~x+y

if I try :

 crop(p1,p2)
 unable to find an inherited method for function ‘crop’ for signature ‘"SpatialPointsDataFrame"’

Now using rasterize

r <- rasterize(p1, r, 'name', fun=min)
crop(r,p2)

class       : RasterLayer 
dimensions  : 18, 36, 648  (nrow, ncol, ncell)
resolution  : 10, 10  (x, y)
extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 
data source : in memory
names       : layer 
values      : 1, 997  (min, max)
agstudy
  • 119,832
  • 17
  • 199
  • 261