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.
Without -srcnodata or changing it to -3.39999999999999996e+38 didn't help and yielded the same faulty mosaic with data gaps.
However, changing the -srcnodata to "-3.3999999521443642e+38" (as displayed in QGIS) worked, and now I can obtain a proper mosaic. I understand that this is just a workaround and rather should be fixed by converting the source data. Thank you @user30184 !
Best Answer
Maybe this would be the complete work-flow (as partially mentioned already by @Aaron and @user2856):
gdal:rasterize_over
) to overwrite the raster's values with the values from your shapefile