计算两点之间跨多边形的相交距离
Calculate intersected distance across polygons between two points
我的形状文件包括不同多边形的所有坐标（纬度、经度），例如此处显示的黑色、红色和绿色。我想计算 AB 线跨越不同多边形的距离。我正在使用 geopandas 并且也在探索 gdal，但还没有想出任何函数可以允许这种计算或者我需要从头开始实现。谢谢！
剪辑和空间连接怎么样？这是一个玩具示例。 gdf_clip_sj 中的 index_right 列包含来自 poly_gdf.
的相应多边形的索引
import geopandas as gpd
import pandas as pd
from shapely.geometry import Point, Polygon, LineString
ls_geom = [LineString([Point(0, 0), Point(5, 5)]),
LineString([Point(5, 5), Point(20, 20)])]
ls_gdf = gpd.GeoDataFrame(geometry=ls_geom)
poly_geom = [Polygon([Point(1, 2), Point(1, 4), Point(3, 2)]),
Polygon([Point(5, 10), Point(5, 20), Point(15, 10)])]
poly_gdf = gpd.GeoDataFrame(geometry=poly_geom)
gdf_clip = ls_gdf.clip(poly_gdf)
gdf_clip_sj = gpd.sjoin(gdf_clip, poly_gdf)
gdf_clip_sj['length'] = gdf_clip_sj['geometry'].length
gdf_clip_sj 的输出：
geometry index_right length
0 LINESTRING (2.00000 2.00000, 2.50000 2.50000) 0 0.707107
1 LINESTRING (10.00000 10.00000, 12.50000 12.50000) 1 3.535534
即：
gdf_cat = pd.concat([poly_gdf, ls_gdf]).plot(facecolor='none')
gdf_clip_cat = pd.concat([gdf_clip_sj, poly_gdf])
gdf_clip_cat.plot(facecolor='none')
编辑：值得注意的是，我最初的答案是使用笛卡尔数学，但由于您引用了 lat/long crs，您可能需要使用 Rob 的优秀答案中所述的半正弦距离。如果你不想安装任何额外的包，你可以使用 derricw in this interesting post.
给出的公式
import numpy as np
def haversine_np(lon1, lat1, lon2, lat2):
"""
Calculate the great circle distance between two points
on the earth (specified in decimal degrees)
All args must be of equal length.
"""
lon1, lat1, lon2, lat2 = map(np.radians, [lon1, lat1, lon2, lat2])
dlon = lon2 - lon1
dlat = lat2 - lat1
a = np.sin(dlat/2.0)**2 + np.cos(lat1) * np.cos(lat2) * np.sin(dlon/2.0)**2
c = 2 * np.arcsin(np.sqrt(a))
km = 6367 * c
return km
我们可以通过将LineString坐标提取到coords列，然后通过havesine_np函数apply来实现。 FWIW，此示例中的索引仅在每行只有 2 个点时才有效。
gdf_clip_sj['coords'] = gdf_clip_sj.apply(lambda x: [y for y in x['geometry'].coords], axis=1)
gdf_clip_sj['hd'] = gdf_clip_sj['coords'].apply(lambda x: haversine_np(x[0][0], x[0][1], x[1][0], x[1][1]))
gdf_clip_sj 的输出：
geometry index_right length coords hd
0 LINESTRING (2.00000 2.00000, 2.50000 2.50000) 0 0.707107 [(2.0, 2.0), (2.5, 2.5)] 78.546912
1 LINESTRING (10.00000 10.00000, 12.50000 12.50000) 1 3.535534 [(10.0, 10.0), (12.5, 12.5)] 389.112802
在 QGIS 中快速抽查：
- 已经生成了一些几何图形。乌克兰上空的网格和对角线
- 使用https://geopandas.org/en/stable/docs/reference/api/geopandas.GeoSeries.intersection.html寻找交点
- 用https://geopy.readthedocs.io/en/stable/#module-geopy.distance计算直线起点到交点的距离
import shapely
import geopandas as gpd
import geopy.distance
world = gpd.read_file(gpd.datasets.get_path("naturalearth_lowres"))
dims = (4, 3)
a, b, c, d = world.loc[world["name"].eq("Ukraine")].total_bounds
g = 0.4
# generate some polygons and a line that cuts through some of the polygons
gdf_grid = gpd.GeoDataFrame(
geometry=[
shapely.geometry.box(minx + g, miny + g, maxx - g, maxy - g)
for minx, maxx in zip(
np.linspace(a, c, dims[0]), np.linspace(a, c, dims[0])[1:]
)
for miny, maxy in zip(
np.linspace(b, d, dims[1]), np.linspace(b, d, dims[1])[1:]
)
],
crs="epsg:4326",
)
gdf_line = gpd.GeoDataFrame(
geometry=[
shapely.geometry.LineString(
[shapely.geometry.Point(a, b), shapely.geometry.Point(c, d)]
)
],
crs="epsg:4326",
)
# start point of linestring, point we sill calc distance from
l_start = gdf_line.loc[0, "geometry"].coords[0]
# associate every polygon with the line being measures
df_x = gdf_grid.merge(gdf_line.rename(columns={"geometry": "line"}), how="cross")
# intersection will give a line string of the points where line intersects polygon
# from each of these points calc distance to start point
gdf_grid["distances"] = df_x.intersection(gpd.GeoSeries(df_x["line"])).apply(
lambda g: [round(geopy.distance.geodesic(l_start, p).km, 2) for p in g.coords]
)
# visualize ....
m = gdf_grid.explore(height=200, width=400)
gdf_line.explore(m=m, color="red")
数据帧
这条线没有穿过一些多边形，因此没有距离。
distances
geometry
0
[120.44, 658.42]
POLYGON ((27.68400190604637 45.69330801043112, 27.68400190604637 48.41419129088147, 22.48560835133485 48.41419129088147, 22.48560835133485 45.69330801043112, 27.68400190604637 45.69330801043112))
1
[]
POLYGON ((27.68400190604637 49.21419129088147, 27.68400190604637 51.9350745713318, 22.48560835133485 51.9350745713318, 22.48560835133485 49.21419129088147, 27.68400190604637 49.21419129088147))
2
[752.25, 937.0]
POLYGON ((33.68239546075789 45.69330801043112, 33.68239546075789 48.41419129088147, 28.48400190604637 48.41419129088147, 28.48400190604637 45.69330801043112, 33.68239546075789 45.69330801043112))
3
[1176.08, 1360.15]
POLYGON ((33.68239546075789 49.21419129088147, 33.68239546075789 51.9350745713318, 28.48400190604637 51.9350745713318, 28.48400190604637 49.21419129088147, 33.68239546075789 49.21419129088147))
4
[]
POLYGON ((39.68078901546941 45.69330801043112, 39.68078901546941 48.41419129088147, 34.48239546075789 48.41419129088147, 34.48239546075789 45.69330801043112, 39.68078901546941 45.69330801043112))
5
[1453.4, 1985.1]
POLYGON ((39.68078901546941 49.21419129088147, 39.68078901546941 51.9350745713318, 34.48239546075789 51.9350745713318, 34.48239546075789 49.21419129088147, 39.68078901546941 49.21419129088147))
geometry
0
LINESTRING (22.08560835133486 45.29330801043113, 40.08078901546941 52.3350745713318)
可视化
与悬停的多边形相交的两点线相关的两个距离
多条线，不规则多边形
- 问题范围已扩大：如果多边形形状不规则且共享边界怎么办？如何检查哪条相交线属于哪个多边形？
- 国家边界不规则且共享。将国家/地区用作多边形
- 已创建 3 条线（水平线、垂直线和对角线）。使用完全相同的方法，遍历线并将距离作为列添加到多边形地理数据框
gdf_poly = world.loc[(world["continent"]=="Europe") & (~world["iso_a3"].isin(["-99","RUS"]))].copy().reset_index(drop=True)
a,b,c,d = gdf_poly.total_bounds
gdf_line = gpd.GeoDataFrame(
data={"line": ["diagonal", "vertical", "horizontal"]},
geometry=[
shapely.geometry.LineString([(a, b), (c, d)]),
shapely.geometry.LineString([((a + c) / 2, b), ((a + c) / 2, d)]),
shapely.geometry.LineString([(a, (b + d) / 2), (c, (b + d) / 2)]),
],
crs="epsg:4326",
)
# multiple lines, put distances into columns
for _, r in gdf_line.iterrows():
l_start = r["geometry"].coords[0]
gdf_poly[r["line"]] = gdf_poly.intersection(
gpd.GeoSeries([r["geometry"]] * len(gdf_poly), crs=gdf_line.crs)
).apply(
lambda g: [
round(geopy.distance.geodesic(l_start, p).km, 2)
for gg in ([g] if isinstance(g, shapely.geometry.LineString) else g.geoms)
for p in gg.coords
]
)
m = gdf_poly.explore(height=400, width=800)
gdf_line.explore(m=m, color="red")
- 德国与垂直线相交四次，因此它们是垂直线的四个距离
我的形状文件包括不同多边形的所有坐标（纬度、经度），例如此处显示的黑色、红色和绿色。我想计算 AB 线跨越不同多边形的距离。我正在使用 geopandas 并且也在探索 gdal，但还没有想出任何函数可以允许这种计算或者我需要从头开始实现。谢谢！
剪辑和空间连接怎么样？这是一个玩具示例。 gdf_clip_sj 中的 index_right 列包含来自 poly_gdf.
的相应多边形的索引import geopandas as gpd
import pandas as pd
from shapely.geometry import Point, Polygon, LineString
ls_geom = [LineString([Point(0, 0), Point(5, 5)]),
LineString([Point(5, 5), Point(20, 20)])]
ls_gdf = gpd.GeoDataFrame(geometry=ls_geom)
poly_geom = [Polygon([Point(1, 2), Point(1, 4), Point(3, 2)]),
Polygon([Point(5, 10), Point(5, 20), Point(15, 10)])]
poly_gdf = gpd.GeoDataFrame(geometry=poly_geom)
gdf_clip = ls_gdf.clip(poly_gdf)
gdf_clip_sj = gpd.sjoin(gdf_clip, poly_gdf)
gdf_clip_sj['length'] = gdf_clip_sj['geometry'].length
gdf_clip_sj 的输出：
geometry index_right length
0 LINESTRING (2.00000 2.00000, 2.50000 2.50000) 0 0.707107
1 LINESTRING (10.00000 10.00000, 12.50000 12.50000) 1 3.535534
即：
gdf_cat = pd.concat([poly_gdf, ls_gdf]).plot(facecolor='none')
gdf_clip_cat = pd.concat([gdf_clip_sj, poly_gdf])
gdf_clip_cat.plot(facecolor='none')
编辑：值得注意的是，我最初的答案是使用笛卡尔数学，但由于您引用了 lat/long crs，您可能需要使用 Rob 的优秀答案中所述的半正弦距离。如果你不想安装任何额外的包，你可以使用 derricw in this interesting post.
给出的公式import numpy as np
def haversine_np(lon1, lat1, lon2, lat2):
"""
Calculate the great circle distance between two points
on the earth (specified in decimal degrees)
All args must be of equal length.
"""
lon1, lat1, lon2, lat2 = map(np.radians, [lon1, lat1, lon2, lat2])
dlon = lon2 - lon1
dlat = lat2 - lat1
a = np.sin(dlat/2.0)**2 + np.cos(lat1) * np.cos(lat2) * np.sin(dlon/2.0)**2
c = 2 * np.arcsin(np.sqrt(a))
km = 6367 * c
return km
我们可以通过将LineString坐标提取到coords列，然后通过havesine_np函数apply来实现。 FWIW，此示例中的索引仅在每行只有 2 个点时才有效。
gdf_clip_sj['coords'] = gdf_clip_sj.apply(lambda x: [y for y in x['geometry'].coords], axis=1)
gdf_clip_sj['hd'] = gdf_clip_sj['coords'].apply(lambda x: haversine_np(x[0][0], x[0][1], x[1][0], x[1][1]))
gdf_clip_sj 的输出：
geometry index_right length coords hd
0 LINESTRING (2.00000 2.00000, 2.50000 2.50000) 0 0.707107 [(2.0, 2.0), (2.5, 2.5)] 78.546912
1 LINESTRING (10.00000 10.00000, 12.50000 12.50000) 1 3.535534 [(10.0, 10.0), (12.5, 12.5)] 389.112802
在 QGIS 中快速抽查：
- 已经生成了一些几何图形。乌克兰上空的网格和对角线
- 使用https://geopandas.org/en/stable/docs/reference/api/geopandas.GeoSeries.intersection.html寻找交点
- 用https://geopy.readthedocs.io/en/stable/#module-geopy.distance计算直线起点到交点的距离
import shapely
import geopandas as gpd
import geopy.distance
world = gpd.read_file(gpd.datasets.get_path("naturalearth_lowres"))
dims = (4, 3)
a, b, c, d = world.loc[world["name"].eq("Ukraine")].total_bounds
g = 0.4
# generate some polygons and a line that cuts through some of the polygons
gdf_grid = gpd.GeoDataFrame(
geometry=[
shapely.geometry.box(minx + g, miny + g, maxx - g, maxy - g)
for minx, maxx in zip(
np.linspace(a, c, dims[0]), np.linspace(a, c, dims[0])[1:]
)
for miny, maxy in zip(
np.linspace(b, d, dims[1]), np.linspace(b, d, dims[1])[1:]
)
],
crs="epsg:4326",
)
gdf_line = gpd.GeoDataFrame(
geometry=[
shapely.geometry.LineString(
[shapely.geometry.Point(a, b), shapely.geometry.Point(c, d)]
)
],
crs="epsg:4326",
)
# start point of linestring, point we sill calc distance from
l_start = gdf_line.loc[0, "geometry"].coords[0]
# associate every polygon with the line being measures
df_x = gdf_grid.merge(gdf_line.rename(columns={"geometry": "line"}), how="cross")
# intersection will give a line string of the points where line intersects polygon
# from each of these points calc distance to start point
gdf_grid["distances"] = df_x.intersection(gpd.GeoSeries(df_x["line"])).apply(
lambda g: [round(geopy.distance.geodesic(l_start, p).km, 2) for p in g.coords]
)
# visualize ....
m = gdf_grid.explore(height=200, width=400)
gdf_line.explore(m=m, color="red")
数据帧
这条线没有穿过一些多边形，因此没有距离。
| distances | geometry | |
|---|---|---|
| 0 | [120.44, 658.42] | POLYGON ((27.68400190604637 45.69330801043112, 27.68400190604637 48.41419129088147, 22.48560835133485 48.41419129088147, 22.48560835133485 45.69330801043112, 27.68400190604637 45.69330801043112)) |
| 1 | [] | POLYGON ((27.68400190604637 49.21419129088147, 27.68400190604637 51.9350745713318, 22.48560835133485 51.9350745713318, 22.48560835133485 49.21419129088147, 27.68400190604637 49.21419129088147)) |
| 2 | [752.25, 937.0] | POLYGON ((33.68239546075789 45.69330801043112, 33.68239546075789 48.41419129088147, 28.48400190604637 48.41419129088147, 28.48400190604637 45.69330801043112, 33.68239546075789 45.69330801043112)) |
| 3 | [1176.08, 1360.15] | POLYGON ((33.68239546075789 49.21419129088147, 33.68239546075789 51.9350745713318, 28.48400190604637 51.9350745713318, 28.48400190604637 49.21419129088147, 33.68239546075789 49.21419129088147)) |
| 4 | [] | POLYGON ((39.68078901546941 45.69330801043112, 39.68078901546941 48.41419129088147, 34.48239546075789 48.41419129088147, 34.48239546075789 45.69330801043112, 39.68078901546941 45.69330801043112)) |
| 5 | [1453.4, 1985.1] | POLYGON ((39.68078901546941 49.21419129088147, 39.68078901546941 51.9350745713318, 34.48239546075789 51.9350745713318, 34.48239546075789 49.21419129088147, 39.68078901546941 49.21419129088147)) |
| geometry | |
|---|---|
| 0 | LINESTRING (22.08560835133486 45.29330801043113, 40.08078901546941 52.3350745713318) |
可视化
与悬停的多边形相交的两点线相关的两个距离
多条线，不规则多边形
- 问题范围已扩大：如果多边形形状不规则且共享边界怎么办？如何检查哪条相交线属于哪个多边形？
- 国家边界不规则且共享。将国家/地区用作多边形
- 已创建 3 条线（水平线、垂直线和对角线）。使用完全相同的方法，遍历线并将距离作为列添加到多边形地理数据框
gdf_poly = world.loc[(world["continent"]=="Europe") & (~world["iso_a3"].isin(["-99","RUS"]))].copy().reset_index(drop=True)
a,b,c,d = gdf_poly.total_bounds
gdf_line = gpd.GeoDataFrame(
data={"line": ["diagonal", "vertical", "horizontal"]},
geometry=[
shapely.geometry.LineString([(a, b), (c, d)]),
shapely.geometry.LineString([((a + c) / 2, b), ((a + c) / 2, d)]),
shapely.geometry.LineString([(a, (b + d) / 2), (c, (b + d) / 2)]),
],
crs="epsg:4326",
)
# multiple lines, put distances into columns
for _, r in gdf_line.iterrows():
l_start = r["geometry"].coords[0]
gdf_poly[r["line"]] = gdf_poly.intersection(
gpd.GeoSeries([r["geometry"]] * len(gdf_poly), crs=gdf_line.crs)
).apply(
lambda g: [
round(geopy.distance.geodesic(l_start, p).km, 2)
for gg in ([g] if isinstance(g, shapely.geometry.LineString) else g.geoms)
for p in gg.coords
]
)
m = gdf_poly.explore(height=400, width=800)
gdf_line.explore(m=m, color="red")
- 德国与垂直线相交四次，因此它们是垂直线的四个距离