Intersezione di due LineString Geopandas

Question

Diciamo che ho il seguente per GeoDataFrames di linee, una delle quali rappresenta strade e una delle quali rappresenta le linee di contorno.

>>> import geopandas as gpd 
>>> import geopandas.tools 
>>> import shapely 
>>> from shapely.geometry import * 
>>> 
>>> r1=LineString([(-1,2),(3,2.5)]) 
>>> r2=LineString([(-1,4),(3,0)]) 
>>> Roads=gpd.GeoDataFrame(['Main St','Spruce St'],geometry=[r1,r2], columns=['Name']) 
>>> Roads 
     Name     geometry 
0 Main St LINESTRING (-1 2, 3 2.5) 
1 Spruce St LINESTRING (-1 4, 3 0) 
>>> 

>>> c1=LineString(Point(1,2).buffer(.5).exterior) 
>>> c2=LineString(Point(1,2).buffer(.75).exterior) 
>>> c3=LineString(Point(1,2).buffer(.9).exterior) 
>>> Contours=gpd.GeoDataFrame([100,90,80],geometry=[c1,c2,c3], columns=['Elevation']) 
>>> Contours 
    Elevation           geometry 
0  100 LINESTRING (1.5 2, 1.497592363336099 1.9509914... 
1   90 LINESTRING (1.75 2, 1.746388545004148 1.926487... 
2   80 LINESTRING (1.9 2, 1.895666254004977 1.9117845... 
>>> 

Se io tracciare questi, appaiono come questo:

Ci sono 3 linea di contorno e 2 strade. Voglio trovare l'elevazione in ogni punto lungo ogni strada. Fondamentalmente voglio intersecare strade e contorni (che dovrebbero darmi 12 punti) e preservare gli attributi da entrambi i geodataframes (nome della strada e altitudine).

posso generare i 12 punti in quanto tali utilizzando un incrocio dei sindacati dei due geodataframes:

>>> Intersection=gpd.GeoDataFrame(geometry=list(Roads.unary_union.intersection(Contours.unary_union))) 
>>> Intersection 
             geometry 
0 POINT (0.1118644118110415 2.13898305147638) 
1 POINT (0.2674451642029509 2.158430645525369) 
2 POINT (0.72 2.636396103067893) 
3 POINT (0.4696699141100895 2.530330085889911) 
4 POINT (0.5385205980649126 2.192315074758114) 
5 POINT (0.6464466094067262 2.353553390593274) 
6 POINT (1.353553390593274 1.646446609406726) 
7 POINT (1.399321982208571 2.299915247776072) 
8  POINT (1.530330085889911 1.46966991411009) 
9 POINT (1.636396103067893 1.7) 
10 POINT (1.670759586114587 2.333844948264324) 
11 POINT (1.827239686607525 2.353404960825941) 
>>> 

Tuttavia, come faccio ora ottenere il nome della strada e l'elevazione per ciascuno di questi 12 punti? Un join spaziale non si comporta come mi aspetterei e restituisce solo 4 punti (tutti i 12 dovrebbero intersecarsi con i file di linea poiché sono stati creati in questo modo per definizione).

>>> gpd.tools.sjoin(Intersection, Roads) 
             geometry index_right  Name 
2 POINT (0.72 2.636396103067893)   1 Spruce St 
3 POINT (0.4696699141100895 2.530330085889911)   1 Spruce St 
5 POINT (0.6464466094067262 2.353553390593274)   1 Spruce St 
6 POINT (1.353553390593274 1.646446609406726)   1 Spruce St 
>>> 

Qualche suggerimento su come posso farlo?

MODIFICA: Sembra che il problema abbia a che fare con la creazione dei punti di intersezione. Se bufferizzo le strade e i contorni di una quantità molto piccola, l'intersezione funziona come previsto. Vedere di seguito:

>>> RoadsBuff=gpd.GeoDataFrame(Roads, geometry=Roads.buffer(.000005)) 
>>> ContoursBuff=gpd.GeoDataFrame(Contours, geometry=Contours.buffer(.000005)) 
>>> 
>>> Join1=gpd.tools.sjoin(Intersection, RoadsBuff).drop('index_right',1).sort_index() 
>>> Join2=gpd.tools.sjoin(Join1, ContoursBuff).drop('index_right',1).sort_index() 
>>> 
>>> Join2 
              geometry  Name Elevation 
0 POLYGON ((1.636395933642091 1.363596995290097,... Spruce St   80 
1 POLYGON ((1.530329916464109 1.469663012468079,... Spruce St   90 
2 POLYGON ((1.353553221167472 1.646439707764716,... Spruce St  100 
3 POLYGON ((0.5385239436706243 2.192310454047735... Main St  100 
4 POLYGON ((0.2674491823047923 2.158426108877007... Main St   90 
5 POLYGON ((0.1118688004427904 2.138978561144256... Main St   80 
6 POLYGON ((0.6464467873602107 2.353546141571978... Spruce St  100 
7 POLYGON ((0.4696700920635739 2.530322836868614... Spruce St   90 
8 POLYGON ((0.3636040748855915 2.636388854046597... Spruce St   80 
9 POLYGON ((1.399312865255344 2.299919147068011,... Main St  100 
10 POLYGON ((1.670752113626148 2.333849053114361,... Main St   90 
11 POLYGON ((1.827232214119086 2.353409065675979,... Main St   80 
>>> 

Quanto sopra è il risultato desiderato, anche se non sono sicuro del perché devo tamponare le linee per farli intersecare i punti che sono stati creati dalla intersezione delle linee.

Answer 1

Si noti che le operazioni unary_union e vengono eseguite sulle geometrie all'interno dello GeoDataFrame, quindi si perdono i dati memorizzati nel resto delle colonne. Penso che in questo caso devi farlo a mano accedendo a ciascuna geometria nei frame di dati. Il seguente codice:

import geopandas as gpd 
from shapely.geometry import LineString, Point 

r1=LineString([(-1,2),(3,2.5)]) 
r2=LineString([(-1,4),(3,0)]) 
roads=gpd.GeoDataFrame(['Main St','Spruce St'],geometry=[r1,r2], columns=['Name']) 

c1=LineString(Point(1,2).buffer(.5).exterior) 
c2=LineString(Point(1,2).buffer(.75).exterior) 
c3=LineString(Point(1,2).buffer(.9).exterior) 
contours=gpd.GeoDataFrame([100,90,80],geometry=[c1,c2,c3], columns=['Elevation']) 

columns_data = [] 
geoms = [] 
for _, n, r in roads.itertuples(): 
    for _, el, c in contours.itertuples(): 
     intersect = r.intersection(c) 
     columns_data.append((n,el)) 
     geoms.append(intersect) 

all_intersection = gpd.GeoDataFrame(columns_data, geometry=geoms, 
        columns=['Name', 'Elevation']) 

print all_intersection 

produce:

 Name Elevation           geometry 
0 Main St  100 (POINT (0.5385205980649125 2.192315074758114),... 
1 Main St   90 (POINT (0.2674451642029509 2.158430645525369),... 
2 Main St   80 (POINT (0.1118644118110415 2.13898305147638), ... 
3 Spruce St  100 (POINT (0.6464466094067262 2.353553390593274),... 
4 Spruce St   90 (POINT (0.4696699141100893 2.53033008588991), ... 
5 Spruce St   80 (POINT (0.7 2.636396103067893), ... 

Avviso ogni geometria ha due punti, che è possibile accedere in seguito, se si desidera punto per informazioni sul punto, oppure è possibile creare una riga per ogni punto di introduzione un ciclo for che scorre oltre i punti, qualcosa come:

for p in intersect: 
    columns_data.append((n,el)) 
    geoms.append(p) 

Ma in questo caso si può dipendere dal fatto che ogni intersezione produce un multi-geometria.

Informazioni sull'altro approccio che utilizza la funzione sjoin, non è stato possibile testarlo poiché la versione di geopandas che sto utilizzando non fornisce il modulo tools. Prova a inserire buffer(0.0) per vedere cosa succede.