GDAL RasterizeLayer brennt nicht alle Polygone in Raster?

12

Ich versuche, ein Shapefile mit GDALs RasterizeLayer in ein Raster zu brennen. Ich erstelle ein Interessenbereich-Raster aus einem anderen Shapefile bei einer bestimmten Pixelgröße. Dieser AOI dient dann als Basis für alle folgenden Rasterisierungen (gleiche Anzahl von Spalten und Zeilen, gleiche Projektion und Geotransformation).

Das Problem tritt jedoch auf, wenn ich die Formen auf der Grundlage derselben Pixelgröße und Projektionen in ein eigenes Raster brenne. Der unten stehende Link (hat nicht genügend Wiederholungen, um das Bild zu veröffentlichen) zeigt das ursprüngliche Shapefile in Tan und das dunkle Rosa, in dem RasterizeLayer Daten gebrannt hat. Das Hellrosa ist ein Knotenwert für die dunkelrosa Rasterdaten. Das Grau ist der AOI, auf dessen Grundlage der Shapefile-Brennvorgang abgeschlossen wurde.

Angesichts der Ausmaße der Shapefile-Polygone würde ich erwarten, dass in den unteren beiden Ecken Brennwerte sowie die beiden Pixel unter den angezeigten Daten angezeigt werden. Offensichtlich ist dies jedoch nicht der Fall.

Bild für problematische Rasterbrände

Wie folgt ist der Code, den ich verwendet habe, um diese zu generieren. Alle Formen wurden mit QGIS erstellt und in derselben Projektion. (Es sollte beachtet werden, dass das Gitter in dem gezeigten Bild nur eine Vorstellung von der Pixelgröße gab, die ich verwendet habe.)

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

Handelt es sich um einen Fehler in GDAL, oder brennt RasterizeLayer Daten, die nicht nur auf dem Vorhandensein oder Fehlen eines Polygons in einem bestimmten Pixelbereich beruhen?

Die von mir verwendeten Dateien finden Sie hier .

Lerche
quelle
Können Sie einen Link zu 'activity_3.shp' und 'AOI_Raster.tif' bereitstellen? Ich möchte sehen, ob ich mich auf meinem Ende wieder herstellen kann.
Rich

Antworten:

10

Ich habe diese Woche mit GDALRasterizeLayers gespielt und habe eine ziemlich gute Vorstellung davon, was es tut. Standardmäßig wird ein Pixel gerastert, wenn sich die Pixelmitte innerhalb des Polygons befindet. Wenn sich nichts in der Mitte befindet, wird es nicht gerastert, auch wenn sich Teile eines Polygons innerhalb der Pixelgrenzen befinden. Verwenden Sie die Option "ALL_TOUCHED", damit das Rastern wie gewünscht funktioniert:

gdal.RasterizeLayer(raster_dataset, [1], shape_layer, None, None, [1], ['ALL_TOUCHED=TRUE'])
Mike T
quelle
JA! Es ist anscheinend ['ALL_TOUCHED=TRUE']leider so, dass nur Polygonebenen fixiert werden. Meine Punkt-Shapefile-Ebenen sind immer noch super-wackelig und erscheinen einen Pixel über der Stelle, an der sie platziert sind.
Lark
Es sieht am Ende so aus . Es ist in der gleichen Projektion wie die anderen, und ich hatte gehofft, dass dies es irgendwie auch auf magische Weise beheben würde, aber es scheint nur hartnäckig einen Pixel von der Stelle wegzubrennen, an der es sich tatsächlich befindet.
Lark
Das sieht sicherlich fehlerhaft aus, wenn der Brennpunkt um dx / 2 und dy / 2 versetzt ist. Ich frage mich, ob dieser Fehler mit dem neuesten Kofferraum immer noch besteht.
Mike T
Es tut nicht! Es funktioniert in 1.9.0. Vielen Dank!
Lark
1
Hier gibt es auch ein gutes Rezept: gis.stackexchange.com/a/16916/9942
j08lue