What is the map unit specified by my CRS?

Solution for What is the map unit specified by my CRS?
is Given Below:

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], [email protected], 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!

You need to look at the grouping of the WKT:

#> 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]]]]

There are 3 *UNIT blocks in there:

DATUM["World Geodetic System 1984",
      ELLIPSOID["WGS 84",6378137,298.257223563,
      LENGTHUNIT["metre",1]]

and 2 Axis:

AXIS["longitude",east,
        ORDER[1],
        ANGLEUNIT["degree",0.0174532925199433,
            ID["EPSG",9122]]],

The first applies to the Ellipsoid (and is in metres) the other 2 are the ones you are interested in as they define how your data is stored (and presented) and are in degrees (as is expected with an unprojected WGS CRS). So your buffer will be in degrees and as you say will be absurd, provided nothing breaks it will return all of your data.

To work in metres you will need to reproject your data to a planar local CRS (preferably one that is EquiDistant for your dataset).

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).