Fill raster holes in R or Grass GIS

Question

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)

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 !

Answer 1

Here a solution without polygonization: (it's not elegant, but it works). However you have to reclassify your holes/island into values (i.e. 999) and all other non island to NA. Like this:

x <- raster(x=matrix(rep(NA,36), nrow=6), xmn=-1000, xmx=1000, ymn=-100, ymx=900)
x[c(8, 15, 16, 17, 22, 25, 26, 30, 31)] <- 999

plot(x)

Now I use the clump() function to check if there are island and the cool thing about that function is, that it also returns IDs of these island:

#Get Islands with IDs
cl <- clump(x,directions=8)
plot(cl)

I then create a dataframe out of the frequencies of the islands (this is just to get the ID of each island)

freqCl <- as.data.frame(freq(cl))

#remove the (row) which corresponds to the NA values (this is important for the last step)
freqCl <- freqCl[-which(is.na(freqCl$value)),]

Check if any island touches a border:

#Check if the island touches any border and therefore isn't a "real island" (first and last column or row)

noIslandID <- c()
#First row
if(any(rownames(freqCl) %in% cl[1,])){
  eliminate   <- rownames(freqCl)[rownames(freqCl) %in% cl[1,]]
  noIslandID  <- append(noIslandID, eliminate)
}
#Last row
if(any(rownames(freqCl) %in% cl[nrow(cl),])){
  eliminate  <- rownames(freqCl)[rownames(freqCl) %in% cl[nrow(cl),]]
  noIslandID <- append(noIslandID, eliminate)
}
#First col
if(any(rownames(freqCl) %in% cl[,1])){
  eliminate   <- rownames(freqCl)[rownames(freqCl) %in% cl[,1]]
  noIslandID  <- append(noIslandID, eliminate)
}
#Last col
if(any(rownames(freqCl) %in% cl[,ncol(cl)])){
  eliminate   <- rownames(freqCl)[rownames(freqCl) %in% cl[,ncol(cl)]]
  noIslandID  <- append(noIslandID, eliminate)
}

Eliminate those island which touch a border:

noIslandID <- unique(noIslandID)
IslandID   <- setdiff(rownames(freqCl), noIslandID)

In a last step, assign a 1 to every "real island" from the inital raster:

for(i in 1:length(IslandID)) {
  x[cl[]==as.numeric(IslandID[i])] <- 1
}

plot(x)

Answer 2

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)