[GIS] Geopandas intersects/contains fails between GeoDataFrames

geodataframegeopandaspythonspatial relationship

I am struggling to make the GeoDataFrame functions: intersects/contains work between GeoDataFrames instances.

Let t be a GeoDataFrame with contents:

    id                                          geom
0  230  POINT (154860.3254511474 -468176.1141802836)
1  231  POINT (149247.4170036312 -436634.9002944119)
2  232  POINT (192530.4468407042 -465274.3331522844)
3  233  POINT (138125.4083732415 -419237.1333624315)
4  234  POINT (195314.0952219047 -468429.8881146177)
5  235  POINT (196320.0493430012 -479448.5877129403)
6  236  POINT (199780.1489409541 -446899.1169124786)
7  237  POINT (161437.1203666388 -440485.6615666826)
8  238  POINT (192413.1253873391 -465423.3888242305)
9  239  POINT (196039.9555616449 -473696.9492107867)

I extract the 5 last points in another GeoDataFrame named prt:

    id                                          geom
5  235  POINT (196320.0493430012 -479448.5877129403)
6  236  POINT (199780.1489409541 -446899.1169124786)
7  237  POINT (161437.1203666388 -440485.6615666826)
8  238  POINT (192413.1253873391 -465423.3888242305)
9  239  POINT (196039.9555616449 -473696.9492107867)

executing the command:

t.intersects(prt)
Out:
0    False
1    False
2    False
3    False
4    False
5     True
6     True
7     True
8     True
9     True

However when I create a polygon from the points of prt and repeat the process it fails:

import shapely.geometry as sh
polygon = sh.polygon.Polygon([[p.x, p.y] for p in prt.geometry.tolist()])

import geopandas as gpd
pds = gpd.GeoDataFrame([polygon], columns=['geom'], geometry = 'geom', crs = {'init':'epsg:3309'})


pds.geom_equals(polygon)
Out[873]: 
0    True
dtype: bool

t.intersects(polygon)
Out[875]: 
0    False
1    False
2     True
3    False
4     True
5     True
6     True
7     True
8     True
9     True
dtype: bool

t.intersects(pds)
Out[874]: 
0    False
1    False
2    False
3    False
4    False
5    False
6    False
7    False
8    False
9    False
dtype: bool

Can you please provide me some help/guidelines? I am pretty new in Python, so any advice is welcome.

Best Answer

With the solution of How to find which points intersect with a polygon in geopandas? and your data

....
polygon = sh.polygon.Polygon([[p.x, p.y] for p in prt.geometry.tolist()])
polygon.wkt
'POLYGON ((196320.0493430012 -479448.5877129403, 199780.1489409541 -446899.1169124786, 161437.1203666388 -440485.6615666826, 192413.1253873391 -465423.3888242305, 196039.9555616449 -473696.9492107867, 196320.0493430012 -479448.5877129403))'
t.intersects(polygon)
0    False
1    False
2     True
3    False
4     True
5     True
6     True
7     True
8     True
9     True
dtype: bool
pds = gpd.GeoDataFrame([polygon], columns=['geometry'], geometry = 'geometry', crs = {'init':'epsg:3309'}
t.intersects(pds.unary_union)
index_left
0    False
1    False
2     True
3    False
4     True
5     True
6     True
7     True
8     True
9     True
dtype: bool
t.intersects(pds.geometry.ix[0]) 
index_left
0    False
1    False
2     True
3    False
4     True
5     True
6     True
7     True
8     True
9     True
dtype: bool

And

from geopandas.tools import sjoin
sjoin(t, pds, how='inner')
    FID               geometry                         index_right
2    2  POINT (192530.4468407042 -465274.3331522844)            0
4    4  POINT (195314.0952219047 -468429.8881146177)            0
5    5  POINT (196320.0493430012 -479448.5877129403)            0
6    6  POINT (199780.1489409541 -446899.1169124786)            0
7    7  POINT (161437.1203666388 -440485.6615666826)            0
8    8  POINT (192413.1253873391 -465423.3888242305)            0
9    9  POINT (196039.9555616449 -473696.9492107867)            0

Related Solutions

[GIS] Geopandas Intersects Speed

For intersections in GeoPandas, it is better to use a spatial-join (see More Efficient Spatial join in Python without QGIS, ArcGIS, PostGIS, etc or How to efficiently determine which of thousands of polygons intersect with a linestring [duplicate]

import geopandas as gpd
parcels = gpd.read_file('parcels.shp')
roads = gpd.read_file('roads.shp')
intersections= gpd.sjoin(parcels, roads, how="left", op='intersects')
intersections.head()
    parcel                  geometry                         index_right  road
0  Parcel 1  POLYGON ((-0.6824583866837387 0.78233034571062...  0.0       Road 1  
1  Parcel 2  POLYGON ((-0.09859154929577452 0.3239436619718...  NaN       Nan
2  Parcel 3  POLYGON ((-0.9103713188220229 -0.1062740076824...  1.0       Road 2  
3  Parcel 4  POLYGON ((0.2266325224071704 0.620998719590268...  0.0       Road 1

With your solution

road = roads.unary_union
parcels['road_intersection'] = parcels.intersects(road)
parcels
    parcel                  geometry                         index_left  road_intersection  
0  Parcel 1  POLYGON ((-0.6824583866837387 0.78233034571062...   0         True 
1  Parcel 2  POLYGON ((-0.09859154929577452 0.3239436619718...   1         False
2  Parcel 3  POLYGON ((-0.9103713188220229 -0.1062740076824...   2         True
3  Parcel 4  POLYGON ((0.2266325224071704 0.620998719590268...   3         True

GeoPandas: How to Intersect MultiLinestring-Based Geodataframes

EDIT: The following answer only works for the 1-D case. To extend it to 2-D, you will need to parameterize your links by the along-road distance, and replace the x-coordinate with the parameterized length. However, I'm fairly confident this is doable much simpler with Geopandas.

It would be too hard to give hints in the comments, so here's a script that should give you what you want. It's not written for efficiency--probably you could get geopandas to do what you want with some finegaling, but here ya go. It's also not written very generally, but that could be done if you have more than one attribute.

import geopandas as gpd
from shapely.geometry import LineString
import numpy as np

slgeom = [[(0,0),(7,0)],[(7,0),(13,0)],[(13,0),(15,0)],[(15,0),(19,0)]]
geoms = []
for s in slgeom:
    geoms.append(LineString(s))
properti = ['a','b','c','d']
split_lines = gpd.GeoDataFrame(geometry=geoms)
split_lines['property'] = properti

olgeom = [[(0,0),(5,0)],[(5,0),(7,0), (10,0)],[(10,0),(13,0),(15,0)],[(15,0),(19,0)]]
geoms = []
for o in olgeom:
    geoms.append(LineString(o))
susc = [1,2,3,4]
original_lines = gpd.GeoDataFrame(geometry=geoms)
original_lines['susc'] = susc


# Do split lines
xs1 = []
attrib1 = []
for g, a in zip(split_lines.geometry.values, split_lines.property.values):
    x = g.coords.xy[0].tolist()
    xs1.extend(x)
    try:
        attrib1[-1] = a
    except:
        pass
    attrib1.extend([a for l in range(len(x))])

# Do originals
xs2 = []
attrib2 = []
for g, a in zip(original_lines.geometry.values, original_lines.susc.values):
    x = g.coords.xy[0].tolist()
    xs2.extend(x)
    try:
        attrib2[-1] = a
    except:
        pass
    attrib2.extend([a for l in range(len(x))])

# Create numpy array for sorting
allxs = list(set(xs1 + xs2))
x_forsort = []
a1_forsort = []
a2_forsort = []
for x in allxs:
    try:
        idx = xs1.index(x)
        a1_forsort.append(attrib1[idx])
    except:
        a1_forsort.append(None)
    try:
        idx = xs2.index(x)
        a2_forsort.append(attrib2[idx])
    except:
        a2_forsort.append(None)
forsort = np.transpose(np.array([allxs, a1_forsort, a2_forsort]))

# Now sort based on x value (1st column)
sorteds = forsort[forsort[:,0].argsort()]

# Work through the sorted lists to create segments with the appropriate attributes
# Store results in a dictionary
output = dict()
output['geometry'] = [] 
output['attrib_1'] = []
output['attrib_2'] = []
for i in range(len(sorteds)-1):

    # Store shapely linestring
    output['geometry'].append(LineString([(sorteds[i,0],0),(sorteds[i+1,0],0)]))

    # Store attributes
    if i == 0:
        output['attrib_1'].append(sorteds[i,1])
        output['attrib_2'].append(sorteds[i,2])
    else:
        if sorteds[i,1] is None:
            output['attrib_1'].append(output['attrib_1'][-1])
        else:
            output['attrib_1'].append(sorteds[i,1])

        if sorteds[i,2] is None:
            output['attrib_2'].append(output['attrib_2'][-1])
        else:
            output['attrib_2'].append(sorteds[i,2])

# Convert back to geopandas dataframe
out_gdf = gpd.GeoDataFrame(output)
out_gdf.crs = original_lines.crs

Result:

out_gdf
Out[185]: 
  attrib_1  attrib_2                 geometry
0        a         1    LINESTRING (0 0, 5 0)
1        a         2    LINESTRING (5 0, 7 0)
2        b         2   LINESTRING (7 0, 10 0)
3        b         3  LINESTRING (10 0, 13 0)
4        c         3  LINESTRING (13 0, 15 0)
5        d         4  LINESTRING (15 0, 19 0)

Best Answer

Related Solutions

[GIS] Geopandas Intersects Speed

GeoPandas: How to Intersect MultiLinestring-Based Geodataframes

Related Question