0

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 ?raster::extract)

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)
2
  • Both the raster and vector need to be in a distance base projection (eg., meters). The developer of raster contributes to this forum so, he will likely chime in here if there is an actual package issue. However, in the meantime I would recommend calculating a buffer using rgeos::rBuffer (for sp objects) or sf::st_buffer (for sf objects) then using exactextract::exact_extract for returning the mean raster values. It will be a massive speed improvement and offers the option of a weighted mean using the fractional intersection of cells (ie., area weighted). May 14 at 14:11
  • Dear Jeffrey, thank you very much for your answer. exact_extract() seems to provide the behavior that I was expecting. I included this suggestion in the main body of the question. May 17 at 12:43
2

Adding a buffer may not affect the result if the buffer is small, relative to the cell size of the raster.

You provide no reproducible example, or information on the data your using such as show(raster) or head(points) making it very difficult to help you. Can you please edit your question and include some of that?

You seem to be barking up the tree with changing the coordinate reference system (CRS). raster::extract should work irrespective of whether the CRS is planar or angular (lonlat). You appear to be setting a CRS to your data on the basis of what you want it to be. That is not valid. If need be, it should only be set to what it is.

(You can put a post-it with "Rolls-Royce" on your bicycle --- but don't expect a chauffeur to appear and take you for a spin.)

See spTransform to change a crs.

You also say that "the units are in degrees and as explained in the extract function help file, this does not make the buffer option to work". That is not correct. Where does it say that?

1
  • Dear Robert, thank you very much for your answer. As suggested, I included a reproducible example and the information of my data. Regarding the point about the CRS you are of course right. My point was that maybe I was expressing the number in the buffer option incorrectly since they were in degrees, like you stated, not adding enough buffer space to intersect the pixels with values. May 17 at 12:46

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.