Intersection of Two LineStrings Geopandas

AJG519

Let's say I have the following to GeoDataFrames of linestrings, one of which represents roads and one of which represents contour lines.

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

If I plot these, they look like this:

enter image description here

There are 3 contour line and 2 roads. I want to find the elevation at each point along each road. Basically I want to intersect roads and contours (which should give me 12 points) and preserve the attributes from both geodataframes (road name and elevation).

I can generate the 12 points as such by using an intersection of the unions of the two 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.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)
>>> 

However, how do I now get the road name and elevation for each of those 12 points? A spatial join does not behave as I would expect and only returns 4 points (all 12 should intersect with the line files since they were created that way by definition).

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

Any suggestions as to how I can do this?

EDIT: It appears that the issue has to do with how the intersection points are created. If I buffer the roads and contours by a very small amount, the intersection works as expected. See below:

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

The above is the desired output although I'm not sure as to why I have to buffer the lines to get them to intersect the points that were created from the intersection of the lines.

eguaio

Notice that operations unary_union and intersection are made over the geometries inside the GeoDataFrame, so you lose the data stored in the rest of the columns. I think in this case you have to do it by hand by accessing each geometry in the data frames. The following 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 

produces:

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

Notice each geometry has two points, that you can access later if you want point by point information, or you can create a row for each point introducing a for loop that iterates over the points, something like:

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

But in this case you depend on knowing that each intersection produces a multi-geometry.

About your other approach using the sjoin function, I couldn't test it because the version of geopandas I'm using does not provide the tools module. Try to put buffer(0.0) to see what happens.

この記事はインターネットから収集されたものであり、転載の際にはソースを示してください。

侵害の場合は、連絡してください[email protected]

編集
0

コメントを追加

0

関連記事

分類Dev

Intersection of two Counters

分類Dev

Intersection of two arrays in BASH

分類Dev

Checking whether intersection of two lists is valid, returning intersection itself not necessary

分類Dev

How to fill between x and two functions (intersection)?

分類Dev

Intersection of two lists maintaining duplicate values in Kotlin

分類Dev

Efficient, or fast, size of the set intersection of two vectors

分類Dev

find intersection point of two vectors independent from direction

分類Dev

How to apply the intersection between two lists in C++?

分類Dev

Finding the intersection in two lists of tuples regardless of tuple order

分類Dev

Intersection of two array which contain exact reverse order of elements

分類Dev

intersection of objects except for one attribute between two arrays

分類Dev

Trying to sort two list of numbers and using uniq to get the intersection

分類Dev

GeoPandas Label Polygons

分類Dev

Problems to import Geopandas

分類Dev

MEAN & Geospatial queries - Find LineStrings Intersecting on Another One given Its Name

分類Dev

LineStringsのPythonリストからGeoJSONを作成する方法

分類Dev

Three.js / Intersection

分類Dev

Array of objects intersection

分類Dev

Pandas intersection of groups

分類Dev

Intersection of sets as columns in pandas

分類Dev

What is the intersection of the keys of dictionaries

分類Dev

KineticJS get all intersection

分類Dev

MongoDB index intersection

分類Dev

SQL intersection - error in syntax

分類Dev

Geopandasのカラーバー

分類Dev

Subsetting NYC shapefile to Manhattan using Geopandas

分類Dev

Jupyter/Geopandas plot order breaks figsize

分類Dev

LineStringsが表示されていません。ラベルだけがします

分類Dev

WorldWind Sphere Line Intersection Bug?

Related 関連記事

  1. 1

    Intersection of two Counters

  2. 2

    Intersection of two arrays in BASH

  3. 3

    Checking whether intersection of two lists is valid, returning intersection itself not necessary

  4. 4

    How to fill between x and two functions (intersection)?

  5. 5

    Intersection of two lists maintaining duplicate values in Kotlin

  6. 6

    Efficient, or fast, size of the set intersection of two vectors

  7. 7

    find intersection point of two vectors independent from direction

  8. 8

    How to apply the intersection between two lists in C++?

  9. 9

    Finding the intersection in two lists of tuples regardless of tuple order

  10. 10

    Intersection of two array which contain exact reverse order of elements

  11. 11

    intersection of objects except for one attribute between two arrays

  12. 12

    Trying to sort two list of numbers and using uniq to get the intersection

  13. 13

    GeoPandas Label Polygons

  14. 14

    Problems to import Geopandas

  15. 15

    MEAN & Geospatial queries - Find LineStrings Intersecting on Another One given Its Name

  16. 16

    LineStringsのPythonリストからGeoJSONを作成する方法

  17. 17

    Three.js / Intersection

  18. 18

    Array of objects intersection

  19. 19

    Pandas intersection of groups

  20. 20

    Intersection of sets as columns in pandas

  21. 21

    What is the intersection of the keys of dictionaries

  22. 22

    KineticJS get all intersection

  23. 23

    MongoDB index intersection

  24. 24

    SQL intersection - error in syntax

  25. 25

    Geopandasのカラーバー

  26. 26

    Subsetting NYC shapefile to Manhattan using Geopandas

  27. 27

    Jupyter/Geopandas plot order breaks figsize

  28. 28

    LineStringsが表示されていません。ラベルだけがします

  29. 29

    WorldWind Sphere Line Intersection Bug?

ホットタグ

アーカイブ