Find nearest point on each line for all linestring intersections in geopandas df

user_10997067 0
Viewed 641

I have a geopandas dataframe containing several line strings created from lat, lon point data. For all line intersections, I need to find the nearest point within each line string to that intersection.

Thus, if two lines in the dataframe intersect, I need to nearest point to that intersection in each linestring. I have used itertools to find all possible intersection points similar to the accepted answer in this post: https://gis.stackexchange.com/questions/137909/intersecting-lines-to-get-crossings-using-python-with-qgis

Is there a simpler way to find the closest point in each linestring to the intersection point, for all line intersections in a geopandas dataframe?

My dataframe looks like so, with each file name containing its own linestring:

                                                            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...
shapely geopandas pandas python
1 Answers

Let's create n random lines :

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, 0, 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)]})

This gives :

>>> 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)

Let's get our intersections :

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)

Here I just got a little misunderstanding, sometimes the inter = i.intersection(j) returns a LineString object, I have no idea how two different lines can get as an intersection output another line (unless they are the same). I leave this up to you.

And now, we can create our df with the resulting points :

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

We can see that we are more working with segments rather than pure lines, since the number of intersections (ie points) is 2220. And I do not agree to consider that we are lucky enough to have 7880 parallel lines.

Then, we import our best friend for the operation :

from shapely.ops import nearest_points

And we compute the desired output :

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

Hope, this answers your question, let me know if you got a better answer.

user_122132020
Related