6

I've searched, but can't find clear explanations how can I simply delete inner rings from existing polygons with python or extract outter rings and overwrite the geometry. Seems GDAL/OGR should be the simplest, but still I need some help...

Simple vision of mine would be:

dataSrc=driver.Open(shp)
layer=dataSrc.GetLayer()
for i in xrange(layer.GetFeatureCount()):
    feature = layer.GetFeature(i)
    geometry = feature.GetGeometryRef()
    if geometry.GetArea()<200:
        layer.DeleteFeature(i)
    if geometry.GetGeometryCount()>1:
        geometry=geometry.GetGeometryRef(0) #overwrite geom, getting the first ring outer
        feature.SetGeometry(geometry)
        layer.SetFeature(feature)

geometry.Destroy()

This doesn't loop through, it's just how I would like to see.

As I know the first ring or the first part of the geometry should be the outer ring, so I assumed I can just ask for it like that.Suppose that's not a problem.

So, probably the question is just about Updating Geometry. I open dataset withouth specific- Reading or Writing, so I assume it should let me write in it not just read. But if I do not want to have a new dataset, can I just update write like that or i should anywhays better create a new dataset?

Or maybe there is even more simple way to solve this with python?

3

I've found the solution with arcgisscripting, if anyone gona need to automate it.

It's really simple,.. just required some time to catch those holes :)

rows=gp.UpdateCursor(shp) #everything is done with one Cursor, but couple of arrays
row=rows.Next()
try:
    while row:
        if row.shape.area < 200:  #it's just deleting small polygons
            rows.DeleteRow(row)
        else:                     #part of cleaning from inner rings (donuts)
            geom=row.shape
            array=geom.GetPart(0)
            pnt=array.Next()
            newarray = gp.CreateObject("Array") #writing geometry to newArray       
            while pnt:
                newarray.add(pnt)
                pnt=array.Next()
                newpnt=newarray.Next()
                if not pnt:                     #this is the HOLE! 
                    break                       #when detect - break the loop!
            row.shape=newarray
            rows.UpdateRow(row)
        row=rows.Next()

    del row
    del rows
except:
    gp.GetMessages()
2

The part about deleting the features. Seems it's often not allowed. Found, there is a possibilty to test it: http://www.gdal.org/ogr/classOGRLayer.html#aeedbda1a62f9b89b8e5f24332cf22286

layer.TestCapability("DeleteField")

and so I get False..

And accessing just the first ring goes well, I've checked with area calculations.

The solution then is just to CREATE a new file with features I need from old file.. Correct me if I am misstaken.

  • You are correct. Updatability of layers varies over formats, but creation of a new dataset always works. – sgillies Sep 17 '12 at 14:50
  • Seems I at such point I am comming back to arcgisscripting.. – najuste Sep 17 '12 at 15:02
1

You can use the ArcGIS Union and Dissolve tools to fill gaps. Make sure to use the NO_GAPS option in the union tool.

import arcpy

in_shp = r'C:/path/to/your/shapefile.shp'
out_shp = r'C:/path/to/your/out_shapefile.shp'

arcpy.Union_analysis(in_shp, "in_memory/union", gaps = 'NO_GAPS')
arcpy.Dissolve_management("in_memory/union", out_shp)

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.