Extracting a value from a raster for a specific point based on the closest cell value using r

Question

I am trying to assign values to some site data which falls just outside the area for which I have weather data. I am trying to extract based on the nearest cell value and if possible a cell value within 40km.

My raster (r) looks like this :

class(r)

class       : RasterBrick 
dimensions  : 201, 464, 93264, 23376  (nrow, ncol, ncell, nlayers)
resolution  : 0.25, 0.25  (x, y)
extent      : -40.5, 75.5, 25.25, 75.5  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=longlat +datum=WGS84 
data source : \\ueahome\eressci5\zuw13bqu\data\NTProfile\Desktop\EOBS European data\rr_0.25deg_reg_v10.0.nc 
names       : X1950.01.01, X1950.01.02, X1950.01.03, X1950.01.04, X1950.01.05, X1950.01.06, X1950.01.07, X1950.01.08, X1950.01.09, X1950.01.10, X1950.01.11, X1950.01.12, X1950.01.13, X1950.01.14, X1950.01.15, ... 
Date        : 1950-01-01, 2013-12-31 (min, max)
varname     : rr 

I am extracting based on Latitude and Longitude data using the following code

vals <- extract(r, matrix(c(issues[22,3], issues[22,2]), ncol = 2), buffer = 40000)

However unfortunately I am getting the following output:

*can't attach a picture as I don't have enought reputation

X1950.01.01 X1950.01.02 X1950.01.03 X1950.01.04 X1950.01.05 X1950.01.06 X1950.01.07 X1950.01.08 X1950.01.09 X1950.01.10
1   0   4.8 4.6 0   0.2 0   0   0   0   0
2   0   4.7 4.5 0   1   0   0   0   0   0
3   0   4.7 4.5 0   1.1 0   0   0   0   0
4   0   4.6 4.3 0   1.2 0   0   0   0   0
5   NA  NA  NA  NA  NA  NA  NA  NA  NA  NA
6   NA  NA  NA  NA  NA  NA  NA  NA  NA  NA
7   NA  NA  NA  NA  NA  NA  NA  NA  NA  NA
8   0   4.1 3.9 0   0.7 0   0   0   0   0
9   0   4   3.7 0   0.9 0   0   0   0   0
10  0   4.1 3.8 0   1   0   0   0   0   0
11  0   4.1 3.8 0   1.1 0   0   0   0   0
12  NA  NA  NA  NA  NA  NA  NA  NA  NA  NA
13  NA  NA  NA  NA  NA  NA  NA  NA  NA  NA
14  NA  NA  NA  NA  NA  NA  NA  NA  NA  NA
15  NA  NA  NA  NA  NA  NA  NA  NA  NA  NA
16  NA  NA  NA  NA  NA  NA  NA  NA  NA  NA

*nb I have checked for this site and none of these rows are actually the closest.

How do I chose the cell value that is closest to the point without reducing the buffer size until one cell value appears (I have too many such site to do this for each site)?

Thanks in advance

Answer 1

When finding the cell of raster r that has the minimum distance to a point falling outside r, we don't need to calculate the distance to interior cells of r.

Specifically:

1. for points with x < xmin(r) and y < ymin(r), the closest of r is the lower left cell;
2. for points with x > xmax(r) and y < ymin(r), the closest of r is the lower right cell;
3. for points with x < xmin(r) and y > ymax(r), the closest of r is the upper left cell;
4. for points with x > xmax(r) and y > ymax(r), the closest of r is the upper right cell;
5. for points with xmin(r) < x < xmax(r) and y < ymin(r), the closest cell of r will lie along the bottom edge;
6. for points with xmin(r) < x < xmax(r) and y > ymax(r), the closest cell of r will lie along the top edge;
7. for points with x < xmin(r) and ymin(r) < y < ymax(r), the closest cell of r will lie along the left edge;
8. for points with xmin(r) < x < xmax(r) and y > ymax(r), the closest cell of r will lie along the top edge; and
9. points where xmin(r) < x < xmax(r) and ymin(r) < y < ymax(r) are already within the raster.

Here's a function that works out within which of these 9 regions each point falls, identifies the raster columns or rows that they correspond to, and identifies the cell at that column or row along the relevant edge/corner.

The function accepts either a 2 x n matrix of coordinates (x then y), or a SpatialPoints object. Set spatial=TRUE if you want the moved points returned as SpatialPoints.

Note that points are assumed to be in a planar coordinate system.

move.points <- function(r, pts, spatial=FALSE) {
  require(raster)
  require(sp)

  if (is(pts, 'SpatialPoints')) pts <- coordinates(pts)
  if (is(!r, 'Raster')) r <- raster(r)

  loc <- colSums(sapply(pts[, 1], '>', bbox(r)[1, ])) * 3 + 
    colSums(sapply(pts[, 2], '>', bbox(r)[2, ]))

  L <- split(as.data.frame(pts), loc)

  new.pts <- lapply(names(L), function(x) {
    switch(x, 
           '0' = xyFromCell(r, ncell(r) - ncol(r) + 1)[rep(1, nrow(L[[x]])), ],
           '1' = xyFromCell(r, cellFromXY(r, cbind(xmin(r), L[[x]][, 2]))),
           '2' = xyFromCell(r, 1)[rep(1, nrow(L[[x]])), ],
           '3' = xyFromCell(r, cellFromXY(r, cbind(L[[x]][, 1], ymin(r)))),
           '4' = {
             xy <- as.matrix(L[[x]])
             dimnames(xy) <- list(NULL, c('x', 'y'))
             xy
           },
           '5' = xyFromCell(r, cellFromXY(r, cbind(L[[x]][, 1], ymax(r)))),
           '6' = xyFromCell(r, ncell(r))[rep(1, nrow(L[[x]])), ],
           '7' = xyFromCell(r, cellFromXY(r, cbind(xmax(r), L[[x]][, 2]))),
           '8' = xyFromCell(r, ncol(r))[rep(1, nrow(L[[x]])), ]
    )
  })

  new.pts <- unsplit(mapply(function(x, y) {
    row.names(x) <- row.names(y)
    as.data.frame(x)
  }, new.pts, L, SIMPLIFY=FALSE), loc)

  colnames(new.pts) <- colnames(pts)
  if(isTRUE(spatial)) new.pts <- SpatialPoints(new.pts)  
  return(new.pts)
}

And an example for a 1000 x 1000 cell raster with 1000 points randomly placed (~890 cells fall outside the raster's extent):

r <- raster(matrix(runif(1e6), ncol=1e3))
pts <- cbind(lon=runif(1e3, -1, 2), lat=runif(1e3, -1, 2))
system.time(new.pts <- move.points(r, pts, spatial=FALSE))

#   user  system elapsed 
#   0.03    0.00    0.04 

For a 2000 x 2000 cell raster, with 1 million points (~890000 points outside the raster's extent):

r <- raster(matrix(runif(4e6), ncol=2e3))
pts <- cbind(lon=runif(1e6, -1, 2), lat=runif(1e6, -1, 2))
system.time(new.pts <- move.points(r, pts, spatial=FALSE))

#   user  system elapsed 
#  12.71    0.23   13.12 

And here's a plot showing how the points from the first example above have been moved.

plot(pts, asp=1, pch=20, panel.first=abline(h=0:1, v=0:1, lwd=2, col='gray'))
segments(pts[, 1], pts[, 2], new.pts[, 1], new.pts[, 2], col='#00000050')
plot(extent(r), add=TRUE, lwd=3)

Once you've done this, you can use extract(r, new.pts).