TopologicalError: The operation 'GEOSIntersection_r' could not be performed
Asked Answered
B

2

8

Hi Guys, I am trying to map the district shapefile into assembly constituencies. I have shape files for [Both].Basically I have to map all the variables given at district level in census data to assembly constituency level. So I am following a pycon talk. Everything is working fine but I am getting error in get_intersection function.Error for that is TopologicalError: The operation 'GEOSIntersection_r' could not be performed. Likely cause is invalidity of the geometry <shapely.geometry.polygon.Polygon object at 0x7f460250ce10>.

I have tried to use both pygeos and rtree. There were some links which said that problem is in pygeos. So,I used rtree.But of no avail.Please help Thanks in advance.

Code tried by me is

constituencies=gpd.GeoDataFrame.from_file('/content/AC_All_Final.shp')
districts=gpd.GeoDataFrame.from_file('/content/2001_Dist.shp')
districts['AREA'] = districts.geometry.area
constituencies['AREA'] = constituencies.geometry.area
merged = gpd.sjoin(districts, constituencies).reset_index().rename(
    columns={'index': 'index_left'})

def get_intersection(row):
    left_geom = districts['geometry'][row['index_left']]
    right_geom = constituencies['geometry'][row['index_right']]
    return left_geom.intersection(right_geom)

***Error is at this point***
merged['geometry'] = merged.apply(get_intersection, axis=1)
merged['AREA'] = merged.geometry.area
Error trace is given below:
TopologyException: Input geom 1 is invalid: Ring Self-intersection at or near point 77.852561819157373 14.546596140487276 at 77.852561819157373 14.546596140487276
---------------------------------------------------------------------------
TopologicalError                          Traceback (most recent call last)
<ipython-input-17-8123669e025c> in <module>()
      4     return left_geom.intersection(right_geom)
      5 
----> 6 merged['geometry'] = merged.apply(get_intersection, axis=1)
      7 merged['AREA'] = merged.geometry.area

7 frames
/usr/local/lib/python3.6/dist-packages/shapely/topology.py in _check_topology(self, err, *geoms)
     36                     "The operation '%s' could not be performed. "
     37                     "Likely cause is invalidity of the geometry %s" % (
---> 38                         self.fn.__name__, repr(geom)))
     39         raise err
     40 

TopologicalError: The operation 'GEOSIntersection_r' could not be performed. Likely cause is invalidity of the geometry <shapely.geometry.polygon.Polygon object at 0x7f460250ce10>
Barraza answered 18/9, 2020 at 12:25 Comment(0)
C
13

The error message tells you exactly what is going on. Some of your geometries are not valid, so you have to make them valid before doing your apply. The simple trick, which works in most of the cases is using buffer(0).

merged['geometry'] = merged.buffer(0)

Since the issue is with geometry validity and is raised by GEOS, it does not matter if you use shapely/rtree backend or pygeos.

Chlorella answered 18/9, 2020 at 21:45 Comment(3)
@martinfleis.Still I am getting the same error after using your code. Error is TopologicalError: The operation 'GEOSIntersection_r' could not be performed. Likely cause is invalidity of the geometry <shapely.geometry.polygon.Polygon object at 0x7f5d03798518>Barraza
Then you have to check what is going on with that particular polygon. You can filter problematic ones via merged[~merged.is_valid].Chlorella
Thanks a lot @martinflies.Was able to solve using thisBarraza
F
11

Starting from shapely 1.8, there will be a make_valid method

However, currently shapely 1.8 isn't a stable release on pypi yet and you'd need to install an unstable one

pip3 install shapely==1.8a2 

Then you can make a shape valid as per

from shapely.validation import make_valid

valid_shape = make_valid(invalid_shape)

Note that the type of shape might change, for example from a Polygon to a MultiPolygon.

However, I'd argue it better to (1) properly avoid invalid shapes or (2) choose make_valid, since it's suggested by the shapely team, rather than the .buffer(0) work around.

Fontainebleau answered 16/7, 2021 at 10:2 Comment(1)
This one worked for me, thanks! (the buffer(0) method didn't)Bust

© 2022 - 2024 — McMap. All rights reserved.