rr-rasterr-spgrass

Fill raster holes in R or Grass GIS


Sample data

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)

enter image description here

The problem

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 !


Solution

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