Kreuzung von zwei Linienstrings Geopandas

Question

Angenommen, ich habe folgende GeoDataFrames von Liniensträngen, von denen einer Straßen darstellt und von denen einer Konturlinien darstellt.

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

Wenn ich diese Handlung, sie sehen wie folgt aus:

Es gibt 3 Konturlinie und zwei Straßen. Ich möchte die Höhe an jedem Punkt entlang jeder Straße finden. Grundsätzlich möchte ich Straßen und Konturen schneiden (was mir 12 Punkte geben soll) und die Attribute beider Geodatenrahmen (Straßenname und -höhe) beibehalten.

kann ich die 12 Punkte als solche erzeugen, indem sie einen Schnittpunkt der Gewerkschaften der beiden geodataframes mit:

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

aber wie bekomme ich jetzt den Straßennamen und Höhe für jede dieser 12 Punkte? Ein räumlicher Join verhält sich nicht wie erwartet und gibt nur 4 Punkte zurück (alle 12 sollten sich mit den Liniendateien schneiden, da sie per Definition auf diese Weise erstellt wurden).

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

Irgendwelche Vorschläge, wie ich das tun kann?

EDIT: Es scheint, dass das Problem damit zu tun hat, wie die Schnittpunkte erstellt werden. Wenn ich die Straßen und Konturen um einen sehr kleinen Betrag puffe, funktioniert die Kreuzung wie erwartet. Siehe unten:

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

Die oben ist die gewünschte Ausgabe obwohl ich nicht sicher bin, warum ich die Zeilen in den Puffer haben, um sie um die Punkte zu schneiden, die aus dem Schnittpunkt der Linien erstellt wurden.

Answer 1

Beachten Sie, dass die Operationen unary_union und intersection über die Geometrien innerhalb der GeoDataFrame vorgenommen werden, so dass Sie die in den restlichen Spalten gespeicherten Daten verlieren. Ich denke in diesem Fall müssen Sie es von Hand tun, indem Sie auf jede Geometrie in den Datenrahmen zugreifen. Der folgende Code:

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 

produziert:

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

Hinweis jede Geometrie zwei Punkte hat, dass Sie später darauf zugreifen können, wenn Sie Punkt für Punkt Informationen wünschen, oder Sie können eine Zeile für jeden Punkt erstellen Einführung ein for-Schleife, dass iteriert über die Punkte, so etwas wie:

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

Aber in diesem Fall, dass Sie zu wissen, abhängen, dass jede Kreuzung eine Multi-Geometrie erzeugt.

Über Ihre anderen Ansatz mit der sjoin Funktion konnte ich es nicht testen, weil die Version von geopandas Ich verwende nicht das Modul tools. Versuchen Sie buffer(0.0) zu sehen, was passiert.