3

In R, how could I make the distinction between inner and outer NA's in a raster with some shape having NA's both around but also inside?

In the example below, how could I for exemple select only the NA's outside the R logo (i.e., how can I make everything which is included in the circle of the logo appear as white)?

library(raster)
r <- raster(system.file("external/rlogo.grd", package="raster"))
r[r>240] = NA
par(mfrow=c(1,2))
plot(r, main='r')
plot(is.na(r), main="is.na(r)")

enter image description here

ztl
  • 2,512
  • 1
  • 26
  • 40

2 Answers2

2

You don't really have many options. This type of analysis usually requires some more elaborate methods. Here however is a simple workaround using the clumpfunction:

#your inital code
library(raster)
r <- raster(system.file("external/rlogo.grd", package="raster"))
rna <- rc <- r
rna[r>240] = NA
par(mfrow=c(2,2))

#reclass values <=240 to NA (needed for clump function. 
#Here, NAs are used to seperate clumps)
rc[r<=240] <- NA
rc <- clump(rc)

#what you get after applying the clump function
#are homogenous areas that are separated by NAs. 
#Let's reclassify all areas with an ID > 8. 
#In this case, these are the areas inside the ring.
rc_reclass <- rc
rc_reclass[rc_reclass>8]  <- 100

#let's do some plotting
plot(r, main='r')
plot(is.na(rna), main="is.na(r)")
plot(rc, main="clumps")
plot(rc_reclass, main="clumps reclass")

enter image description here

maRtin
  • 6,336
  • 11
  • 43
  • 66
  • Thanks - I also think `clump` is the good approach. Technically, the desired output was here `rc_reclass = clump(is.na(r))` and then `rc_reclass[rc_reclass>1]=NA` to only select the outer `NA`'s, but thus you provide an adequate solution. Thanks! – ztl Jun 07 '17 at 15:44
2

I agree with @maRtin, it's a bit tricky. Not only don't you have a dedicated NoData value, but also the image is a bit smudgy.

Nevertheless, I think I found a way which is a bit better than clump, which uses the spatial domain to separate the areas:

First, I'm getting the focal values of the pixel neighbourhoods:

#make copy
r2 <- r

# focal values
fv <- getValuesFocal(r2,ngb = c(3,3))

Then I I first I exclude all pixels, which have a neighbourhood value of greater than 242.8. This was purely trial and error, but it gives a good result.

ix <- rowMeans(fv,na.rm = T) > 242.8
r2[ix] <- NA

You could actually already deem this acceptable. The only problem is, that there is a small border around the value area which should be NA.

enter image description here

So somehow I need to get rid of the remaining NA pixels. This I try to do with an iterative exclusion. For every iteration, I see if there's pixels which have still NA values around and a max value lower than a certain threshold. Again, there's a lot of playing around involved and I guess you could achieve a better result than this, but I guess this would be a way to go.

while (TRUE){

  fv <- getValuesFocal(r2,ngb = c(3,3))
  ix <- apply(fv,1,function(x) max(x,na.rm=T)) > 243 & rowSums(is.na(fv)) > 0

  if (any(ix)){

    r2[ix] <- NA


  } else {
    break
  }
}

After a couple of iterations, I get this:

enter image description here

There are clearly already some pixels gone which shouldn't be, maybe it could be done with a bit more fiddling around.

Another interesting thought would be looking at all three channels. If you load the image with brick, you can get the RGB channels. I've tried a few things like max, mean, mode, sd, etc. but to no avail.

Val
  • 6,585
  • 5
  • 22
  • 52