0

我有一个栅格网格，我想根据包“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
4

1 回答 1

0

裁剪不会像您想象的那样裁剪到国家多边形。事实上，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))

我相信这就是你所追求的。 裁剪网格

于 2014-05-23T13:00:56.677 回答