10

I'm trying to intersect two polygons based on the following link, intersecting two shapefiles from Python or command line

But problem occured that shows a "ValueError: Geometry column cannot contain mutiple geometry types when writing to file."

I'm new on these free open source GIS tool manipulation, can anyone give me some idea? My source code is like this:

    from shapely.geometry import shape,Polygon,MultiPolygon,mapping
    import geopandas as gpd

    g1 = gpd.GeoDataFrame.from_file("./origin_test.shp")
    g2 = gpd.GeoDataFrame.from_file("./ref_test.shp")

    data=[]
    for index, orig in g1.iterrows():
        for index2, ref in g2.iterrows():
            if ref['geometry'].intersects(orig['geometry']):
               owdspd=orig['wdspd']
               data.append({'geometry':ref['geometry'].intersection(orig['geometry']),'wdspd':owdspd})

    data = data.set_geometry()
    df = gpd.GeoDataFrame(data,columns=['geometry','wdspd'])
    df.to_file('./intersection.shp')

1 Answer 1

19

1) The problem is that the intersection of two polygons is not always a polygon

enter image description here

import geopandas as gp
poly1 = gp.read_file("poly_origin.shp")
poly2 = gp.read_file("poly_test.shp")
data = []
for index, orig in poly1.iterrows():
    for index2, ref in poly2.iterrows():      
        if ref['geometry'].intersects(orig['geometry']): 
         owdspd=orig['id']
         data.append({'geometry':ref['geometry'].intersection(orig['geometry']),'wdspd':owdspd})

for geom in data: 
   print geom
{'geometry': <shapely.geometry.polygon.Polygon object at 0x10aab9c10>, 'wdspd': 1}
{'geometry': <shapely.geometry.polygon.Polygon object at 0x10aab9e50>, 'wdspd': 2}
{'geometry': <shapely.geometry.linestring.LineString object at 0x10aacc250>, 'wdspd': 3}

The result is 2 Polygons and a LineString therefore the error

"ValueError: Geometry column cannot contain mutiple geometry types when writing to file"

2) A solution is to convert all the geometries to the same type

from shapely.geometry import LineString
for i in data:
 geom = i['geometry']
 if geom.geom_type=='Polygon':
      i['geometry']=  LineString(list(geom.exterior.coords))

Now, you can convert to a GeoDataFrame

df = gp.GeoDataFrame(data,columns=['geometry','wdspd'])
df                          geometry                   wdspd
0  LINESTRING (96.42123294751703 -44.632201588700...      1
1  LINESTRING (282.773470482277 -127.548992317643...      2
2  LINESTRING (144.3726591760299 -221.86142322097...      3

3) Or to use GeoPandas Overlay

inter = gp.overlay(poly1, poly2, how='intersection')
inter
     id  id_2                geometry
 0   2     2  POLYGON ((282.773470482277 -127.5489923176438,...
 1   1     1  POLYGON ((51.0880537778328 -116.6195085418168,...

Export the result

df.to_file("linestring_intersection.shp")

enter image description here

8
  • Thank you sooo much, it's really helpful. According to your idea(method 2 and 3), i tried to output the data as shape file in both way with following code: df = gpd.GeoDataFrame(data,columns=['geometry','wdspd']) df.to_file('./wdspd_intersection.shp') f.head() But it failed. Could you mind tell me how to export the result? Thanks...
    – King-Zhao
    Mar 2, 2017 at 10:33
  • look above in Export the result
    – gene
    Mar 2, 2017 at 10:43
  • Yeah, thanks @gene . I've tried, the linestring is a line, i tried the overlay way, it shows an error like this self.session.start(self, **kwargs) File "fiona\ogrext.pyx", line 942, in fiona.ogrext.WritingSession.start (fiona/ogrext.c:16386) ValueError: Null layer:
    – King-Zhao
    Mar 2, 2017 at 10:48
  • 1
    In addition, how to exclude geometry types other than Ploygon? @gene
    – King-Zhao
    Mar 2, 2017 at 11:50
  • 1
    It would be nice to exclude line geometry types... probably possible with a simple test but I do not know how.
    – r4gt4g
    Jan 27, 2020 at 18:06

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.