I am rather new to spatial work with R and recently ran into the following issue:
I have gridded climate data (maximum temperature, monthly) and want to extract observations at some locations, using a circular buffer with a radius of 20KM. Let me first set it up and plot it...
library(raster)
# load NetCDF file into SpatialRaster
ncin <- raster::brick(x='cru_ts4.05.1901.2020.tmx.dat.nc')
#> Warning in .varName(nc, varname, warn = warn): varname used is: tmx
#> If that is not correct, you can set it to one of: tmx, stn
# get shape file
poly <- readRDS("C:/Users/m/shapes/gadm36_NGA_1_sf.rds") %>%
st_transform(proj4string(ncin))
ncin <- raster::crop(ncin, poly, snap = "out")
#clip raster outside the poly border
# result is a raster brick with dimensions (lat=20, lon=25, t=1440)
load ("C:/Users/m/NGA_panel.Rda") # load household data, with coordinates
# plot households as dots on the grid; see whether all parts display correctly!
plot(ncin$X2020.12.16, main="HH-locations")
# plot(poly$geometry, add = TRUE)
points(NGA_panel$LON_DD_MOD, NGA_panel$LAT_DD_MOD, pch=19, cex=0.5, col =1)
Here is the resulting graph. The black dots represent households. For each of them, I want to extract a value from the grid, using a circular buffer with radius 20 KM. The value extracted is the mean of raster values within the buffer.
My CRS is +proj=longlat +datum=WGS84 +no_defs, which according to the documentation of raster::buffer() does not necessarily have meters as its default map unit. st_crs(ncin) yields the following details:
st_crs(ncin) # does the current crs take meters as length-unit?
#> Coordinate Reference System:
#> User input: +proj=longlat +datum=WGS84 +no_defs
#> wkt:
#> GEOGCRS["unknown",
#> DATUM["World Geodetic System 1984",
#> ELLIPSOID["WGS 84",6378137,298.257223563,
#> LENGTHUNIT["metre",1]],
#> ID["EPSG",6326]],
#> PRIMEM["Greenwich",0,
#> ANGLEUNIT["degree",0.0174532925199433],
#> ID["EPSG",8901]],
#> CS[ellipsoidal,2],
#> AXIS["longitude",east,
#> ORDER[1],
#> ANGLEUNIT["degree",0.0174532925199433,
#> ID["EPSG",9122]]],
#> AXIS["latitude",north,
#> ORDER[2],
#> ANGLEUNIT["degree",0.0174532925199433,
#> ID["EPSG",9122]]]]
Does LENGTHUNIT["metre",1] above mean that a buffer radius will correctly be scaled in meters? Continuing with the extraction gives me...
NGA_spatial <- SpatialPointsDataFrame(
NGA_panel[108:107], proj4string=ncin@crs, NGA_panel) # spatialize dataset (i.e. make SPDF)
NGA_spatial <- NGA_spatial[sample(nrow(NGA_spatial), 100), ] #### LIMIT to 100 random rows FOR TESTING!
# extract values for households, using a circular buffer with radius 20KM:
data <- raster::extract(ncin,
NGA_spatial, # SPDF with centroids for buffer
buffer = 20000, # buffer size, units depend on CRS (meters or degrees?)
fun=mean, # what value to extract
df=TRUE) # create a data frame?
summary(data[sample(ncol(data), 3)]) # summary statistics of 3 random columns
#> X2005.09.16 X2000.02.15 X2001.02.15
#> Min. :28.80 Min. :29.20 Min. :31.20
#> 1st Qu.:29.77 1st Qu.:30.30 1st Qu.:32.17
#> Median :30.70 Median :33.30 Median :34.10
#> Mean :31.21 Mean :32.46 Mean :33.63
#> 3rd Qu.:32.25 3rd Qu.:33.80 3rd Qu.:34.55
#> Max. :37.10 Max. :35.00 Max. :36.00
Clearly there is some variation in my extracted values, but I am uncertain whether this is extracted by a buffer of 20000 meters (as I want it to be) or a buffer of 20000 degrees (which would be absurd).
Can someone perhaps tell me whether this correctly extracts in meters and whether there is a way to control for that (e.g. by plotting the buffer somehow)? Thanks a lot!
The documentation of ?raster::extract for argument "buffer" states that
If the data are not projected (latitude/longitude), the unit should be meters. Otherwise it should be in map-units (typically also meters).
+proj=longlat +datum=WGS84 is clearly (?) latitude/longitude, so you should use meters.
This is illustrated in below, taken (with minor changes) from the examples in ?extract:
library(raster)
r <- raster(ncols=36, nrows=18, vals=1:(18*36))
xy <- cbind(0, seq(20, 80, by=20))
xy
# [,1] [,2]
#[1,] 0 20
#[2,] 0 40
#[3,] 0 60
#[4,] 0 80
extract(r, xy, buffer=1000000)
#[[1]]
#[1] 234 235 270 271
#[[2]]
#[1] 162 163 198 199
#[[3]]
#[1] 89 90 91 92 126 127
#[[4]]
# [1] 13 14 15 16 17 18 19 20 21 22 23 24 51 52 53 54 55 56 57 58
As you can see, you get increasingly more values (i.e., cells within the buffer) as you move pole-wards. This is as expected as the longitudinal distance decreases, roughly proportional to the cos(latitude).