此流程是否提取多边形层内所有像素 > 2 的面积，对吗？

Question

目标：提取多边形内值 > 2 的像素面积（以 km2 为单位）。值未反映真实的国家/地区。

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

Answer 1

你声明你

need a CRS that is valid for multiple countries in the Americas and allows me to estimate area in km2- epsg:3857

似乎是基于一种常见的误解，即您不能使用 longitude/latitude 数据来确定区域大小（here 是一些讨论）。

事实上，longitude/latitude是一个很好的测量面积的坐标参考系统。您可以使用一些投影（平面坐标参考系统），但大多数投影都会扭曲区域。因此，如果要使用一个，则需要使用等积投影（例如圆柱等积投影）。

不要使用墨卡托投影（"+proj=merc +a=6378137 +b=6378137，epsg:3857）。墨卡托保留形状，这就是它用于网络制图的原因。它还使格陵兰岛比非洲大；你不能用它来计算面积。更多讨论

一般最好不要投影光栅数据（会有质量损失）。所以这里有一些非常相似的工作流程可以避免这种情况。首先使用 terra，然后使用 raster 和 exactextractr 来计算您想要的内容。

示例数据

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)

计算每个单元格的面积并结合使用的值

a <- cellSize(r, unit="km")
ra <- c(r, a)
names(ra) <- c("values", "area")

提取、子集和计算总和

e <- extract(ra, p, exact=TRUE)
e <- e[e$values>2, ]
sum(e$area * e$fraction)
# [1] 44069.83

或者

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

使用 raster 你可以做类似的事情

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

显式调用 exactextractr

ext <- exactextractr::exact_extract(s, ps) 
ext <- ext[[1]]
ext <- ext[ext$values > 2, ] 
sum(ext$area * ext$coverage_fraction)
#[1] 44056.61

这是一个很好的方法，您可以在其中使用 exactextractr

w <- rr > 2
ext <- exactextractr::exact_extract(aa, ps, weights=w, fun="weighted_sum") 
ext
# [1] 44056.61

此流程是否提取多边形层内所有像素 > 2 的面积，对吗？

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

r

area

r-raster