我需要将 shapefile 转换为光栅格式。
我在 R 包“raster”中使用了函数“rasterize”,但结果看起来不正确。
tst <- rasterize(shpfile, r, fun="count")
Found 5 region(s) and 5 polygon(s)
没有出现记录的网格:
range(tst[],na.rm=TRUE)
[1] Inf -Inf
Warning messages:
1: In min(x) : no non-missing arguments to min; returning Inf
2: In max(x) : no non-missing arguments to max; returning -Inf
sum(tst[],na.rm=TRUE)
[1] 0
我写的 R 脚本:
# download the GIS shape file
download.file("http://esp.cr.usgs.gov/data/little/abiebrac.zip",
destfile = "abiebrac.zip")
# unzip
unzip("abiebrac.zip")
# Read shapefile into R
library(maptools)
shpfile <- readShapePoly("abiebrac")
extent(shpfile)
# mapping
plot(shpfile)
library(maps)
map("world",xlim=c(-180,-50),ylim=c(7,83),add=FALSE)
plot(shpfile,add=TRUE,lwd=10)
# rasterize
library(raster)
GridSize <- 0.5 # degree
r <- raster(ncol= round(abs(-180-(-50))/GridSize),
nrow=round(abs(83-7)/GridSize))
extent(r) <- extent(c(-180, -50, 7, 83))
tst <- rasterize(shpfile, r, fun="count")
# summary
sum(tst[],na.rm=TRUE)
range(tst[],na.rm=TRUE)
# mapping
plot(tst,col="red",main="abiebrac")
map("world",xlim=c(-180,-50),ylim=c(7,83),add=TRUE)