x <- raster(x=matrix(data=1:36, nrow=6), xmn=-1000, xmx=1000, ymn=-100, ymx=900)
x[c(8, 15, 16, 17, 22, 25, 26, 30, 31)] <- NA
plot(x)
How do I distinguish (algorithmically) the holes in the raster i.e., the area bounded by cells c(15:17, 22) from the other gaps that are not holes (i.e., the rest of the empty cells)?
This would make it possible to do operations only on the hole / island regions of the raster, fill holes with a custom value etc etc.
The actual rasters have around 30000 holes and therefore speed is important. I am interested in both R and Grass GIS solutions. Many thanks for your help, much appreciated !
What about polygonizing? For speed I don't know what it will worth, but you could:
x[!is.na(values(x))]<-1
plot(x)
x[is.na(values(x))]<-0
hole <- rasterToPolygons(x, fun=NULL, n=4, na.rm=TRUE, digits=12, dissolve=T)
Now you have to make your singlepart polygons to multipart:
hole2 <- SpatialPolygons(lapply(hole@polygons[[1]]@Polygons, function(xx) Polygons(list(xx),ID=round(runif(1,1,100000000)))))
plot(hole2, add=T)
Now you find the "real" holes, which are the ones not touching the border
around <- as(extent(x), "SpatialLines")
touch_border <- gIntersects(around, hole2, byid=T)
extract(x, hole2[!touch_border,],cellnumbers=T)
That gives you the cells number by hole. It found the cell "8" which you're not saying it's a hole so I'm not sure if it's correct, but it must be very close!
If it's slow in R, do the same algorithm in GRASS (mostly the rasterToPolygons
call)