2013-02-21 82 views
2

我需要將shapefile轉換爲光柵格式。在R中將GIS形狀文件轉換爲柵格?

我在R包「光柵」中使用了「光柵化」功能,但結果看起來不正確。

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) 

回答

6

我不知道你爲什麼使用「計數」的樂趣論點,但在這種情況下,因爲沒有重疊,它正在產生NA結果。您還需要在spatialPolygonDataFrame對象中定義一個屬性字段以將值分配給柵格。您也可以直接從sp對象中拉出範圍。

此代碼似乎產生你想要的。

require(raster) 
require(rgdal) 
require(sp) 

setwd("D:/TMP") 
shpfile <- readOGR(getwd(), "abiebrac") 
r <- raster(extent(shpfile)) 
    res(r)=0.05 
    r <- rasterize(shpfile, field="ABIEBRAC_", r) 

plot(r) 
    plot(shpfile,lwd=10,add=TRUE) 
+1

我用https://github.com/brry/misc/blob/master/shp2raster.R自動化了這個 – 2016-11-25 11:02:23