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

The above code (run time of several hours) results in this map (plotted with Panoply). Obviously there is something wrong here (compare to the map on http://protectedplanet.net).

Protected Planet (1 km^2)

Excluding all polygons with area below 100000 km^2 gives the following map.

a <- a[a$GIS_AREA >= 100000,]

Protected Planet (100000 km^2)

The second map is closer to the one on http://protectedplanet.net. From my understanding, all raster cells that are colored in the second map should also appear in the first map. Is that correct? And if so, what is wrong with my code?

flohump
  • 61
  • 1
  • 4
  • I do not think there is anything wrong with your code. I think the problem is that the WDPA data set has a number of bad (self-intersecting) polygons. I will have a look at rasterize and see if I can fix this. – Robert Hijmans Aug 31 '15 at 08:09
  • That would be great if you could have a look at the rasterize function. Looking forward to your reply. – flohump Sep 02 '15 at 08:13
  • I don't think the exclusions are helping. You might be excluding some rings/holes that are necessary for the geometry to be valid. – Benjamin Feb 09 '17 at 01:49
  • I've never needed more than 15GB of RAM to run a test case... Maximum not working example (MNWE)... – Benjamin Feb 09 '17 at 04:22
  • R fonction rasterize is kind of slow... I use now `v.to.rast.value` for SAGA to do it. It's annoying because you have to work outside R (However, RSAGA exists) but it's way faster so it compensate. – Bastien Feb 09 '17 at 15:57

0 Answers0