R土地重叠向量和返回土地面积

Question

我遵循terra::intersect

的示例。

library(terra)
f <- system.file("ex/lux.shp", package="terra")
v <- vect(f)
e <- ext(5.6, 6, 49.55, 49.7)
x <- intersect(v, e)
p <- vect(c("POLYGON ((5.8 49.8, 6 49.9, 6.15 49.8, 6 49.6, 5.8 49.8))",
"POLYGON ((6.3 49.9, 6.2 49.7, 6.3 49.6, 6.5 49.8, 6.3 49.9))"), crs=crs(v))
values(p) <- data.frame(pid=1:2, area=expanse(p))
y <- intersect(v, p)

和我最终想要它总结的是通过p和v连接创建的所有多边形的土地面积。我的意思是，对于p中的每个多边形，在Diekirch, Grevenmacher和Luxembourg中有多少公顷的多边形。

我试试这个:

lapply(y, FUN=expanse, unit="ha")
expanse(v)
yy <- union(v,p)
lapply(yy, FUN=expanse, unit="ha")

，但在这两种情况下，宽阔返回与

相同的值

expanse(v)

Answer 1

你可以做

e <- expanse(y, unit="km")
aggregate(cbind(area=e), y[,c("pid", "NAME_1"), drop=TRUE], sum)
#  pid       NAME_1       area
#1   1     Diekirch 226.271445
#2   2     Diekirch   8.147181
#3   2 Grevenmacher 286.019968
#4   1   Luxembourg 194.340843
#5   2   Luxembourg  32.281530

或者像这样(我更喜欢上面的方法，因为它避免了聚合多边形的需要)

z <- aggregate(y, c("pid", "NAME_1"))
cbind(z[, c("pid", "NAME_1"), drop=TRUE], 
area=expanse(z, unit="km"))
#  pid       NAME_1       area
#1   1     Diekirch 226.271445
#2   1   Luxembourg 194.340843
#3   2     Diekirch   8.147181
#4   2 Grevenmacher 286.019968
#5   2   Luxembourg  32.281530

Answer 2

以下是sf和tidyverse的方法:

library(terra)
library(tidyverse)
library(sf)
f <- system.file("ex/lux.shp", package="terra")
v <- vect(f)
e <- ext(5.6, 6, 49.55, 49.7)
x <- intersect(v, e)
p <- vect(c("POLYGON ((5.8 49.8, 6 49.9, 6.15 49.8, 6 49.6, 5.8 49.8))",
"POLYGON ((6.3 49.9, 6.2 49.7, 6.3 49.6, 6.5 49.8, 6.3 49.9))"), crs=crs(v))
values(p) <- data.frame(pid=1:2, area=expanse(p))
y <- terra::intersect(v, p)
# sanity check: view data, cropping polygons, and cropped layer
plot(v)
plot(p, add = TRUE, col = "red")
plot(y, add = TRUE, col = "blue")
# convert to sf (data frame), mutate an area, and sum area per name group
st_as_sf(y) %>% 
mutate(area = st_area(geometry)) %>% 
group_by(NAME_1) %>% 
summarise(area_sum = sum(area)) %>% 
ungroup()

返回每个名称组的面积总和:

Simple feature collection with 3 features and 2 fields
Geometry type: GEOMETRY
Dimension:     XY
Bounding box:  xmin: 5.8 ymin: 49.6 xmax: 6.5 ymax: 49.9
Geodetic CRS:  WGS 84
# A tibble: 3 × 3
NAME_1         area_sum                             geometry
<chr>             [m^2]                       <GEOMETRY [°]>
1 Diekirch     233635561. MULTIPOLYGON (((6.285424 49.87085, …
2 Grevenmacher 285068244. POLYGON ((6.42977 49.72977, 6.5 49.…
3 Luxembourg   225869770. MULTIPOLYGON (((6.217318 49.68268, …