0

Dear all am trying to perform a union between two shapefiles ( for the moment ) - line 88-

 oLayer.Union(oLayer2,oLayer3)

and I'm getting stack with this error

Traceback (most recent call last):
  File "/home/Uka/PycharmProjects/histcult/histoire.py", line 88, in <module>
    oLayer.Union(oLayer2,oLayer3)
  File "/usr/lib/python3/dist-packages/osgeo/ogr.py", line 2777, in Union
    return _ogr.Layer_Union(self, *args, **kwargs)
RuntimeError: Attempt to write non-polygon (LINESTRING) geometry to POLYGON type shapefile.

and this is my code :

import os
from osgeo import ogr
import osgeo.osr as osr

ogr.UseExceptions()
osr.UseExceptions()

# setting up the shapefile  driver
shpDriver = ogr.GetDriverByName("ESRI Shapefile")

# creating the data source
OutputDS = shpDriver.CreateDataSource("test_sample/output_histoire.shp")
OutputDS2 = shpDriver.CreateDataSource("test_sample/output_histoire2.shp")
OutputDS3 = shpDriver.CreateDataSource("test_sample/output.shp")

# create the spatial reference, WGS84
srs = osr.SpatialReference()
srs.ImportFromEPSG(2154)

# creating output layer the layer
oLayer = OutputDS.CreateLayer("culture", srs, ogr.wkbPolygon)
oLayer2 = OutputDS.CreateLayer("culture2", srs, ogr.wkbPolygon)
oLayer3 = OutputDS.CreateLayer("output", srs, ogr.wkbPolygon)

#Output Field names
cultureField = ogr.FieldDefn("cultPr", ogr.OFTString)
cultureField.SetWidth(30)
oLayer.CreateField(cultureField)
oLayer.CreateField(ogr.FieldDefn("ID_C", ogr.OFTInteger))


#input layers
cultureUrl = "test_sample/culture.shp"
parcelleUrl = "test_sample/parcelles2016test.shp"


cultureDtSource = shpDriver.Open(cultureUrl, 0) # 0 means read-only. 1 means writeable.
parceDtSource = shpDriver.Open(parcelleUrl, 0) #
cultureLayer = cultureDtSource.GetLayer(0)
parcelleLayer = parceDtSource.GetLayer(0)
parcelleDef = parcelleLayer.GetLayerDefn()




# Check to see if shapefile is found.
if cultureDtSource is None  :
    print ('Could not open %s' % (cultureUrl))
elif parceDtSource is None :
    print('Could not open %s' % (parcelleUrl))
else:
    print ('Opened %s' % (cultureUrl))
    print('Opened %s' % (parcelleUrl))

    cultureFC = cultureDtSource.GetLayer().GetFeatureCount()
    parcelleFC = parceDtSource.GetLayer().GetFeatureCount()
    print ("Number of features in %s: %d" % (os.path.basename(cultureUrl), cultureFC))
    print("Number of features in %s: %d" % (os.path.basename(parcelleUrl), parcelleFC))
    for i in range(parcelleDef.GetFieldCount()): # iterate on parcelle features to get the names
        print (parcelleDef.GetFieldDefn(i).GetName())
    nbrIntersection = 0

    for featureP in parcelleLayer: # iterate on parcelle features to get the the geometry


        geomP = featureP.GetGeometryRef()
        for featureC in cultureLayer: # iterate on culture  features to get check the intersection
            geomC = featureC.GetGeometryRef()

            intersection = geomC.Intersection(geomP)
            if not intersection.IsEmpty():
                # create the feature
                outputFeature = ogr.Feature(oLayer.GetLayerDefn())
                # Set the attributes
                outputFeature.SetField("cultPr", featureC.GetField("cult"))
                # Set the feature geometry
                print(intersection.GetGeometryName())
                outputFeature.SetGeometry(intersection)
                # Create the feature in the layer (shapefile)
                oLayer.CreateFeature(outputFeature)
                # Dereference the feature
                feature = None
                nbrIntersection = nbrIntersection + 1
        cultureLayer.ResetReading()
    parcelleLayer.ResetReading()
    parcelleLayer.Erase(oLayer,oLayer2);
    # add union between oLayer and oLayer2
    oLayer.Union(oLayer2,oLayer3)
    print('you got ', nbrIntersection, ' intersections')
    # Save and close the data source
    OutputDS = None
    OutputDS2 = None
    OutputDS3 = None
  • I guess that intersection returns both polygons and linestrings and you can't save them into one shapefile. Either sort output geometries by geometry type, filter out lines because I guess you are not interested in them, or use some format that supports mixed geometry types for output. – user30184 Oct 11 '17 at 20:10
  • i used the line 77 to see what are theintersection reslut's types : print(intersection.GetGeometryName()) and the result was that they are all polygons – Abdelaziz Oct 12 '17 at 6:55
  • Then keep on searching what part of the code does return lines because that is the issue Attempt to write non-polygon (LINESTRING) geometry to POLYGON type shapefile. Quick test is to write output into GML and find if the result set really contains linestrings as I suppose. – user30184 Oct 12 '17 at 8:37
  • after digging i didn't find any LINESTRING , so i created my own UNION methode and used it . thanks – Abdelaziz Oct 19 '17 at 8:55
  • If GDAL claims that it has found a linestring but there isn't any then it can be a GDAL bug. Create and put available a small zip that contains source data and python code and write a mail to gdal-dev mailing list. – user30184 Oct 19 '17 at 9:22

Your Answer

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

Browse other questions tagged or ask your own question.