Пересечение двух геопанд LineStrings

Скажем, у меня есть следующие строки линий GeoDataFrames, одна из которых представляет дороги, а другая - горизонтальные линии.

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

Если я построю их, они будут выглядеть так:

введите описание изображения здесь

Есть 3 контурных линии и 2 дороги. Я хочу найти высоту в каждой точке на каждой дороге. В основном я хочу пересечь дороги и контуры (что должно дать мне 12 баллов) и сохранить атрибуты из обоих фреймов геоданных (название дороги и отметка).

Я могу сгенерировать 12 точек как таковые, используя пересечение объединений двух кадров геоданных:

>>> 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.3636038969321072 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.363603896932107)
10   POINT (1.670759586114587 2.333844948264324)
11   POINT (1.827239686607525 2.353404960825941)
>>>

Однако как мне теперь получить название дороги и высоту для каждой из этих 12 точек? Пространственное соединение ведет себя не так, как я ожидал, и возвращает только 4 точки (все 12 должны пересекаться с файлами линий, поскольку они были созданы таким образом по определению).

>>> gpd.tools.sjoin(Intersection, Roads)
                                       geometry  index_right       Name
2  POINT (0.3636038969321072 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
>>>

Есть предложения, как я могу это сделать?

РЕДАКТИРОВАТЬ: Похоже, проблема связана с тем, как создаются точки пересечения. Если я немного буферизирую дороги и контуры, перекресток будет работать должным образом. См. ниже:

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

Вышеуказанное является желаемым результатом, хотя я не уверен, почему мне нужно буферизовать линии, чтобы заставить их пересекать точки, которые были созданы из пересечения линий.


python geopandas shapely
person AJG519    schedule 20.11.2015    source источник

Ответы (1)


arrow_upward
8
arrow_downward

Обратите внимание, что операции unary_union и intersection выполняются над геометриями внутри GeoDataFrame, поэтому вы теряете данные, хранящиеся в остальных столбцах. Я думаю, что в этом случае вам придется делать это вручную, обращаясь к каждой геометрии во фреймах данных. Следующий код:

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

производит:

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

Обратите внимание, что каждая геометрия имеет две точки, к которым вы можете получить доступ позже, если вам нужна информация по точкам, или вы можете создать строку для каждой точки, вводя цикл for, который повторяется по точкам, что-то вроде:

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

Но в этом случае вам нужно знать, что каждое пересечение создает мульти-геометрию.

Что касается вашего другого подхода с использованием функции sjoin, я не смог его протестировать, потому что версия geopandas, которую я использую, не предоставляет модуль tools. Попробуй поставить buffer(0.0), что получится.

person eguaio    schedule 14.12.2015