两个多边形之间重叠？
问题描述 投票：0回答：1
我有两个带有纬度/经度坐标的示例数据集，称为 data1 和 data2，如下所示：
data1
longitude latitude
406700.4 7201767
398192.5 7206899
402995.9 7223166
424210.3 7264295
449045.0 7315044
476329.8 7304172
516543.0 7281898
543151.1 7292038
553048.0 7299822
551770.3 7275758
561142.6 7261667
560292.0 7251473
544158.5 7229430
566519.0 7224846
578221.2 7242260
589978.4 7237221
614511.6 7241832
629119.1 7229805
643408.0 7229532
675159.5 7261319
679089.6 7273392
688517.1 7283675
699653.7 7305478
703623.1 7319123
697325.5 7335428
732159.0 7341021
722027.8 7336020
729827.0 7312423
741539.1 7289712
742604.1 7284240
data2
548055.7 7291742
550953.6 7313012
553370.8 7302842
543710.0 7306752
542931.4 7305144
544490.5 7308543
564527.1 7320696
578368.4 7297569
594972.5 7299350
635337.6 7300907
649867.0 7308755
642136.1 7323236
644287.0 7301646
647413.1 7298133
660147.9 7294922
661158.7 7304247
677788.7 7303612
688676.9 7303413
686350.5 7301771
690122.9 7303681
683197.7 7308678
648936.3 7312211
637416.5 7316094
631286.5 7335698
621747.9 7359180
610528.1 7380706
601504.8 7402010
635005.1 7320464
680882.5 7274095
677057.4 7259603
我已经使用以下代码计算了两个数据集的起始范围并从这些起始范围生成了多边形：
library(ks)
library(dplyr)
#Home range 1
h1 <- Hpi(data1, pilot = "samse", binned = T)
kernPI1 <- kde(data1, H = h1)
cont = contourLevels(kernPI1, cont = 95)
line1 = contourLines(x = kernPI1$eval.points[[1]], y =
kernPI1$eval.points[[2]],
z = kernPI1$estimate,
level = cont)
x1<-line1[[1]][["x"]]
y1<-line1[[1]][["y"]]
dataline1<-data.frame(x1,y1)
polygon1 <- dataline1 %>%
st_as_sf(coords = c("x1", "y1"), crs = 5321) %>%
summarise(geometry = st_combine(geometry)) %>%
st_cast("POLYGON")
# home range 2
h2 <- Hpi(data2, pilot = "samse", binned = T)
kernPI2 <- kde(data2, H = h2)
cont = contourLevels(kernPI2, cont = 95)
line2 = contourLines(x = kernPI2$eval.points[[1]], y =
kernPI2$eval.points[[2]],
z = kernPI2$estimate,
level = cont)
x2<-line2[[1]][["x"]]
y2<-line2[[1]][["y"]]
dataline2<-data.frame(x2,y2)
polygon2 <- dataline2 %>%
st_as_sf(coords = c("x2", "y2"), crs = 5321) %>%
summarise(geometry = st_combine(geometry)) %>%
st_cast("POLYGON")
我想找到这两个多边形之间的重叠面积，可以是 m2 中的数字，也可以是百分比。我如何在 R 中编写此代码？
1个回答
0
投票
投票
让我们使用以下 Python 脚本来查看原始数据，不考虑对其进行任何修正，因为它们具有经度和纬度属性：
import matplotlib.pyplot as plt
# Define the points for data1
data1_points = [
(406700.4, 7201767), (398192.5, 7206899), (402995.9, 7223166),
(424210.3, 7264295), (449045.0, 7315044), (476329.8, 7304172),
(516543.0, 7281898), (543151.1, 7292038), (553048.0, 7299822),
(551770.3, 7275758), (561142.6, 7261667), (560292.0, 7251473),
(544158.5, 7229430), (566519.0, 7224846), (578221.2, 7242260),
(589978.4, 7237221), (614511.6, 7241832), (629119.1, 7229805),
(643408.0, 7229532), (675159.5, 7261319), (679089.6, 7273392),
(688517.1, 7283675), (699653.7, 7305478), (703623.1, 7319123),
(697325.5, 7335428), (732159.0, 7341021), (722027.8, 7336020),
(729827.0, 7312423), (741539.1, 7289712), (742604.1, 7284240)
]
# Define the points for data2
data2_points = [
(548055.7, 7291742), (550953.6, 7313012), (553370.8, 7302842),
(543710.0, 7306752), (542931.4, 7305144), (544490.5, 7308543),
(564527.1, 7320696), (578368.4, 7297569), (594972.5, 7299350),
(635337.6, 7300907), (649867.0, 7308755), (642136.1, 7323236),
(644287.0, 7301646), (647413.1, 7298133), (660147.9, 7294922),
(661158.7, 7304247), (677788.7, 7303612), (688676.9, 7303413),
(686350.5, 7301771), (690122.9, 7303681), (683197.7, 7308678),
(648936.3, 7312211), (637416.5, 7316094), (631286.5, 7335698),
(621747.9, 7359180), (610528.1, 7380706), (601504.8, 7402010),
(635005.1, 7320464), (680882.5, 7274095), (677057.4, 7259603)
]
# Extract x and y coordinates for each dataset
data1_x = [point[0] for point in data1_points]
data1_y = [point[1] for point in data1_points]
data2_x = [point[0] for point in data2_points]
data2_y = [point[1] for point in data2_points]
# Plot the points
plt.figure(figsize=(10, 6))
plt.scatter(data1_x, data1_y, color='blue', label='Data1')
plt.scatter(data2_x, data2_y, color='green', label='Data2')
# Set labels and legend
plt.xlabel('Longitude')
plt.ylabel('Latitude')
plt.title('Scatter Plot of Data Points')
plt.legend()
# Show plot
plt.grid(True)
plt.show()
import matplotlib.pyplot as plt
from scipy.spatial import ConvexHull
import numpy as np
# Define the points for data1
data1_points = np.array([
(406700.4, 7201767), (398192.5, 7206899), (402995.9, 7223166),
(424210.3, 7264295), (449045.0, 7315044), (476329.8, 7304172),
(516543.0, 7281898), (543151.1, 7292038), (553048.0, 7299822),
(551770.3, 7275758), (561142.6, 7261667), (560292.0, 7251473),
(544158.5, 7229430), (566519.0, 7224846), (578221.2, 7242260),
(589978.4, 7237221), (614511.6, 7241832), (629119.1, 7229805),
(643408.0, 7229532), (675159.5, 7261319), (679089.6, 7273392),
(688517.1, 7283675), (699653.7, 7305478), (703623.1, 7319123),
(697325.5, 7335428), (732159.0, 7341021), (722027.8, 7336020),
(729827.0, 7312423), (741539.1, 7289712), (742604.1, 7284240)
])
# Define the points for data2
data2_points = np.array([
(548055.7, 7291742), (550953.6, 7313012), (553370.8, 7302842),
(543710.0, 7306752), (542931.4, 7305144), (544490.5, 7308543),
(564527.1, 7320696), (578368.4, 7297569), (594972.5, 7299350),
(635337.6, 7300907), (649867.0, 7308755), (642136.1, 7323236),
(644287.0, 7301646), (647413.1, 7298133), (660147.9, 7294922),
(661158.7, 7304247), (677788.7, 7303612), (688676.9, 7303413),
(686350.5, 7301771), (690122.9, 7303681), (683197.7, 7308678),
(648936.3, 7312211), (637416.5, 7316094), (631286.5, 7335698),
(621747.9, 7359180), (610528.1, 7380706), (601504.8, 7402010),
(635005.1, 7320464), (680882.5, 7274095), (677057.4, 7259603)
])
# Compute convex hull for data1
hull_data1 = ConvexHull(data1_points)
hull_points_data1 = data1_points[hull_data1.vertices]
# Compute convex hull for data2
hull_data2 = ConvexHull(data2_points)
hull_points_data2 = data2_points[hull_data2.vertices]
# Plot the convex hulls
plt.figure(figsize=(10, 6))
# Plot data1 convex hull
plt.fill(hull_points_data1[:,0], hull_points_data1[:,1], 'blue', alpha=0.3, label='Data1 Convex Hull')
# Plot data2 convex hull
plt.fill(hull_points_data2[:,0], hull_points_data2[:,1], 'green', alpha=0.3, label='Data2 Convex Hull')
# Set labels and legend
plt.xlabel('Longitude')
plt.ylabel('Latitude')
plt.title('Convex Hulls of Data Points')
plt.legend()
# Show plot
plt.grid(True)
plt.show()
给予：
和
from shapely.geometry import Polygon
from shapely.geometry import MultiPoint
from scipy.spatial import ConvexHull
import numpy as np
# Define the points for data1
data1_points = np.array([
(406700.4, 7201767), (398192.5, 7206899), (402995.9, 7223166),
(424210.3, 7264295), (449045.0, 7315044), (476329.8, 7304172),
(516543.0, 7281898), (543151.1, 7292038), (553048.0, 7299822),
(551770.3, 7275758), (561142.6, 7261667), (560292.0, 7251473),
(544158.5, 7229430), (566519.0, 7224846), (578221.2, 7242260),
(589978.4, 7237221), (614511.6, 7241832), (629119.1, 7229805),
(643408.0, 7229532), (675159.5, 7261319), (679089.6, 7273392),
(688517.1, 7283675), (699653.7, 7305478), (703623.1, 7319123),
(697325.5, 7335428), (732159.0, 7341021), (722027.8, 7336020),
(729827.0, 7312423), (741539.1, 7289712), (742604.1, 7284240)
])
# Define the points for data2
data2_points = np.array([
(548055.7, 7291742), (550953.6, 7313012), (553370.8, 7302842),
(543710.0, 7306752), (542931.4, 7305144), (544490.5, 7308543),
(564527.1, 7320696), (578368.4, 7297569), (594972.5, 7299350),
(635337.6, 7300907), (649867.0, 7308755), (642136.1, 7323236),
(644287.0, 7301646), (647413.1, 7298133), (660147.9, 7294922),
(661158.7, 7304247), (677788.7, 7303612), (688676.9, 7303413),
(686350.5, 7301771), (690122.9, 7303681), (683197.7, 7308678),
(648936.3, 7312211), (637416.5, 7316094), (631286.5, 7335698),
(621747.9, 7359180), (610528.1, 7380706), (601504.8, 7402010),
(635005.1, 7320464), (680882.5, 7274095), (677057.4, 7259603)
])
# Compute convex hull for data1
hull_data1 = ConvexHull(data1_points)
hull_points_data1 = data1_points[hull_data1.vertices]
# Compute convex hull for data2
hull_data2 = ConvexHull(data2_points)
hull_points_data2 = data2_points[hull_data2.vertices]
# Convert hull points to Shapely Polygon objects
polygon_data1 = Polygon(hull_points_data1)
polygon_data2 = Polygon(hull_points_data2)
# Compute the intersection of the polygons
intersection_polygon = polygon_data1.intersection(polygon_data2)
# Check if there is an intersection
if intersection_polygon.is_empty:
intersection_area = 0.0
else:
intersection_area = intersection_polygon.area
print("Area of intersection:", intersection_area)
印刷：
Area of intersection: 7294242975.459987
现在需要根据坐标的平均位置和比例将其转换为有用的测量区域单位。
出于某种奇怪的原因，R 代码预期执行与所提供的 Python 代码相同的任务，但报告没有交叉区域，如问题评论中所报告的那样：
library(sf)
# Define the points for data1
data1_points <- matrix(c(
406700.4, 7201767, 398192.5, 7206899, 402995.9, 7223166,
424210.3, 7264295, 449045.0, 7315044, 476329.8, 7304172,
516543.0, 7281898, 543151.1, 7292038, 553048.0, 7299822,
551770.3, 7275758, 561142.6, 7261667, 560292.0, 7251473,
544158.5, 7229430, 566519.0, 7224846, 578221.2, 7242260,
589978.4, 7237221, 614511.6, 7241832, 629119.1, 7229805,
643408.0, 7229532, 675159.5, 7261319, 679089.6, 7273392,
688517.1, 7283675, 699653.7, 7305478, 703623.1, 7319123,
697325.5, 7335428, 732159.0, 7341021, 722027.8, 7336020,
729827.0, 7312423, 741539.1, 7289712, 742604.1, 7284240
), ncol = 2, byrow = TRUE)
# Define the points for data2
data2_points <- matrix(c(
548055.7, 7291742, 550953.6, 7313012, 553370.8, 7302842,
543710.0, 7306752, 542931.4, 7305144, 544490.5, 7308543,
564527.1, 7320696, 578368.4, 7297569, 594972.5, 7299350,
635337.6, 7300907, 649867.0, 7308755, 642136.1, 7323236,
644287.0, 7301646, 647413.1, 7298133, 660147.9, 7294922,
661158.7, 7304247, 677788.7, 7303612, 688676.9, 7303413,
686350.5, 7301771, 690122.9, 7303681, 683197.7, 7308678,
648936.3, 7312211, 637416.5, 7316094, 631286.5, 7335698,
621747.9, 7359180, 610528.1, 7380706, 601504.8, 7402010,
635005.1, 7320464, 680882.5, 7274095, 677057.4, 7259603
), ncol = 2, byrow = TRUE)
# Convert the points to Simple Features objects with CRS 5321
data1_sf <- st_as_sf(data.frame(data1_points), coords = c("X1", "X2"), crs = 5321)
data2_sf <- st_as_sf(data.frame(data2_points), coords = c("X1", "X2"), crs = 5321)
# Compute the convex hulls
hull_data1 <- st_convex_hull(data1_sf)
hull_data2 <- st_convex_hull(data2_sf)
# Compute the intersection of the convex hulls
intersection_polygon <- st_intersection(hull_data1, hull_data2)
# Check if intersection_polygon is empty
if (nrow(intersection_polygon) == 0) {
print("No intersection found")
} else {
# Plot the intersection polygon
plot(intersection_polygon)
}
给出：
~ $ Rscript /home/o/oOo/O@@/stackoverflow/matplotlib/intersection-3.r
Linking to GEOS 3.10.2, GDAL 3.4.1, PROJ 8.2.1; sf_use_s2() is TRUE
[1] "No intersection found"
最新问题
- 根据 FormType 中的另一个选项值定义表单选项允许的值
- 了解 U-Boot 内存占用
- 如何将受保护范围应用于Google表格中的所有用户？
- 我该怎么做才能找到 pod 被逐出的原因
- Django 管理界面：了解带有提交参数的表单保存行为
- 更改linux内核安装目录
- 无法在 Compose 中获取上下文
- GitHub 访问令牌从存储库克隆、推送和拉取的最低权限？
- Compose Popover 错误中的上下文
- 正确重置数据，这些数据在Python中会随着时间的推移而严格增加
- 不同类型的类似 Vec 的容器，可以在堆上共享同一块内存
- 如何使用 Javascript 刷新另一个选项卡
- 带有嵌套序列化器的 Django RESTFRAMEWORK 序列化器
- 如何制作一个匹配正则表达式任何开始部分的正则表达式？
- 在 VB.NET 中制作一个 DLL，我可以将其用作 DllImport
- 为什么 Spring boot 应用程序的 SSL yaml 值被忽略？
- &mut Foo 和 *mut Foo 有什么区别？
- Cron 作业和 docker
- 为什么要给StringBuilder设置容量？
- ng 构建需要永远并且永远不会完成构建
© www.soinside.com 2019 - 2024. All rights reserved.