I am having problem in rasterizing a polygon shapefile: "island" polygons (I mean polygons which are totally included into a larger one) are ignored and do not appear in the raster. My code is:
seg <- readOGR(dsn="my path", layer="seg_10_9_5s") map <- raster(map2007.tif) # import of raster of same extension than shape # the shape has only a field with a unique id for each polygon length(unique(seg@data$id)) # give me the right number of polygons (1568) print(seg@data) # everything seem OK nrow(seg@data) # 1568 idpol <- unique(seg@data$id) # give me a vector of 1568 id values # rasterize the polygons segr <- rasterize(seg,map,"id") writeRaster(segr, filename="seg.tif", mat="GTiff", overwrite=TRUE,datatype='INT2S') # I look for the values of grid cells freq <- freq(segr) # freq has only 1492 rows! f <- freq[,1] f # I compare the list of value in vector and raster setdiff(idpol,f) # id of polygons which are absent from the raster
In QGIS I display the shape and seg.tif, the polygons identified as lacking in raster were not rasterized, the cells get the value of the surrounding polygon.
Maybe the problem is with the original shapefile? but it seems OK.
I put the shape and the map in https://dl.dropboxusercontent.com/u/22467480/misterio.zip
Can anyone see what I'm missing here?