Is this flow to extract the area of all pixels > 2 within a polygon layer, correct?

user_5329164 0
Viewed 143

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
r-raster area r
1 Answers

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
user_6352450
Related