I'm trying to extract a mean value from multiple pixels of raster file "r" for a vector shapefile "v", by using the r extract function (see
Since some of the file "v" shapes falls outside the raster map (e.g. cities on a coast), and given that the NA values are already transformed to zeros, I extract by using a buffer, so that the latter would still get a reasonable value as a mean.
Assuming a buffer of 1000 meters, the code is then:
extracted_values <- extract(r, v, fun=mean, buffer = 1000)
EDIT: Following Robert advice I add a reproducible example that mimics my data situation, and also compares with the solution proposed by Jeffrey:
# Generate a raster and change the crs to the ones of my data r <- raster(ncol=1000, nrow=1000, xmx=-80, xmn=-150, ymn=20, ymx=60) values(r) <- runif(ncell(r)) r <- projectRaster(r, crs = "+proj=longlat +a=6367470 +b=6367470 +no_defs ",method = "bilinear") # Set a part of the raster equal to zero r[1:800000] = 0 # generate sample polygon crdref <- CRS('+proj=longlat +datum=WGS84') lon <- c(-116.8, -114.2, -112.9, -111.9, -114.2, -115.4, -117.7) lat <- c(41.3, 42.9, 42.4, 39.8, 37.6, 38.3, 37.6) lonlat <- cbind(lon, lat) v <- spPolygons(lonlat, crs=crdref) v <- spTransform(x = v, CRSobj = crs(r)) # Now these represents the starting data situation # convert to 3857 to have distances in meters r <- projectRaster(r, crs = "+init=epsg:3857") v <- spTransform(x = v, CRSobj = crs(r)) # I add a "buffered" version to confirm visually that now the buffered polygon intersects the raster pixels >0 v_buf <- gBuffer(v,width=1500000) plot(r) plot(v, add = TRUE) plot(v_buf, add=TRUE) # Results of extract: no_buf <- extract(r, v, fun=mean) print(no_buf) buf <- extract(r, v, fun=mean, buffer=1500000) print(buf) buf_higher <- extract(r, v, fun=mean, buffer=1600000) print(buf_higher) buf_weights <- extract(r, v, fun=mean, buffer=1500000, weights= T, normalizeWeights = T) print(buf_weights) # with extract_exact and using the v_buf polygon (already with buffer) buf_exact <- exact_extract(r, v_buf, 'weighted_mean', weights = 'area') print(buf_exact) # Does it change increasing the buffer? v_buf2 <- gBuffer(v,width=1600000) buf_exact2 <- exact_extract(r, v_buf2, 'weighted_mean', weights = 'area') print(buf_exact2)
It seems however that the buffer option is not affecting the final result, whether by excluding it, or by increasing/decreasing the number of meters.
EDIT: I add here the result of show(r), edited for privacy in the source and names fields. In any case the original raster file is taken from a GRIB file from the "ERA5-Land hourly data from 1981 to present":
class : RasterLayer band : 1 (of 156 bands) dimensions : 1801, 3600, 6483600 (nrow, ncol, ncell) resolution : 0.1, 0.1 (x, y) extent : -0.05, 359.95, -90.05, 90.05 (xmin, xmax, ymin, ymax) crs : +proj=longlat +a=6367470 +b=6367470 +no_defs source : FILENAME.grib names : NAMES
The reason of the issue could probably be related to the raster coordinate system. Originally the file r is CRS arguments: +proj=longlat +a=6367470 +b=6367470 +no_defs (notice that I always change the crs of v to match the one of the raster before extracting; the original crs of v is +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 ). If I understood correctly, the units are in degrees and as explained in the extract function help file, this does not make the buffer option to work. I tried then to reproject the raster r to crs = "+init=epsg:3857", which (I could be wrong) should be in meters, but again, varying the buffer option number seems not to affect the extracted mean value.
EDIT: I cancelled this paragraph because as Robert correctly pointed out, there is no such affirmation in the help file (probably I got confused with one of the other answers found online, I apologize for that and I will be more precise next time). In any case I observed the issue both with the original CRS and both with the reprojected one in meters.
So I want to ask: why the buffer option is not affecting the mean result for the extract function?
And, more in general, if I want to express the buffer in meters which crs would be best suited for the task? EDIT: cancelled this last part because not specific enough as it probably depends on the problem's characteristics such at which scale one is operating.
EDIT: following the advice by Jeffrey the procedure works, in the sense that increasing the buffer causes the extracted mean value to change (as correctly, the polygons with buffer extra space now cover pixels with values different than zero). Here is the revised code (notice that here the area is weighted):
library(exactextractr) library(rgeos) r <- projectRaster(r, crs = "+init=epsg:3857") v <- spTransform(x = v, CRSobj = crs(r)) v <- gBuffer(v,width=1000, byid = T) extracted_values <- exact_extract(r, v, 'weighted_mean', weights = 'area', force_df = T, append_cols = T)