我在地图上有一组空间坐标，以及穿过它们的多线串。我需要计算出线条在每种颜色的多边形内所花费的长度。
caveats:
最小代表：
poly1 <-
# create list of matrices and the first point same as last point
list(
matrix(
c(0, 0,
4, 0,
5, 1,
4, 2,
3, 2,
1, 1,
0, 0),
ncol=2, byrow=T
)
)
poly1 <- sf::st_polygon(poly1)
poly2 <-
# create list of matrices and the first point same as last point
list(
matrix(
c(4, 1,
7, 0,
7, 1,
6, 3,
3, 2,
2, 1,
4, 1),
ncol=2, byrow=T
)
)
poly2 <- sf::st_polygon(poly2)
poly3 <-
# create list of matrices and the first point same as last point
list(
matrix(
c(7, 1,
10, 1,
12, 2,
11, 4,
8, 3,
7, 2,
7, 1),
ncol=2, byrow=T
)
)
poly3 <- sf::st_polygon(poly3)
line <-
# create list of matrices and the first point same as last point
list(
matrix(
c(0, 1,
2, 0,
4, 1,
6, 3,
8, 2,
10, 1,
12, 1),
ncol=2, byrow=T
)
)
line <-
sf::st_multilinestring(line)
ggplot() +
geom_sf(data = poly1, fill = "green",alpha=.5) +
geom_sf(data = poly2, fill = "blue",alpha=.5) +
geom_sf(data = poly3, fill = "green", alpha=.5)+
geom_sf(data = line, color="black", size=2) +
ggthemes::theme_map()
期望的输出：
poly1_cropped <- st_intersection( line, poly1)
poly2_cropped <- st_intersection( line, poly2)
poly3_cropped <- st_intersection( line, poly3)
ggplot() +
geom_sf(data = poly1, fill = "green",alpha=.5) +
geom_sf(data = poly2, fill = "blue",alpha=.5) +
geom_sf(data = poly3, fill = "green", alpha=.5)+
geom_sf(data = line, color="black", size=2) +
ggthemes::theme_map() +
geom_sf(data = poly1_cropped, color="green", size=2) +
geom_sf(data = poly2_cropped, color="blue", size=2) +
geom_sf(data = poly3_cropped, color="green", size=2)
然后是一个数据框，用于量化线条何时穿过每个形状，如下所示：
shape | color | shapes_overlapping | length
poly1 green 0 3
poly1 green 1 .5
poly2 blue 1 .5
poly2 blue 0 2
poly3 green 0 2.5
我有一个适用于您的示例的解决方案，您希望可以根据实际数据进行修改。
我做的第一件事就是结合所有 3 个POLYGON物体变成Simple feature collection，最初有一个MULTIPOLYGON但然后分成每个单独的POLYGON using st_cast().
# combine polygons into a single simple feature collection
c(poly1, poly2, poly3) %>%
# make into simple feature
st_sfc %>%
st_sf %>%
# split into individual polygons
st_cast('POLYGON') %>%
{. ->> int1}
int1
# Simple feature collection with 3 features and 0 fields
# Geometry type: POLYGON
# Dimension: XY
# Bounding box: xmin: 0 ymin: 0 xmax: 12 ymax: 4
# CRS: NA
# geometry
# 1 POLYGON ((0 0, 4 0, 5 1, 4 ...
# 2 POLYGON ((4 1, 7 0, 7 1, 6 ...
# 3 POLYGON ((7 1, 10 1, 12 2, ...
接下来，查找是否有任何多边形重叠，如果重叠，则每个重叠处有多少层。我们用st_intersection()在这里，但与您找到两个几何图形之间的交集的示例不同，我们将其应用于单个几何图形sf对象，返回它的自交集，包括n.overlaps（层数）和origins（重叠区域中的原始多边形）。更多详细信息可以在以下位置找到sf page here。我们还做了一个新的id对于每个多边形和层数来说都是唯一的列。
# now, we need to find if any polygons overlap, and how many layers there are
int1 %>%
st_intersection %>%
filter(
st_geometry_type(.) == 'POLYGON'
) %>%
mutate(
id = row_number()
) %>%
{. ->> int2}
int2
# Simple feature collection with 4 features and 3 fields
# Geometry type: POLYGON
# Dimension: XY
# Bounding box: xmin: 0 ymin: 0 xmax: 12 ymax: 4
# CRS: NA
# n.overlaps origins geometry id
# 1 1 1 POLYGON ((4 0, 0 0, 1 1, 3 ... 1
# 1.1 2 1, 2 POLYGON ((4 2, 5 1, 4.75 0.... 2
# 2 1 2 POLYGON ((7 1, 7 0, 4.75 0.... 3
# 3 1 3 POLYGON ((7 2, 8 3, 11 4, 1... 4
为了证明这一点，我们可以绘制int2和颜色（fill）每个多边形由我们id column.
int2_plot <- int2 %>%
ggplot()+
geom_sf(aes(fill = factor(id)))+
geom_sf(data = line)
int2_plot
然后，我们使用ggplot2::ggplot_build()提取绘图的组成部分 - 特别是与每个多边形对应的颜色id。我们用plotrix::color.id()将十六进制颜色代码转换为颜色名称。这可能需要也可能不需要，具体取决于您如何解释颜色。
# now, get colours of each polygon from ggplot
int2_plot %>%
ggplot_build %>%
.$data %>%
.[[1]] %>%
tibble %>%
select(fill, group, id = group) %>%
rowwise %>%
mutate(
colour = plotrix::color.id(fill),
) %>%
select(id, colour) %>%
{. ->> plot_colours}
plot_colours
# # A tibble: 4 x 2
# # Rowwise:
# id colour
# <int> <chr>
# 1 1 salmon
# 2 2 chartreuse3
# 3 3 turquoise3
# 4 4 mediumpurple1
现在我们找到长度line与每个多边形相交。我们用st_intersection()找到LINESTRING与多边形重叠的对象，然后st_length()来计算它的长度。因为这个示例数据没有CRS，这个长度是无单位的，并作为一个返回numeric向量。如果您的数据有CRS, st_length()将返回一个units具有数字和距离单位的向量，例如3.726 [m];可以使用以下方式检索该号码as.numeric().
# now, what's the length of line intersecting each polygon (colour)?
int2 %>%
mutate(
line_overlap = st_intersection(int2, line) %>% st_length
) %>%
{. ->> int3}
int3
# Simple feature collection with 4 features and 4 fields
# Geometry type: POLYGON
# Dimension: XY
# Bounding box: xmin: 0 ymin: 0 xmax: 12 ymax: 4
# CRS: NA
# n.overlaps origins geometry id line_overlap
# 1 1 1 POLYGON ((4 0, 0 0, 1 1, 3 ... 1 3.7267800
# 1.1 2 1, 2 POLYGON ((4 2, 5 1, 4.75 0.... 2 0.7071068
# 2 1 2 POLYGON ((7 1, 7 0, 4.75 0.... 3 2.1213203
# 3 1 3 POLYGON ((7 2, 8 3, 11 4, 1... 4 2.9814240
目前，我们有一行用于重叠区域，其中包含重叠中涉及的多边形的向量（origins）。为了为每个多边形创建一个单独的行，我们将origins to a list，然后可以使用以下方法将每个列表元素（原始多边形）分成单独的行unnest().
# now, duplicate rows where n.overlaps > 1 - so there's a row for each polygon in the overlap
int3 %>%
rowwise %>%
mutate(
origins = list(origins)
) %>%
tidyr::unnest(cols = c('origins')) %>%
{. ->> int4}
int4
# Simple feature collection with 5 features and 4 fields
# Geometry type: POLYGON
# Dimension: XY
# Bounding box: xmin: 0 ymin: 0 xmax: 12 ymax: 4
# CRS: NA
# # A tibble: 5 x 5
# n.overlaps origins geometry id line_overlap
# <int> <int> <POLYGON> <int> <dbl>
# 1 1 1 ((4 0, 0 0, 1 1, 3 2, 2 1, 4 1, 4.75 0.75, 4 0)) 1 3.73
# 2 2 1 ((4 2, 5 1, 4.75 0.75, 4 1, 2 1, 3 2, 4 2)) 2 0.707
# 3 2 2 ((4 2, 5 1, 4.75 0.75, 4 1, 2 1, 3 2, 4 2)) 2 0.707
# 4 1 2 ((7 1, 7 0, 4.75 0.75, 5 1, 4 2, 3 2, 6 3, 7 1)) 3 2.12
# 5 1 3 ((7 2, 8 3, 11 4, 12 2, 10 1, 7 1, 7 2)) 4 2.98
最后，我们加入之前拉出的每个多边形的颜色。
# now, join in the colours of each polygon
int4 %>%
left_join(plot_colours) %>%
{. ->> int5}
int5
# Simple feature collection with 5 features and 5 fields
# Geometry type: POLYGON
# Dimension: XY
# Bounding box: xmin: 0 ymin: 0 xmax: 12 ymax: 4
# CRS: NA
# # A tibble: 5 x 6
# n.overlaps origins geometry id line_overlap colour
# <int> <int> <POLYGON> <int> <dbl> <chr>
# 1 1 1 ((4 0, 0 0, 1 1, 3 2, 2 1, 4 1, 4.75 0.75, 4 0)) 1 3.73 salmon
# 2 2 1 ((4 2, 5 1, 4.75 0.75, 4 1, 2 1, 3 2, 4 2)) 2 0.707 chartreuse3
# 3 2 2 ((4 2, 5 1, 4.75 0.75, 4 1, 2 1, 3 2, 4 2)) 2 0.707 chartreuse3
# 4 1 2 ((7 1, 7 0, 4.75 0.75, 5 1, 4 2, 3 2, 6 3, 7 1)) 3 2.12 turquoise3
# 5 1 3 ((7 2, 8 3, 11 4, 12 2, 10 1, 7 1, 7 2)) 4 2.98 mediumpurple1
如果您不希望结果为sf对象，我们可以使用st_drop_geometry()仅返回属性组件。
# keep the attribute data only
int5 %>%
st_drop_geometry
# # A tibble: 5 x 5
# n.overlaps origins id line_overlap colour
# * <int> <int> <int> <dbl> <chr>
# 1 1 1 1 3.73 salmon
# 2 2 1 2 0.707 chartreuse3
# 3 2 2 2 0.707 chartreuse3
# 4 1 2 3 2.12 turquoise3
# 5 1 3 4 2.98 mediumpurple1