Monday, 28 November 2016

'.Simplify' causes 'null geometry' errors (Python, OGR)


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()


No comments:

Post a Comment

arcpy - Changing output name when exporting data driven pages to JPG?

Is there a way to save the output JPG, changing the output file name to the page name, instead of page number? I mean changing the script fo...