2015-11-20 38 views
11

比方说,我对线条的GeoDataFrame有以下几点,其中一条代表道路,其中一条代表轮廓线。两条LineString的交叉点Geopandas

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

如果我描绘出这些,他们是这样的:

enter image description here

有3轮廓线和2路。我想找到每条道路上每个点的高程。基本上我想交叉道路和轮廓(应该给我12分),并保留来自两个地理数据框(道路名称和高程)的属性。

我可以用两个geodataframes的工会的交集产生的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.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) 
>>> 

不过,我现在该如何获得每个这12点的路名和高度?空间连接的行为并不像我期望的那样,只返回4个点(所有12个点应该与行文件相交,因为它们是按照定义创建的)。

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

有关我如何做到这一点的任何建议?

编辑: 看来这个问题与如何创建交点有关。如果我将道路和轮廓缓冲一小部分,则交叉口按预期工作。请看下图:

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

以上是所需的输出,虽然我不知道,为什么我有缓冲线,让他们相交从线的交叉创建的点。

回答

5

请注意,操作unary_unionintersection是在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.7 2.636396103067893), ... 

注意每个几何有两点,以后可以访问,如果你想逐点信息,或者你可以为每个点引入创建一行一个for循环,在迭代点,像:

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

但在这种情况下,你需要知晓每个路口产生多几何。

关于使用sjoin函数的其他方法,我无法测试它,因为我使用的geopandas版本没有提供tools模块。尝试将buffer(0.0)看看会发生什么。