¿Excluir la extensión al poligonizar un archivo raster usando Python?

Question

Quiero poligonizar un archivo raster. El archivo rasterizado tiene el siguiente aspecto.

Aquí está el archivo rasterizado

Aquí están mis códigos en Python para poligonizar el archivo raster. La fuente de los códigos de Poligonizar una banda rasterizada

#raster to polygons

from osgeo import gdal, ogr

#  get raster datasource

src_ds = gdal.Open( "C:\\Users\\select3.tif" )

srcband = src_ds.GetRasterBand(1)

#  create output datasource

dst_layername = "output_select3"
drv = ogr.GetDriverByName("ESRI Shapefile")
dst_ds = drv.CreateDataSource( dst_layername + ".shp" )
dst_layer = dst_ds.CreateLayer(dst_layername, srs = None )

gdal.Polygonize( srcband, None, dst_layer, -1, [], callback=None )
#The following two steps are very important, otherwise will generate error like (Coordinate not defined)
dst_ds.Destroy()
src_ds=None

El resultado es el siguiente.

Y aquí está el shapefile de salida

No estoy satisfecho con esto, porque sólo quiero tener las áreas poligonales (Marca de verificación) y excluir la extensión (marca de la cruz) de los polígonos. ¿Alguien podría decirme cómo arreglar esto en python?

Answer 1

Esto puede lograrse aplicando una máscara como segundo argumento en la función Polygonize, como se indica en el GDAL documentación. La máscara tiene que ser una capa de trama separada, que tiene 0 donde no quiere que el algoritmo poligonee.

Con sus datos, siga estos pasos para implementarlos:

1) Ejecute la calculadora de trama en su trama original ("select3.tif") con esta expresión para establecer el valor de las áreas que desee en 1, y todo el resto en 0:

2) Edita tu código con estas líneas:

from osgeo import gdal, ogr

#  get raster datasource
src_ds = gdal.Open( "C:\\Users\\select3.tif" )
srcband = src_ds.GetRasterBand(1)

mask_ds = gdal.Open( "C:\\Users\\mask.tif" ) #Path to the mask layer generated above
maskband = mask_ds.GetRasterBand(1) 

#  create output datasource
dst_layername = "output_select3"
drv = ogr.GetDriverByName("ESRI Shapefile")
dst_ds = drv.CreateDataSource( dst_layername + ".shp" )
dst_layer = dst_ds.CreateLayer(dst_layername, srs = None )

gdal.Polygonize( srcband, maskband, dst_layer, -1, [], callback=None ) #Mask's band as second argument in Polygonize function

dst_ds.Destroy()
src_ds=None

Resultado:

Answer 2

Normalmente no se necesita una máscara para esa tarea, pero los valores RGB de las áreas del polígono son iguales a 0. Así que hay que cambiarlos a un valor (por ejemplo 125) y cambiar los otros a 0. Dado que todas las bandas tienen el mismo valor, es suficiente cambiar una banda para poligonizar el raster como lo hiciste.

from osgeo import gdal, ogr

src_ds = gdal.Open( "C:\\Users\\select3.tif" )
srcband = src_ds.GetRasterBand(1)

# Invert pixel value to use srcband as a mask.
numpy_band = srcband.ReadAsArray() ##
numpy_band[numpy_band<125] = 125   ##
numpy_band[numpy_band>125] = 0     ##
srcband.WriteArray(numpy_band)     ##

dst_layername = "C:\\Users\\mask.tif"
drv = ogr.GetDriverByName("ESRI Shapefile")
dst_ds = drv.CreateDataSource( dst_layername + ".shp" )
dst_layer = dst_ds.CreateLayer(dst_layername, srs = None )

gdal.Polygonize( srcband, srcband, dst_layer, -1, [], callback=None )  ##
dst_ds.Destroy()
src_ds=None