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 !

I am doing tasks like this usually with the same command as you: `x[c(15, 16, 17, 22)] <- Value` I am however also interested in seeing faster solutions. — maRtin, Dec 17 '16 at 11:37
Of course, that is how the holes were created in the first place; just for demo purposes. The trick is to deduce, algorithmically (not manually), that cells 15-17,22 are holes while the other empty cells are not ;-). — shekeine, Dec 17 '16 at 11:46

maRtin · Answer 1 · 2016-12-28T10:35:19.810

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)

Bastien · Accepted Answer · 2016-12-23T01:17:29.253

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)

edited Dec 23 '16 at 01:17

answered Dec 22 '16 at 19:57

Bastien

3,007
20
38

Many thanks for the attempt. However, cells 25, 26, 30, 31 (included in the polygon named `hole ` are not holes: see this for explanation [http://gis.stackexchange.com/questions/27255/how-to-identify-feature-vertices-that-are-part-of-a-donut-hole-in-arcgis-10. cheers – shekeine Dec 22 '16 at 22:39
can you adjust your reproducible example to represent that problem. It may help finding a solution. – Bastien Dec 22 '16 at 23:41
The question asks precisely that and the example accordingly contains both holes and open NA/NULL areas. Read: "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)?". Many thanks for trying however. – shekeine Dec 22 '16 at 23:49
I think I understand now. I don't have the answer (yet) but if I understand well, you have to transform the polygon `hole[1,]` from multipart to singlepart. That would give you 4 polygons. Than, you intersects those 4 polygones with the boundarie of your raster (`as(extent(x), "SpatialPolygons")`) and keep only the ones that do not touch. Then you have your "holes". – Bastien Dec 23 '16 at 00:27
@SheKeine, Changed it so to get closer following your last comment. – Bastien Dec 23 '16 at 01:18
@bastian Well, this is fast enough with the help of Grass e.g., I use Grass' `r.to.vect` instead of R's `rasterToPolygons`. Many thanks! – shekeine Dec 28 '16 at 09:03

score 0 · Answer 3 · answered Dec 28 '19 at 19:15

In this case, a cell is in a hole if it is bounded by at least one non-missing value in each of the four cardinal directions. You can test this using the code below, which also replaces the missing values with 99.

x2 <- x                        ## Create a copy of the raster
notNA <- !is.na(as.matrix(x))  ## Matrix identifying cells that are not NA

nr <- nrow(x) ## Number of rows
nc <- ncol(x) ## Number of columns

for(i in 2:(nr-1)) {       ## Loop over rows, ignoring first and last
    for(j in 2:(nc-1)) {   ## Loop over columns, ignoring first and last
        if(notNA[i,j]) ## Skip cell if not NA
            next
        ## For NA cells, determine if non-NA cell exists in each cardinal direction
        any.north <- any(notNA[1:(i-1),j])
        any.south <- any(notNA[(i+1):nr,j])
        any.west <- any(notNA[i,1:(j-1)])
        any.east <- any(notNA[i,(j+1):nc])
        if(any.north & any.south & any.west & any.east)
            x2[i,j] <- 99 ## If cell is in hole, repalce with new value
    }
}

Here's the new raster:

New raster with holes filled by the value 99

Fill raster holes in R or Grass GIS

Sample data

The problem

3 Answers3