¿GDAL poligonaliza en python creando un polígono en blanco?

12

Tengo problemas para usar la función Polygonize en python. El ejemplo del libro de cocina para esto se puede encontrar aquí .

La parte relevante de mi código es:

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

Sé que la banda tiene información relevante, aquí hay un fragmento de 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.],

Cuando abro la tabla de atributos en QGIS, está vacía: Captura de pantalla QGIS

Editar:

La conversión funciona bien en QGIS usando Raster -> Conversion -> Polygonize tool

Captura de pantalla del ráster que se poligonalizará:

trama para poligonalizar

Y captura de pantalla de la conversión resultante de la herramienta QGIS:

ráster poligonalizado de la herramienta QGIS

Estoy usando la distribución Enthought en Windows 7, GDAL versión 1.10.0-3

El problema es que no puedo poligonalizar un ráster en Python usando GDAL y el ejemplo del libro de cocina, puedo poligonalizar este mismo ráster sin ningún problema en la GUI de QGIS

camdenl
fuente
¿Cómo se ve tu trama? ¿Realmente contiene polígonos? ¿Funciona si usas gdal_polygonize.py en su lugar?
BradHards
Editado para agregar capturas de pantalla del proceso de trabajo en QGIS
camdenl
¿Cuál es el problema real aquí?
Fezter
Problema específico agregado
camdenl
3
Tuve un problema similar (se creó un archivo de formas en blanco), y la creación del campo no ayudó. Lo que estaba haciendo mal era que no había cerrado el archivo de forma en mi código antes de llamar a polygonize. Lo cierras en tu ejemplo, solo estoy publicando esto para referencia de otros.
Stephanie

Respuestas:

19

El problema es que no estaba creando un campo para almacenar la banda ráster. Después de buscar en el archivo gdal_polygonize.py, me di cuenta de que esto no se hace automáticamente al llamar a gdal.Polygonize, que en su lugar usa la función que se encuentra aquí .

Aquí está el paso adicional necesario para crear un campo y escribir una banda en el campo:

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

Luego podemos escribir la banda en este campo, con un índice de 0:

gdal.Polygonize(band, None, outLayer, 0, [], callback=None )
camdenl
fuente
También estoy tratando de usar la función gdal.Polygonize () para obtener mi ráster como polígono en python. ¡Pero en la línea final está mostrando un error de tiempo de ejecución! ¿por qué?
Shiuli Pervin
Funciona bien con archivos ráster georreferenciados. El resultado son demasiados polígonos, pero solo quiero un gran polígono que muestre el contorno de la trama. ¿Alguien tiene alguna idea de cómo funciona al mismo tiempo?
Shiuli Pervin
Todavía obtengo un shapefile vacío, pero tengo filas en el archivo dbf. Por favor aclararme!
Satya Chandra
Acabo de tener este problema, pero en lugar de agregar un campo ficticio, puede ingresar un índice de -1. Vea aquí que el campo solo se agrega si index> = 0.
jon_two