11

On extracting values of a raster to points I find that I have several NA's, and rather than use a buffer and fun arguments of extract function, instead I'd like to extract the nearest non-NA Pixel to a point that overlaps NA.

I am using the basic extract function:

data.extr<-extract(loc.thr, data[,11:10])

5 Answers 5

9

Here's a solution without using the buffer. However, it calculates a distance map separately for each point in your dataset, so it might be ineffective if your dataset is large.

set.seed(2)

# create a 10x10 raster
r <- raster(ncol=10,nrow=10, xmn=0, xmx=10, ymn=0,ymx=10)
r[] <- 1:10
r[sample(1:ncell(r), size = 25)] <- NA

# plot the raster
plot(r, axes=F, box=F)
segments(x0 = 0, y0 = 0:10, x1 = 10, y1 = 0:10, lty=2)
segments(y0 = 0, x0 = 0:10, y1 = 10, x1 = 0:10, lty=2)

# create sample points and add them to the plot
xy = data.frame(x=runif(10,1,10), y=runif(10,1,10))
points(xy, pch=3)
text(x = xy$x, y = xy$y, labels = as.character(1:nrow(xy)), pos=4, cex=0.7, xpd=NA)

# use normal extract function to show that NAs are extracted for some points
extracted = extract(x = r, y = xy)

# then take the raster value with lowest distance to point AND non-NA value in the raster
sampled = apply(X = xy, MARGIN = 1, FUN = function(xy) r@data@values[which.min(replace(distanceFromPoints(r, xy), is.na(r), NA))])

# show output of both procedures
print(data.frame(xy, extracted, sampled))

#          x        y extracted sampled
#1  5.398959 6.644767         6       6
#2  2.343222 8.599861        NA       3
#3  4.213563 3.563835         5       5
#4  9.663796 7.005031        10      10
#5  2.191348 2.354228        NA       2
#6  1.093731 9.835551         2       2
#7  2.481780 3.673097         3       3
#8  8.291729 2.035757         9       9
#9  8.819749 2.468808         9       9
#10 5.628536 9.496376         6       6
3
  • Your solution has successfully worked through the task with my raster which was not very large.
    – Joe
    Commented Dec 19, 2014 at 14:09
  • Hrm - but if you're not in memory.... then this no longer works. Thoughts?
    – jebyrnes
    Commented Sep 16, 2016 at 20:56
  • Nice solution but it makes R freeze for large number of points, 30k. Commented Jan 8, 2018 at 18:37
4

This is a raster-based solution, by first filling the NA pixels with the nearest non-NA pixel value. Note however, that this does not take into account the position of a point within a pixel. Instead, it calculates the distances between pixel centers to determine the nearest non-NA pixel.

First, it calculates for each NA raster pixel the distance and direction to the nearest non-NA pixel. The next step is to calculate the coordinates of this non-NA cell (assumes projected CRS), extract its value and to store this value at the NA location.

Starting data: a projected raster, with identical values as in the answer from koekenbakker:

set.seed(2)
# set projected CRS
r <- raster(ncol=10,nrow=10, xmn=0, xmx=10, ymn=0,ymx=10, crs='+proj=utm +zone=1') 
r[] <- 1:10
r[sample(1:ncell(r), size = 25)] <- NA

# create sample points
xy = data.frame(x=runif(10,1,10), y=runif(10,1,10))

# use normal extract function to show that NAs are extracted for some points
extracted <- raster::extract(x = r, y = xy)

Calculate the distance and direction from all NA pixels to the nearest non-NA pixel:

dist <- distance(r)  
# you can also set a maximum distance: dist[dist > maxdist] <- NA
direct <- direction(r, from=FALSE)

Retrieve coordinates of NA pixels

# NA raster
rna <- is.na(r) # returns NA raster

# store coordinates in new raster: https://stackoverflow.com/a/35592230/3752258 
na.x <- init(rna, 'x')
na.y <- init(rna, 'y')

# calculate coordinates of the nearest Non-NA pixel
# assume that we have a orthogonal, projected CRS, so we can use (Pythagorean) calculations
co.x <- na.x + dist * sin(direct)
co.y <- na.y + dist * cos(direct)

# matrix with point coordinates of nearest non-NA pixel
co <- cbind(co.x[], co.y[]) 

Extract values of nearest non-NA cell with coordinates 'co'

# extract values of nearest non-NA cell with coordinates co
NAVals <- raster::extract(r, co, method='simple') 
r.NAVals <- rna # initiate new raster
r.NAVals[] <- NAVals # store values in raster

Fill the original raster with the new values

# cover nearest non-NA value at NA locations of original raster
r.filled <- cover(x=r, y= r.NAVals)

sampled <- raster::extract(x = r.filled, y = xy)

# compare old and new values
print(data.frame(xy, extracted, sampled))

#           x        y extracted sampled
# 1  5.398959 6.644767         6       6
# 2  2.343222 8.599861        NA       3
# 3  4.213563 3.563835         5       5
# 4  9.663796 7.005031        10      10  
# 5  2.191348 2.354228        NA       3
# 6  1.093731 9.835551         2       2
# 7  2.481780 3.673097         3       3
# 8  8.291729 2.035757         9       9
# 9  8.819749 2.468808         9       9 
# 10 5.628536 9.496376         6       6

Note that point 5 takes another value than the answer of Koekenbakker, since this method does not take into account the position of the point within a pixel (as mentioned above). If this is important, this solution might not be appropriate. In other cases, e.g. if the raster cells are small compared to the point accuracy, this raster-based method should give good results.

2
  • Great explanation, i am glad you clarified this
    – Joe
    Commented May 26, 2016 at 6:17
  • perhaps this may qualify as a different question..but would be nice to have the new xy point with non-NA returned alongside the corresponding data
    – Joe
    Commented Nov 1, 2016 at 1:13
2

For a raster stack, use @koekenbakker's solution above, and turn it into a function. A raster stack's @layers slot is a list of rasters, so, lapply it across and go from there.

#new layer
r2 <- raster(ncol=10,nrow=10, xmn=0, xmx=10, ymn=0,ymx=10)
r2[] <- 1:10
r2[sample(1:ncell(r2), size = 25)] <- NA

#make the stack
r_stack <- stack(r, r2)

#a function for sampling
sample_raster_NA <- function(r, xy){
  apply(X = xy, MARGIN = 1, 
        FUN = function(xy) r@data@values[which.min(replace(distanceFromPoints(r, xy), is.na(r), NA))])

}

#lapply to get answers
lapply(r_stack@layers, function(a_layer) sample_raster_NA(a_layer, xy))

Or to be fancy (speed improvements?)

purrr::map(r_stack@layers, sample_raster_NA, xy=xy)

Which makes me wonder if the whole thing can be sped up even more using dplyr...

0

The function nearestLand in package seegSDM can do it quite well:

set.seed(2)

# create a 10x10 raster
r <- raster(ncol=10,nrow=10, xmn=0, xmx=10, ymn=0,ymx=10)
r[] <- 1:10
r[sample(1:ncell(r), size = 25)] <- NA
extracted = extract(x = r, y = xy)

# find points that fall in NA areas on the raster
plot(r)
segments(x0 = 0, y0 = 0:10, x1 = 10, y1 = 0:10, lty=2)
segments(y0 = 0, x0 = 0:10, y1 = 10, x1 = 0:10, lty=2)
outside_mask <- is.na(extracted)
outside_pts <- xy[outside_mask,]

points(xy, pch = 16)
text(x = xy$x, y = xy$y, labels = as.character(1:nrow(xy)), pos=4, cex=0.7, xpd=NA)

points(outside_pts, pch = 16, col = 'red')


# find the nearest cell
nearest.cells <- nearestLand(outside_pts, r, 90000)

# plot the new points (for those which were reassigned) in green
points(nearest.cells, pch = 16, col = ' green')

# plot where they moved to
arrows(outside_pts[, 1],
       outside_pts[, 2],
       nearest.cells[, 1],
       nearest.cells[, 2], length = 0.1)

#Extracting the raster values using the new coordinates:
extract(x = r, y = nearest.cells)
0

SEEG-Oxford/seegSDM

This github repository has a very useful function copied below. It returns the coordinate of the closest non-na pixel.

nearestLand <- function (points, raster, max_distance) {
  nearest <- function (lis, raster) {
    neighbours <- matrix(lis[[1]], ncol = 2)
    point <- lis[[2]]
    land <- !is.na(neighbours[, 2])
    if (!any(land)) {
      return (c(NA, NA))
    } else{
      coords <- xyFromCell(raster, neighbours[land, 1])   
      if (nrow(coords) == 1) {
        return (coords[1, ])
      }
      dists <- sqrt((coords[, 1] - point[1]) ^ 2 +
                      (coords[, 2] - point[2]) ^ 2)
      return (coords[which.min(dists), ])
    }
  }
  neighbour_list <- extract(raster, points,
                            buffer = max_distance,
                            cellnumbers = TRUE)
  neighbour_list <- lapply(1:nrow(points),
                           function(i) {
                             list(neighbours = neighbour_list[[i]],
                                  point = as.numeric(points[i, ]))
                           })
  return (t(sapply(neighbour_list, nearest, raster)))
}

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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