问题描述
目标：提取多边形内值 > 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
解决方法
你的声明
需要一个对多个国家/地区有效的 CRS 美洲并允许我以 km2 为单位估计面积- epsg:3857
似乎是基于一种常见的误解，即您不能使用经度/纬度数据来确定区域大小（here 是一些讨论）。
事实上，经度/纬度是测量面积的绝佳坐标参考系统。您可以使用一些投影（平面坐标参考系统），但大多数投影会扭曲区域。因此，如果要使用一个，则需要使用等积投影（例如圆柱等积）。
不要不要使用墨卡托投影（"+proj=merc +a=6378137 +b=6378137，epsg:3857）。墨卡托保存形状，这就是它用于网络制图的原因。它还使格陵兰岛比非洲大；你不能用它来计算面积。更多讨论here
通常最好不要投影栅格数据（存在质量损失）。所以这里有一些非常相似的工作流程可以避免这种情况。首先使用 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,exact=TRUE)
sum(ee$area * ee$fraction,na.rm=TRUE)
#[1] 44069.83
使用 raster 你可以做类似的事情
library(raster)
rr <- raster(nrows=10,xmn=0,ymn=0,xmx=10,ymx=10)
values(rr) <- rep(-2:7,10)
ps <- sf::st_as_sfc('POLYGON ((2 2,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,weights=w,fun="weighted_sum")
ext
# [1] 44056.61