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.
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, extent), runif(10, extent, extent)) 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