library(raster)
library(rgdal)
#download data
download.file("http://www.protectedplanet.net/downloads/WDPA_July2015?type=shapefile","WDPA_July2015-shapefile.zip")
#unzip files
unzip("WDPA_July2015-shapefile.zip")
#read in shapefile
a <- shapefile("WDPA_July2015-shapefile-polygons.shx",verbose=TRUE)
#Exclude polygons with areas smaller than 1 km^2
a <- a[a$GIS_AREA >= 1,]
#generate 0.5 degree raster
r <- raster(ncols=720, nrows=360)
#convert polygon cover to raster (as cell share)
x <- rasterize(a,r,getCover=TRUE,progress="text")
#write output as netCDF
writeRaster(x,"WDPA_July2015_raster.nc",xname="lon",yname="lat",overwrite=TRUE)
上記のコード (数時間の実行時間) は、このマップ (Panoply でプロット) になります。明らかに、ここには何か問題があります ( http://protectedplanet.netの地図と比較してください)。
面積が 100000 km^2 未満のすべてのポリゴンを除外すると、次のマップが得られます。
a <- a[a$GIS_AREA >= 100000,]
2 番目のマップはhttp://protectedplanet.netのマップに近いものです。私の理解では、2 番目のマップで色付けされたすべてのラスター セルは、最初のマップにも表示されるはずです。あれは正しいですか？もしそうなら、私のコードの何が問題なのですか?