[GIS] Extract median value for polygons from multiresolution raster data, using R (velox, raster -package)

extractrraster

I am trying to extract median values from rasters for 1043 different polygons (plots).

Raster resolution varies from 1 to 30 meters and the plot size is r = 6 meters. As you can see, at coarser resolution data the polygon doesn't always contain any pixel center.

The polygons area in a SpatialPolygonsDataFrame and I am trying to use the veloxpackage for the extraction. For the smaller resolution it works beautifully;

v1$extract(spdf,fun = median) # v1 = veloxraster, 1 m resolution

returns the median value for the plots. However;

v30$extract(spdf,fun = median) # v30 = veloxraster, 30 m resolution

returns the median for those plots which are divided to multiple pixels (or pixel centers), and returns NA for those polygons, which do not coincide with any pixel center coordinates.

The raster package has an option small = TRUE, where the single pixel value is returned. According to the packages' documentation:

Logical. If TRUE and y represents points and a buffer argument is
used, the function always return a number, also when the buffer does
not include the center of a single cell. The value of the cell in
which the point falls is returned if no cell center is within the
buffer. If y represents polygons, a value is also returned for
relatively small polygons (e.g. those smaller than a single cell of
the Raster* object), or polygons with an odd shape, for which
otherwise no values are returned because they do not cover any raster
cell centers.

What would be an easy way to achieve that with velox? The extraction speed difference is huge (1.5 s / raster vs. 360-380 s / raster), so I'd like to perform these operation with velox if by any means possible. 

# Reproducible example:   
library(velox)
library(raster)

## Make VeloxRaster
mat <- matrix(1:100, 10, 10)
extent <- c(0,1,0,1)
vx <- velox(mat, extent=extent, res=c(0.1,0.1), crs="+proj=longlat +datum=WGS84 +no_defs")

## Make SpatialPolygonsDataFrame
library(sp)
library(rgeos)
set.seed(0)
coords <- cbind(runif(10, extent[1], extent[2]), runif(10, extent[3], extent[4]))
sp <- SpatialPoints(coords)

# Default example
# https://cran.r-project.org/web/packages/velox/README.html
spol_norm <- gBuffer(sp, width=0.2, byid=TRUE)
spdf_norm <- SpatialPolygonsDataFrame(spol_norm, data.frame(id=1:length(spol_norm)), FALSE)

# Smaller buffer
spol_small<- gBuffer(sp, width=0.05, byid=TRUE)
spdf_small <- SpatialPolygonsDataFrame(spol_small, data.frame(id=1:length(spol_small)), FALSE)

plot(spdf_norm); par(new=F)
plot(spdf_small)

## Extract values and calculate mean, see results
(ex.mat.norm <- vx$extract(spdf_norm, fun=median))
(ex.mat.small <- vx$extract(spdf_small, fun=median)) # -> 3 NA's

EDIT: Here is the workaround solution by using the raster package for extracting data for the NA polygons and then replacing the NA's with the extracted. When using tons of data, like I do, this solution can speed up the processing quite significantly. Continuing from previous code block:

# Rasterize matrix
r <- raster(mat)
# Use polygon center coordinates (SpatialPoints)
ex.mat.small.temp <- extract(r,sp)
ex.mat.small.temp # check
# Compare
ex.mat.small == ex.mat.small.temp
# Replace NA's with ex.mat.small.temp
require(purrr)
ex.mat.small[,1] <- map2_dbl(ex.mat.small,ex.mat.small.temp, coalesce)

# Fixed outcome
ex.mat.small

Best Answer

velox maintainer here.

The development version of velox (available on github) now implements a small option for the VeloxRaster_extract method. If small = TRUE, velox will return raster values for all polygons that intersect with the raster extent. Specifically, for all small or oddly shaped polygons that do not intersect with any cell center, velox will look for intersections with entire cell boxes.

Using the development version of velox and setting the small option to TRUE, your example now yields no NA values:

# Install latest version of velox
library(devtools)
install_github("hunzikp/velox")

# Reproducible example:   
library(velox)
library(raster)

## Make VeloxRaster
mat <- matrix(1:100, 10, 10)
extent <- c(0,1,0,1)
vx <- velox(mat, extent=extent, res=c(0.1,0.1), crs="+proj=longlat     +datum=WGS84 +no_defs")

## Make SpatialPolygonsDataFrame
library(sp)
library(rgeos)
set.seed(0)
coords <- cbind(runif(10, extent[1], extent[2]), runif(10, extent[3], extent[4]))
sp <- SpatialPoints(coords)

# Default example
# from https://cran.r-project.org/web/packages/velox/README.html
spol_norm <- gBuffer(sp, width=0.2, byid=TRUE)
spdf_norm <- SpatialPolygonsDataFrame(spol_norm,     data.frame(id=1:length(spol_norm)), FALSE)

# Smaller buffer
spol_small<- gBuffer(sp, width=0.05, byid=TRUE)
spdf_small <- SpatialPolygonsDataFrame(spol_small, data.frame(id=1:length(spol_small)), FALSE)

plot(spdf_norm); par(new=F)
plot(spdf_small)

## Extract values and calculate mean, see results
(ex.mat.norm <- vx$extract(spdf_norm, fun=median))
(ex.mat.small <- vx$extract(spdf_small, fun=median, small = TRUE)) # -> No NA values anymore.
Related Question