2

Question relates to an earlier one: Can't retrieve geometry (Python, OGR, Rtree)

Im trying to simplify a dataset with 180k+ polygons and place the simplified geometries in a separate dataset. With the use of Check Geometry tool in ArcMap the original dataset contains one error:

at feat_id = -1 PROBLEM: could not find spatial index

But I can't select this feature and, for example, remove it. When I continue and apply .Simplify() to each feature iteratively in my script, it introduces 'null geometry' errors. Is my use of .Simplify() as shown below correct or not? I tried both .Simplify(10) and .Simplify(5). The former results in more errors (71) compared to latter (8).

if.os.path.exists('class.shp'):
    driver.DeleteDataSource('class.shp')
outFile = driver.CreateDataSource('class.shp')
outLayer = outFile.CreateLayer('class.shp', srs, geom_type = inLayer.GetLayerDefn().GetGeomType())

selectLayerDef = selectData.GetLayerDefn()
for i in range(0, selectLayerDef.GetFieldCount()):
outLayer.CreateField(selectLayerDef.GetFieldDefn(i))

outLayerDef = outLayer.GetLayerDefn()
for i in range(0, selectData.GetFeatureCount()):
    inFeature = selectData.GetNextFeature()
    outFeature = ogr.Feature(outLayerDef)
    for i in range(0, outLayerDef.GetFieldCount()):
        outFeature.SetField(outLayerDef.GetFieldDefn(i).GetNameRef(),
                            inFeature.GetField(i)
    geom = inFeature.GetGeometryRef()
    geomSimple = geom.Simplify(10)
    outFeature.SetGeometry(geomSimple.Clone())
    outFeature.SetField('FEAT_TYPE', 'LAV')
    outLayer.CreateFeature(outFeature)
    inFeature.Destroy()
    outFeature.Destroy()

UPDATE:

I changed the above code to add a check on if the simplified geometry has issues. When this is the case it would clone the original geometry instead of the simplified geometry. It's certainly not pretty but it seems to work.

xmin, xmax, ymin, ymax = 0
outLayerDef ´outLayer.GetLayerDefn()
for i in range (0, selectData.GetFeatureCount()):
    inFeature = selectData.GetNextFeature()
    outFeature = ogr.Feature(outLayerDef)
    for i in range(0, outLayerDef.GetFieldCount()):
        outFeature.SetField(outLayerDef.GetFieldDefn(i).GetNameRef(),
                            inFeature.GetField(i))
    testFeature = ogr.Feature(outLayerDef)
    for i in range(0, outLayerDef.GetFieldCount()):
        outFeature.SetField(outLayerDef.GetFieldDefn(i).GetNameRef(),
                            inFeature.GetField(i))
    geom = inFeature.GetGeometryRef()
    geomSimple = geom.Simplify(10)
    testFeature.SetGeometry(geomSimple.Clone())
    testGeom = testFeature.GetGeometryRef()
        if testGeom.GetEnvelope() == (xmin, xmax, ymin, ymax):
            outFeature.SetGeometry(inFeature.GetGeometryRef().Clone())
            outFeature.SetField('FEAT_TYPE', 'LAVBEBYG')
            outLayer.CreateFeature(outFeature)
            testFeature.Destroy()
            inFeature.Destroy()
            outFeature.Destroy()
        else:
            outFeature.SetGeometry(geomSimple.Clone())
            outFeature.SetField('FEAT_TYPE', 'LAVBEBYG')
            outLayer.CreateFeature(outFeature)
            testFeature.Destroy()
            inFeature.Destroy()
            outFeature.Destroy()

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.