1

This question already has an answer here:

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?

marked as duplicate by whyzar, nmtoken, mgri, Richard Law, Mapperz Jun 6 '17 at 13:42

This question has been asked before and already has an answer. If those answers do not fully address your question, please ask a new question.

  • Suspect it's an issue with polygon geometry not being 100% correct. I ran the vector layer through v.clean using default values (tool = break / threshold = 0.100) and used that new layer with the rasterize command - wasn't 100% perfect, still had 18 'missing' islands instead of 77, but probably if you play with threshold and v.clean tools you'll fix it sufficiently. – Simbamangu Feb 19 '16 at 7:22
  • 1
    I think it is a bug in rasterize. I will look at that asap... – Robert Hijmans Feb 19 '16 at 7:56
  • This(gis.stackexchange.com/questions/200556/…) question has an answer that should work for you. I had the same problem as you and the person who asked the question. Someone supplied an explanation with a good solution. – shea Jun 5 '17 at 19:23