3

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')
6

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

  • 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 '17 at 10:33
  • look above in Export the result – gene Mar 2 '17 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 '17 at 10:48
  • If possible i would like to export in polygons this could make my work more simple. – King-Zhao Mar 2 '17 at 10:49
  • If the result contains other geometries than polygons in the result, you cannot (or try with the solution 3) – gene Mar 2 '17 at 10:51

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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