r语言 - 此流程是否用于提取多边形图层内所有像素> 2 的面积，正确吗？

Question

目标：提取具有>2。价值观并不能反映真实的乡村地区。

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

您的声明

需要一个对美洲，并允许我以平方公里为单位估计面积-epsg:3857

似乎基于一种常见的误解，即不能使用经度/纬度数据来确定区域大小(这里有一些讨论(。

事实上，经度/纬度是测量面积的一个很好的坐标参考系。您可以使用一些投影(平面坐标参考系(，但大多数投影会扭曲区域。因此，如果你要使用一个，你需要使用等面积投影(例如圆柱形等面积投影(。

不要使用墨卡托投影("+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