4

I have a raster representing a large landscape, where pixel values indicate particular land classes. I have been using the "distance" function in Package 'raster' in R to calculate the distance from all NA pixels to the nearest pixel with a value. In addition to the distance, I'm interested in knowing what the value was associated with the distance calculation. For example, if I have two land class values, e.g. 1= "Urban" and 2= "Agriculture", I'd like to know if a pixel is closest to an "Urban" pixel or an "Agriculture" pixel. I looked into the source code for the distance function, which converts the pixels to points and then calculates all relative distances but only returns the shortest distance. See below for an example of the distance function:

    library(raster)
    r <- raster(ncol=36,nrow=18)
    r[] <- NA
    r[500] <- 1
    r[100] <-2
    dist <- distance(r) #calculates the distance to the nearest non-NA pixel

Is there a function in a different package that does this (I'd like to stay in the R environment)?

  • 1
    In many cases there will be multiple nearest pixels. Wouldn't you want to know about this, and to know how many of those pixels are in each class? – whuber Mar 20 '16 at 23:09
  • @whuber...Insightful comment. It would be helpful though in my particular case I would likely favor one land cover type over another. – KevinB Mar 30 '16 at 23:06
5

You can do something along these lines:

Example data:

library(raster)
r <- raster(ncol=10,nrow=10, ext=extent(c(0,10,0,10)), crs='+proj=utm +zone=1')
r[] <- 1  
r[15:42] <- 2
set.seed(0)
r[sample(ncell(r), 50)] <- NA

Create separate layers for each class:

class1 <- reclassify(r, cbind(2, NA))
class2 <- reclassify(r, cbind(1, NA))

Compute the distance to each class:

dist1 <- distance(class1) 
dist2 <- distance(class2) 
s <- stack(dist1, dist2)
names(s) <- c('c1', 'c2')

Now you can extract by xy coordinates (or alternatively, cell numbers):

extract(s, cbind(3,4))

#     c1       c2
#[1,]  0 2.828427

For this location, the distance to 'c1' is zero (that is, the cell is class1) and the distance to the nearest cell of class2 is 2.8

If you want a RasterLayer, you can do:

x <- which.min(s) 

Note that there will be ties, and that the first layers is taken. This affects the following cells:

equals <- dist1 == dist2

That aside, to only see the nearest value for the NA cells:

y <- mask(x, r, inverse=TRUE)
  • 1
    It would help many readers if you would explain how this code is intended to work. – whuber Mar 20 '16 at 23:10
  • done as requested. – Robert Hijmans Mar 21 '16 at 10:19
  • 1
    Thank you. Your explanation shows a good start, but it doesn't quite seem to achieve an answer to the question. I interpret the question as requesting a raster (or equivalent data structure) whose values indicate the class of a nearest non-NA cell. – whuber Mar 21 '16 at 14:19
  • @RobertH... Thank you. This worked well and was much simpler than converting to points as I had originally thought. – KevinB Mar 30 '16 at 23:08
  • Hrm - but what if your raster has continuous values? – jebyrnes Sep 16 '16 at 20:19

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.