Обрезка растрового слоя с помощью SpatialPolygonDataFrame

У меня есть растровая сетка, которую я хочу обрезать в соответствии с границами карты мира, предоставленной данными пакета «maptools». Проведя некоторые исследования, я обнаружил, что должен использовать функцию crop(), а затем функцию mask(), но получаю сообщение об ошибке.

Вот мой код:

# 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

person Remssssss    schedule 23.05.2014    source источник


Ответы (1)


Обрезка не обрезает полигоны страны, как вы думаете. На самом деле, crop(grid, world) ничего не обрезает, так как размер мира больше, чем размер сетки. Вы можете проверить это, построив график test. Это точно так же, как grid.

Вместо этого после grid[] = runif(1:ncell(grid)) вы можете сделать:

# 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))

Я считаю, что это то, что вы после. обрезанная сетка

person koekenbakker    schedule 23.05.2014