GDAL RasterizeLayer no quema todos los polígonos en ráster?

12

Estoy tratando de grabar un shapefile en un ráster usando RasterizeLayer de GDAL. Precreo un ráster de área de interés desde un archivo de formas diferente, dado un tamaño de píxel específico. Este AOI sirve entonces como base para todas las siguientes rasterizaciones (mismo número de columnas y filas, misma proyección y geotransformación).

Sin embargo, el problema ocurre cuando voy a grabar las formas en su propio ráster, en función del mismo tamaño de píxel y las mismas proyecciones. El siguiente enlace (no tiene suficiente representante para publicar la imagen), muestra el archivo de forma original en bronceado y el rosa oscuro donde RasterizeLayer ha quemado datos. El rosa claro son valores de nodata para los datos ráster de color rosa oscuro. El gris es el AOI en función del cual se completó la grabación del archivo de forma.

Dada la extensión de los polígonos del archivo de forma, esperaría ver valores de grabación en las dos esquinas inferiores, así como los dos píxeles debajo de los datos que se muestran. Obviamente, sin embargo, este no es el caso.

Imagen del problema: quemaduras ráster terminadas

El siguiente es el código que usé para generar estos. Todas las formas fueron creadas usando QGIS, y todas fueron creadas en la misma proyección. (Cabe señalar que la cuadrícula en la imagen mostrada fue solo para dar una idea del tamaño de píxel que estaba usando).

from osgeo import ogr
from osgeo import gdal

aoi_uri = 'AOI_Raster.tif'
aoi_raster = gdal.Open(aoi_uri)

def new_raster_from_base(base, outputURI, format, nodata, datatype):

    cols = base.RasterXSize
    rows = base.RasterYSize
    projection = base.GetProjection()
    geotransform = base.GetGeoTransform()
    bands = base.RasterCount

    driver = gdal.GetDriverByName(format)

    new_raster = driver.Create(str(outputURI), cols, rows, bands, datatype)
    new_raster.SetProjection(projection)
    new_raster.SetGeoTransform(geotransform)

    for i in range(bands):
        new_raster.GetRasterBand(i + 1).SetNoDataValue(nodata)
        new_raster.GetRasterBand(i + 1).Fill(nodata)

    return new_raster

shape_uri = 'activity_3.shp'
shape_datasource = ogr.Open(shape_uri)
shape_layer = shape_datasource.GetLayer()

raster_out = 'new_raster.tif'

raster_dataset = new_raster_from_base(aoi_raster, raster_out, 'GTiff',
                                -1, gdal.GDT_Int32)
band = raster_dataset.GetRasterBand(1)
nodata = band.GetNoDataValue()

band.Fill(nodata)

gdal.RasterizeLayer(raster_dataset, [1], shape_layer, burn_values=[1])

¿Es esto un error en GDAL, o RasterizeLayer está quemando datos basados ​​en algo más que la presencia o la falta de un polígono dentro de un área de píxeles especificada?

Los archivos que estaba usando se pueden encontrar aquí .

Alondra
fuente
¿Puede proporcionar un enlace a 'activity_3.shp' y 'AOI_Raster.tif'? Quiero ver si puedo recrearme por mi parte.
Rico

Respuestas:

10

He estado jugando con GDALRasterizeLayers esta semana y tengo una muy buena idea de lo que está haciendo. Por defecto, rasterizará un píxel si el centro del píxel está dentro del polígono. Si no hay nada en el centro, no se rasterizará, incluso si hay partes de un polígono dentro de los límites de píxeles. Para permitir que la rasterización funcione de la manera que desea, pruebe la opción "ALL_TOUCHED":

gdal.RasterizeLayer(raster_dataset, [1], shape_layer, None, None, [1], ['ALL_TOUCHED=TRUE'])
Mike T
fuente
¡SI! Aparentemente ['ALL_TOUCHED=TRUE'], aunque desafortunadamente, solo se arreglaron las capas poligonales. Las capas de mi archivo de formas puntuales todavía están súper inestables y aparecen un píxel desde donde están colocadas.
Alondra
Se termina pareciéndose a esto . Está en la misma proyección que los demás, y esperaba que esto también lo arreglara mágicamente de alguna manera, pero solo parece quemar tercamente un píxel de donde realmente se encuentra.
Alondra
Eso ciertamente parece digno de error, donde el punto de quemado se compensa con dx / 2 y dy / 2. Me pregunto si ese error aún persiste con la última troncal.
Mike T
¡No es asi! Funciona en 1.9.0. ¡Muchas gracias!
Alondra
1
También hay una buena receta aquí: gis.stackexchange.com/a/16916/9942
j08lue