GDAL-Polygonisierung in Python Erstellen eines leeren Polygons?

12

Ich habe Probleme bei der Verwendung der Polygonize-Funktion in Python. Das Kochbuchbeispiel hierfür finden Sie hier .

Der relevante Teil meines Codes ist:

sourceRaster = gdal.Open('myraster.tif')
band = sourceRaster.GetRasterBand(1)
bandArray = band.ReadAsArray()
outShapefile = "polygonized"
driver = ogr.GetDriverByName("ESRI Shapefile")
if os.path.exists(outShapefile+".shp"):
    driver.DeleteDataSource(outShapefile+".shp")
outDatasource = driver.CreateDataSource(outShapefile+ ".shp")
outLayer = outDatasource.CreateLayer("polygonized", srs=None)
gdal.Polygonize( band, None, outLayer, -1, [], callback=None )
outDatasource.Destroy()
sourceRaster = None

Ich weiß, dass die Band relevante Informationen hat, hier ein Ausschnitt aus bandArray:

array([[ 4.,  4.,  3.,  3.,  3.,  2.,  2.,  2.,  2.,  3.,  3.,  3.,  3.,
         3.,  3.,  3.,  3.,  4.,  4.,  4.,  4.,  4.,  4.,  4.,  4.,  4.,
         4.,  4.,  4.,  4.],

Wenn ich die Attributtabelle in QGIS öffne, ist sie leer: QGIS-Bildschirmaufnahme

Bearbeiten:

Die Konvertierung funktioniert in QGIS mit dem Werkzeug Raster -> Konvertierung -> Polygonisieren einwandfrei

Screenshot des zu polygonisierenden Rasters:

Raster zum Polygonisieren

Und Screenshot der resultierenden Konvertierung vom QGIS-Tool:

polygonisiertes Raster aus dem QGIS-Tool

Ich verwende die Enthought-Distribution unter Windows 7, GDAL Version 1.10.0-3

Das Problem ist, dass ich ein Raster in Python mit GDAL und dem Kochbuchbeispiel nicht polygonisieren kann. Ich kann dasselbe Raster ohne Probleme in der QGIS-GUI polygonisieren

camdenl
quelle
Wie sieht dein Raster aus? Enthält es wirklich Polygone? Funktioniert es, wenn Sie stattdessen gdal_polygonize.py verwenden?
BradHards
Bearbeitet, um Screenshots des Arbeitsprozesses in QGIS
hinzuzufügen
Was ist das eigentliche Problem hier?
Fezter
Spezifisches Problem
hinzugefügt
3
Ich hatte ein ähnliches Problem (es wurde ein leeres Shapefile erstellt), und das Erstellen des Felds half nicht. Was ich falsch gemacht habe war, dass ich das Shapefile in meinem Code nicht geschlossen hatte, bevor ich polygonize aufgerufen hatte. Sie schließen es in Ihrem Beispiel. Ich poste dies nur als Referenz für andere.
Stephanie

Antworten:

19

Das Problem ist, dass ich kein Feld zum Speichern des Rasterbands erstellt habe. Nachdem ich die Datei gdal_polygonize.py durchgesehen hatte, stellte ich fest, dass dies beim Aufruf von gdal.Polygonize, das stattdessen die hier gefundene Funktion verwendet, nicht automatisch erfolgt .

Hier ist der zusätzliche Schritt, der erforderlich ist, um ein Feld zu erstellen und ein Band in das Feld zu schreiben:

newField = ogr.FieldDefn('MYFLD', ogr.OFTInteger)
outLayer.CreateField(newField)

Wir können dann das Band mit einem Index von 0 in dieses Feld schreiben:

gdal.Polygonize(band, None, outLayer, 0, [], callback=None )
camdenl
quelle
Ich versuche auch, die Funktion gdal.Polygonize () zu verwenden, um mein Raster als Polygon in Python abzurufen. Aber am Ende wird ein Laufzeitfehler angezeigt !! Warum?
Shiuli Pervin
Es funktioniert gut mit georeferenzierten Rasterdateien. Das Ergebnis sind zu viele Polygone, aber ich möchte nur ein großes Polygon, das den Umriss des Rasters zeigt. Hat jemand eine Idee, wie es gleichzeitig funktioniert?
Shiuli Pervin
Ich bekomme immer noch ein leeres Shapefile, aber ich habe Zeilen in der DBF-Datei. Bitte klären Sie mich!
Satya Chandra
Ich hatte gerade dieses Problem, aber anstatt ein Dummy-Feld hinzuzufügen, können Sie einen Index von -1 eingeben. Sehen Sie hier, dass das Feld nur hinzugefügt wird, wenn index> = 0 ist.
jon_two