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()
Internally, PostGIS stores geometries in a binary specification, but it is queried and viewed outside as a hex-encoded string. There are two popular variations of well-known binary (WKB):
- EWKB (via
ST_AsEWKB) - an extended WKB specification designed by PostGIS.
- OGC WKB (via
ST_AsBinary) - specified by the OGC and ISO. For a while it was 2D-only, but later extended to support Z, M and ZM geometries.
The two specifications are the same for 2D geometries, but are different for higher-order geometries with Z, M and ZM coordinates.
Older versions of GDAL/OGR (1.x) only understand the EWKB for 3D geometries, so for these I recommend using ST_AsEWKB. (But if you only have 2D geometries, either format are fine). For example:
import psycopg2
from osgeo import ogr
ogr.UseExceptions()
conn = psycopg2.connect('dbname=postgis user=postgres')
curs = conn.cursor()
curs.execute("select ST_AsEWKB('POINT Z (1 2 3)'::geometry) AS g")
b = bytes(curs.fetchone()[0])
print(b.encode('hex')) # 0101000080000000000000f03f00000000000000400000000000000840
g = ogr.CreateGeometryFromWkb(b)
print(g.ExportToWkt()) # POINT (1 2 3)
curs.execute("select ST_AsBinary('POINT Z (1 2 3)'::geometry) AS g")
b = bytes(curs.fetchone()[0])
print(b.encode('hex')) # 01e9030000000000000000f03f00000000000000400000000000000840
g = ogr.CreateGeometryFromWkb(b)
# RuntimeError: OGR Error: Unsupported geometry type
Also, note that older GDAL/OGR versions do not support M coordinates, and these will be parsed but ignored.
With GDAL 2.0 and more recent, ISO WKT/WKB is supported. This means that CreateGeometryFromWkb can read either WKB flavour (without specifying) and ExportToIsoWkt() shows output with a modern WKT syntax.
import psycopg2
from osgeo import ogr
ogr.UseExceptions()
conn = psycopg2.connect('dbname=postgis user=postgres')
curs = conn.cursor()
curs.execute("select ST_AsEWKB('POINT Z (1 2 3)'::geometry) AS g")
b = bytes(curs.fetchone()[0])
print(b.encode('hex')) # 0101000080000000000000f03f00000000000000400000000000000840
g = ogr.CreateGeometryFromWkb(b)
print(g.ExportToIsoWkt()) # POINT Z (1 2 3)
curs.execute("select ST_AsBinary('POINT Z (1 2 3)'::geometry) AS g")
b = bytes(curs.fetchone()[0])
print(b.encode('hex')) # 01e9030000000000000000f03f00000000000000400000000000000840
g = ogr.CreateGeometryFromWkb(b)
print(g.ExportToIsoWkt()) # POINT Z (1 2 3)
Additionally, GDAL 2.1 or later will create/export WKT/WKB with M or ZM coordinates as expected.
Best Answer
Don't expect OGR to be Pythonic. In fact, there is a large list of gotchas documented to describe these. Essentially, all of the variables are references to objects described in C++, so they often don't make sense with workflows in Python. Modules like
copyonly work on Python objects, which is why your example doesn't work as expected. You need to keep references to the features and the geometries alive, not just the geometry, otherwise it will crash.For example, here is something that will calculate the union of current and previous features, showing the area of all three geometries:
But if, for example, outside the loop did
gNew = gCurrent.Union(gPrevious)it will crash since the references to the features have changed.