Crop a RasterLayer with a SpatialPolygonDataFrame

Question

I have a raster grid I want to crop according to land borders of the world map provided by the data of the package 'maptools'. By doing some reasearch, I found that I have to use the crop() function and then the mask() function, but I get an error message.

Here is my code :

# load the worldmap SpatialPolygonDataFrame
library(maptools)
data(wrld_simpl)
ll=CRS("+init=epsg:4326")
world<-spTransform(wrld_simpl, ll)
ext<-extent(-10.417,31.917,34.083,71.083)
# get only region covering europe 
world@bbox<-as.matrix(ext)


# create the origin in WGS84 CRS
ll = CRS("+init=epsg:4326") 
origin = SpatialPoints(cbind(7,40),ll) 
# create the raster grid 
grid = raster() 
# define extent and resolution
e = extent(5,18,40,50)
extent(grid)<-e
res(grid)=c(0.08,0.08) 
# fill the raster values
grid[] = runif(1:20250)

# crop the grid with worldmap to erase sea areas
test<-crop(grid,world)
rasterize(world,test)
Found 246 region(s) and 3768 polygon(s)
Error in if (x1a > rxmx) { : argument is of length zero

Answer 1

Crop does not crop to the country polygons as you would think. In fact, the crop(grid, world) does not really crop anything, as the extent of world is larger than the extent of grid. You can check that by plotting test. It's just the same as grid.

Instead, after grid[] = runif(1:ncell(grid)) you could do:

# we're going to make 3 plots
par(mfrow=c(1,3))

# crop the Worldmap to grid 
worldcrop<-crop(world, grid)
plot(worldcrop)

# rasterize output, give cells value of NAME(seas are NA)
worldcropr = rasterize(worldcrop,grid, field='NAME', fun='first')
plot(worldcropr)

# mask random grid by worldcropr
gridmask = mask(x=grid, mask=worldcropr)
plot(gridmask)

# number of land pixels:
sum(!is.na(gridmask@data@values))

# number of sea pixels:
sum(is.na(gridmask@data@values))

I believe that's what you're after.