GDAL RasterizeLayer ne brûle pas tous les polygones dans le raster?

12

J'essaie de graver un fichier de formes dans un raster en utilisant RasterizeLayer de GDAL. Je pré-crée une raster de zone d'intérêt à partir d'un fichier de formes différent, étant donné une taille de pixel spécifique. Cet AOI sert ensuite de base à toutes les rastérisations suivantes (même nombre de colonnes et de lignes, même projection et géotransformation).

Le problème se produit, cependant, lorsque je vais graver les formes dans leur propre raster, en fonction de la même taille de pixel et des mêmes projections. Le lien ci-dessous (pas assez de représentants pour publier l'image) montre le fichier de formes original en beige et le rose foncé où RasterizeLayer a gravé des données. Le rose clair correspond aux valeurs de nodata pour les données raster rose foncé. Le gris est l'AOI sur la base duquel la gravure du fichier de formes a été terminée.

Compte tenu de l'étendue des polygones de fichiers de formes, je m'attendrais à voir des valeurs de gravure dans les deux coins inférieurs, ainsi que les deux pixels sous les données qui s'affichent. Mais ce n'est évidemment pas le cas.

Image pour Problème - Gravures raster terminées

Voici le code que j'ai utilisé pour les générer. Toutes les formes ont été créées à l'aide de QGIS et ont toutes été créées dans la même projection. (Il convient de noter que le quadrillage dans l'image montrée était juste pour donner une idée de la taille de pixel que j'utilisais.)

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])

Est-ce un bogue dans GDAL, ou est-ce que RasterizeLayer brûle des données en se basant sur autre chose que la simple présence ou absence d'un polygone dans une zone de pixels spécifiée?

Les fichiers que j'utilisais se trouvent ici .

Alouette
la source
Pouvez-vous fournir un lien vers 'activity_3.shp' et 'AOI_Raster.tif'? Je veux voir si je peux recréer de mon côté.
Rich

Réponses:

10

J'ai joué avec GDALRasterizeLayers cette semaine et j'ai une assez bonne idée de ce qu'il fait. Par défaut, il pixellise un pixel si le centre du pixel se trouve dans le polygone. S'il n'y a rien au centre, il ne sera pas tramé, même s'il y a des parties d'un polygone dans les limites des pixels. Pour permettre à la pixellisation de fonctionner comme vous le souhaitez, essayez l'option "ALL_TOUCHED":

gdal.RasterizeLayer(raster_dataset, [1], shape_layer, None, None, [1], ['ALL_TOUCHED=TRUE'])
Mike T
la source
OUI! C'est apparemment ['ALL_TOUCHED=TRUE'], mais malheureusement, que seules les couches de polygones fixes. Mes couches de fichiers de formes ponctuelles sont toujours super bancales et affichent un pixel de l'endroit où elles sont placées.
Lark
Cela finit par ressembler à ça . Il est dans la même projection que les autres, et j'espérais que cela le réparerait aussi par magie, mais il semble juste brûler obstinément un pixel de l'endroit où il se trouve réellement.
Lark
Cela semble certainement digne des bogues, où le point de gravure est compensé par dx / 2 et dy / 2. Je me demande si ce bug persiste toujours avec le dernier tronc.
Mike T
Ce ne est pas! Cela fonctionne dans 1.9.0. Merci beaucoup!
Lark
1
Il y a aussi une assez bonne recette ici: gis.stackexchange.com/a/16916/9942
j08lue