7

I have geodataframe of many LineStrings. The LineStrings intersect but are not split at those intersections. My current solution to achieve this is from here:

network = gp.read_file(filenameNetwork)
newNetwork = gp.GeoDataFrame()
for splittedGeom in network.geometry.unary_union:
    part = gp.GeoDataFrame([[splittedGeom]], columns=['geometry'])
    newNetwork = newNetwork.append(part)

However I loose all columns with this solution. I also tried this, but it takes way too long:

from shapely import ops

streets = streets.reset_index(drop=True)
streets = streets[['geometry', 'costs']]
headers = list(streets.columns)

index = 0
newStreets = gp.GeoDataFrame( columns=['geometry'])
for line in range(len(streets)-1):
    print(line, len(streets))
    linegeom = streets.at[line, 'geometry']
    isNotSplitted = True
    for line2 in range(len(streets)):
        if line2 == line:
            continue
        linegeom2 = streets.at[line2, 'geometry']
        if linegeom2.crosses(linegeom):
            try:
                linegeomsplitted = ops.split(linegeom, linegeom2)
            except:
                continue
            isNotSplitted = False
            for split in range((len(list(linegeomsplitted.geoms)))):
                splittedline = (list(linegeomsplitted.geoms))[split]
                for head in headers:
                    if head == 'geometry':
                        headValue = splittedline
                    else:
                        headValue = streets.at[line, head]
                    newStreets.at[index, head] = headValue
                index += 1
    if isNotSplitted:
        for head in headers:
            headValue = streets.at[line, head]
            newStreets.at[index, head] = headValue
        index += 1
streets = newStreets
streets = streets.drop_duplicates(subset=['geometry'],
                                  keep='first')

Any suggestions?

5
+50

I found a solution.

Using my example:

a) The original shapefile

enter image description here

import geopandas as gpd
df = gpd.read_file("stac-graphe.shp")
df
id   test                geometry
1   test1   LINESTRING (10.244 -273.317, 784.201 -222.924)
2   test2   LINESTRING (210.484 -553.461, 324.991 -4.534)
3   test3   LINESTRING (169.970 -134.276, 126.511 -218.533...
4   test4   LINESTRING (100.000 -433.317, 724.390 -112.341...
5   test5   LINESTRING (232.683 -113.317, 694.146 -445.024...
6   test6   LINESTRING (563.415 -552.341, 559.512 -22.585)

b) Buffer the original geometry to avoid float arithmetic problems (in intersects or within)

df2 = df.copy()
df2.geometry = df2.geometry.buffer(0.01)

c) Use unary_union to split all the self-intersected LineStrings

un = df.geometry.unary_union
geom = [i for i in un]
id = [j for j in range(len(geom))]
unary = gpd.GeoDataFrame({"id":id,"geometry":geom})
unary.head()
id                   geometry
0   LINESTRING (10.244 -273.317, 192.920 -261.423)
1   LINESTRING (192.920 -261.423, 272.484 -256.242)
2   LINESTRING (272.484 -256.242, 418.308 -246.748)
3   LINESTRING (418.308 -246.748, 469.403 -243.421)
4   LINESTRING (469.403 -243.421, 561.095 -237.451)

d) Use a spatial join (with within or intersect) to join the two dataframes and retrieve the original attributes

from geopandas.tools import sjoin
result =sjoin(unary, df2, how="inner",op='within')
result.head()
id_left                   geometry               index_right id_right   test
0   LINESTRING (10.244 -273.317, 192.920 -261.423)   0         1       test1
1   LINESTRING (192.920 -261.423, 272.484 -256.242)  0         1       test1
2   LINESTRING (272.484 -256.242, 418.308 -246.748)  0         1       test1
3   LINESTRING (418.308 -246.748, 469.403 -243.421)  0         1       test1
4   LINESTRING (469.403 -243.421, 561.095 -237.451)  0         1       test1

enter image description here

1

Not a solution, but it may help: Trying to make an intersection of the union gives two iterables for each attribute where they intersect. However, some line parts become points...:

the original:

enter image description here

output:

enter image description here

Maybe with some edits on below code it could work?

import geopandas as gp
network = gp.read_file(filepath)
newNetwork = gp.GeoDataFrame()

geom = network.unary_union
newNetwork = gp.GeoDataFrame(columns=network.columns)

for i in range(len(network)):
    for splittedGeom in network.intersection(network.unary_union)[i]:
        part = gp.GeoDataFrame([list(network.loc[i,network.columns[:-1]])+[splittedGeom]],columns=list(network.columns))
        newNetwork = newNetwork.append(part)

newNetwork.plot()
1
  • thanks, any idea why they become points? – Leo Nov 30 '20 at 7:08

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.