Goal: Extract the area in km2 of pixels that have values > 2 within a polygon. Values are not reflecting real country areas.
library(sf)
library(raster)
library(exactextractr)
# Generate raster
r <- raster::raster(matrix(0:7, ncol=10), xmn=0, ymn=0, xmx=10, ymx=10)
poly <- sf::st_as_sfc('POLYGON ((2 2, 7 6, 4 9, 2 2))')
#In my case I need a CRS that is valid for multiple countries in the Americas and allows me to estimate area in km2- epsg:3857
crs(r) <- "+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext +no_defs"
poly<-st_set_crs(poly,"+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext +no_defs")
#Extract area of pixels that have values > 2. This is in particular what I'm interested in, is my function argument doing what I say it does.
ext<-exact_extract(r,poly,function(values, coverage_fraction)
length(values > 2)) #6 values
#Determine pixel size
res(r) #1 10
res.m2<-10
res.km2<-res.m2/1000000
#Determine area in km2:multiply number of pixels >2 by the pixel area
tot.area<-res.km2*ext
Your statement that you
need a CRS that is valid for multiple countries in the Americas and allows me to estimate area in km2- epsg:3857
seems based on the common misconception that you cannot use longitude/latitude data to determine area sizes (here is some discussion).
In fact, longitude/latitude is a great coordinate reference system to measure area. You can use some projections (planar coordinate reference systems), but most projections distort area. So if you were to use one, you would need to use an equal-area projection (e.g. cylindrical equal-area).
Do not use the Mercator projection ("+proj=merc +a=6378137 +b=6378137, epsg:3857). Mercator conserves shape, and that is why it is used in web-mapping. It also makes Greenland larger than Africa; and you cannot use it to compute area. More discussion here
It is generally best to not project raster data (there is quality loss). So here are some very similar work-flows that avoid that. First with terra and then with raster and exactextractr that compute what you are after.
Example data
library(terra)
p <- vect('POLYGON ((2 2, 7 6, 4 9, 2 2))')
r <- rast(nrows=10, ncols=10, xmin=0, ymin=0, xmax=10, ymax=10)
r <- init(r, -2:7)
Compute area of each cell and combine with the values used
a <- cellSize(r, unit="km")
ra <- c(r, a)
names(ra) <- c("values", "area")
Extract, subset, and computed sum
e <- extract(ra, p, exact=TRUE)
e <- e[e$values>2, ]
sum(e$area * e$fraction)
# [1] 44069.83
Alternatively
x <- ifel(r>2, r, NA)
a <- cellSize(r, unit="km")
ax <- mask(a, x)
ee <- extract(ax, p, exact=TRUE)
sum(ee$area * ee$fraction, na.rm=TRUE)
#[1] 44069.83
With raster you can do something similar
library(raster)
rr <- raster(nrows=10, ncols=10, xmn=0, ymn=0, xmx=10, ymx=10)
values(rr) <- rep(-2:7, 10)
ps <- sf::st_as_sfc('POLYGON ((2 2, 7 6, 4 9, 2 2))')
ps <- as(ps, "Spatial")
crs(ps) <- crs(rr)
aa <- area(rr)
s <- stack(aa, rr)
names(s) <- c("area", "values")
v <- extract(s, ps, exact=TRUE, weights=TRUE, normalizeWeights=FALSE)
v <- as.data.frame(v[[1]])
v <- v[v$values > 2, ]
sum(v$area * v$weight)
# [1] 44056.61
Explicitly calling exactextractr
ext <- exactextractr::exact_extract(s, ps)
ext <- ext[[1]]
ext <- ext[ext$values > 2, ]
sum(ext$area * ext$coverage_fraction)
#[1] 44056.61
Here is a nice way in which you can use exactextractr
w <- rr > 2
ext <- exactextractr::exact_extract(aa, ps, weights=w, fun="weighted_sum")
ext
# [1] 44056.61