在geopandas df中找到所有线串交点的每条线的最近点

编程问答 2022-06-09

问题描述

我有一个geopandas数据框，其中包含从纬度，经点数据创建的几个线串。对于所有线相交，我需要找到每个线串中距该相交的最近点。

因此，如果数据框中的两条线相交，则需要最接近每个线串中该点的交点。我已经使用itertools查找了所有可能的相交点，类似于本文中的已接受答案：https://gis.stackexchange.com/questions/137909/intersecting-lines-to-get-crossings-using-python-with-qgis

对于geopandas数据框中的所有线相交点，是否有一种更简单的方法来找到每个线串中最接近交点的点？

我的数据框看起来像这样，每个文件名都包含自己的线串：

                                                            geometry
file                                                            
2015_may14_10  LInesTRING (-140.43855 59.80302,-140.44101 59...
2015_may14_11  LInesTRING (-140.84909 59.83433,-140.84758 59...
2015_may14_12  LInesTRING (-140.66859 59.79890,-140.66600 59...
2015_may14_15  LInesTRING (-140.19642 59.86655,-140.19795 59...
2015_may14_16  LInesTRING (-141.08783 59.94741,-141.08610 59...

解决方法

让我们创建 n 条随机线：

import geopandas as gpd
from shapely.geometry import LineString,Point,Polygon
from shapely import wkt
import numpy as np
xmin,xmax,ymin,ymax = 0,10000,10000
n = 100
xa = (xmax - xmin) * np.random.random(n) + xmin
ya = (ymax - ymin) * np.random.random(n) + ymin
xb = (xmax - xmin) * np.random.random(n) + xmin
yb = (ymax - ymin) * np.random.random(n) + ymin

lines = gpd.GeoDataFrame({'index':range(n),'geometry':[LineString([(a,b),(k,l)]) for a,b,k,l in zip(xa,ya,xb,yb)]})

这给出：

>>> lines
index   geometry
0   0   LINESTRING (4444.630 3081.439,6132.674 5849.463)
1   1   LINESTRING (7015.940 6378.245,4568.386 757.205)
2   2   LINESTRING (8766.417 6070.131,690.359 7511.385)
3   3   LINESTRING (4245.544 4009.196,8496.307 1557.175)
4   4   LINESTRING (1489.436 9364.784,2109.740 5923.480)
...     ...     ...
95  95  LINESTRING (4783.454 7840.857,1935.396 2435.260)
96  96  LINESTRING (1884.455 4982.662,6257.958 3580.912)
97  97  LINESTRING (7072.811 7843.319,4811.589 2486.040)
98  98  LINESTRING (6933.272 6427.046,7528.579 2064.067)
99  99  LINESTRING (3876.400 5183.790,5360.753 1901.207)

让我们得到我们的交点：

res = []
for i in lines.loc[:,'geometry']:
    for j in lines.loc[:,'geometry']:
        inter = i.intersection(j)
        if inter.geom_type != 'LineString':
            res.append(inter)

这里我有点误解，有时 inter = i.intersection(j) 返回一个 LineString 对象，我不知道两条不同的线如何作为交点输出另一条线（除非它们相同） .我把这留给你。

现在，我们可以使用结果点创建我们的 df ：

points = gpd.GeoDataFrame({'geometry':res})
>>>points

    geometry
0   POINT (4811.366 3682.806)
1   POINT (5149.727 4237.644)
2   POINT (4607.312 3348.202)
3   POINT (6026.639 5675.588)
4   POINT (4514.359 3195.779)
...     ...
2215    POINT (4788.793 3166.070)
2216    POINT (4704.895 3351.608)
2217    POINT (4581.390 3624.734)
2218    POINT (4320.392 4201.921)
2219    POINT (4949.041 2811.691)

2220 rows × 1 columns

我们可以看到我们更多地使用线段而不是纯线，因为交叉点（即点）的数量是 2220。我不同意认为我们很幸运有7880平行线。

然后，我们为操作导入我们最好的朋友：

from shapely.ops import nearest_points

然后我们计算所需的输出：

intersection = []
line = []
my_point = []

for i in points.index:
    for j in lines.index:
        intersection.append(points.loc[i,'geometry'])
        line.append(lines.loc[j,'geometry'])
        my_point.append([p.wkt for p in nearest_points(points.loc[i,'geometry'],lines.loc[j,'geometry'])][1])


result = gpd.GeoDataFrame({'intersection':intersection,'line':line,'nearest_point':my_point})

result.geometry = result.loc[:,'nearest_point'].apply(wkt.loads)
result.drop(columns=['nearest_point'],inplace=True)

>>>result

intersection    line    geometry
0   POINT (4811.365980053641 3682.805619834874)     LINESTRING (4444.630325108094 3081.43918610815...   POINT (4811.366 3682.806)
1   POINT (4811.365980053641 3682.805619834874)     LINESTRING (7015.939846319573 6378.24453843603...   POINT (5677.967 3305.464)
2   POINT (4811.365980053641 3682.805619834874)     LINESTRING (8766.416847858662 6070.13073873083...   POINT (5346.331 6680.480)
3   POINT (4811.365980053641 3682.805619834874)     LINESTRING (4245.544341245415 4009.19558793877...   POINT (4811.366 3682.806)
4   POINT (4811.365980053641 3682.805619834874)     LINESTRING (1489.4355376526 9364.784164867619,...   POINT (2109.740 5923.480)
...     ...     ...     ...
221995  POINT (4949.040525093341 2811.690701237854)     LINESTRING (4783.453909575222 7840.85687296287...   POINT (2745.435 3972.709)
221996  POINT (4949.040525093341 2811.690701237854)     LINESTRING (1884.454611847149 4982.66168904636...   POINT (5294.551 3889.693)
221997  POINT (4949.040525093341 2811.690701237854)     LINESTRING (7072.811488307434 7843.31900543939...   POINT (4949.041 2811.691)
221998  POINT (4949.040525093341 2811.690701237854)     LINESTRING (6933.272054846982 6427.04550331467...   POINT (7381.288 3143.559)
221999  POINT (4949.040525093341 2811.690701237854)     LINESTRING (3876.399925481877 5183.78974899146...   POINT (4949.041 2811.691)

222000 rows × 3 columns

希望这能回答您的问题，如果您有更好的答案，请告诉我。

geopandaspandaspandaspythonshapely