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:
Editar:
La conversión funciona bien en QGIS usando Raster -> Conversion -> Polygonize tool
Captura de pantalla del ráster que se poligonalizará:
Y captura de pantalla de la conversión resultante 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
Respuestas:
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:
Luego podemos escribir la banda en este campo, con un índice de 0:
fuente