Geopandas autodivididas Linestring geodataframe de forma rápida sem perder todos os atributos

Nov 24 2020

Tenho geodataframe de muitas LineStrings. As LineStrings se cruzam, mas não são divididas nessas interseções. Minha solução atual para conseguir isso é daqui :

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)

Porém eu perco todas as colunas com esta solução. Eu também tentei isso, mas demora muito:

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

Alguma sugestão?

Respostas

4 gene Nov 30 2020 at 11:31

Eu encontrei uma solução.

Usando meu exemplo :

a) O shapefile original

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 a geometria original para evitar problemas aritméticos flutuantes (em intersectsou within)

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

c) Use unary_unionpara dividir todas as LineStrings com intersecção automática

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 uma junção espacial (com withinou intersect) para juntar os dois dataframes e recuperar os atributos originais

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

1 GevaertJoep Nov 28 2020 at 20:14

Não é uma solução, mas pode ajudar: tentar fazer uma interseção da união fornece dois iteráveis ​​para cada atributo onde eles se cruzam. No entanto, algumas partes da linha tornam-se pontos ...:

o original:

resultado:

Talvez com algumas edições no código abaixo possa funcionar?

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