GDAL RasterizeLayer nie wypala wszystkich wielokątów do rastra?

12

Próbuję wypalić plik kształtu do rastra przy użyciu RasterizeLayer GDAL. Wstępnie tworzę obszar zainteresowania rastrem z innego pliku kształtu, biorąc pod uwagę określony rozmiar piksela. Ten AOI służy następnie jako podstawa dla wszystkich kolejnych rasteryzacji (ta sama liczba kolumn i wierszy, ta sama projekcja i geotransforma).

Problem pojawia się jednak, gdy idę wypalić kształty do ich własnego rastra, w oparciu o ten sam rozmiar pikseli i rzuty. Link poniżej (nie ma wystarczającej liczby powtórzeń, aby opublikować obraz), pokazuje oryginalny plik kształtu w jasnobrązowym kolorze i ciemnoróżowy, w którym RasterizeLayer wypalił dane. Jasnoróżowy to wartości nodata dla ciemnoróżowych danych rastrowych. Szary to AOI, na podstawie którego zakończono zapisywanie pliku kształtu.

Biorąc pod uwagę zakres wielokątów kształtu, spodziewałbym się zobaczyć wartości wypalenia w dwóch dolnych rogach, a także dwa piksele pod danymi, które się pokazują. Oczywiście jednak tak nie jest.

Obraz problemu zakończonego oparzeniem rastrowym

Poniżej znajduje się kod, którego użyłem do ich wygenerowania. Wszystkie kształty zostały utworzone za pomocą QGIS i wszystkie zostały utworzone w tej samej projekcji. (Należy zauważyć, że siatka na pokazanym zdjęciu miała jedynie na celu wyobrażenie o rozmiarze piksela, którego używałem.)

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

Czy to błąd w GDAL, czy też RasterizeLayer wypala dane w oparciu o coś innego niż tylko obecność lub brak wielokąta w określonym obszarze pikseli?

Pliki, których używałem, można znaleźć tutaj .

Skowronek
źródło
Czy możesz podać link do „activity_3.shp” i „AOI_Raster.tif”? Chcę zobaczyć, czy mogę odtworzyć po mojej stronie.
Bogaty

Odpowiedzi:

10

W tym tygodniu grałem z GDALRasterizeLayers i mam całkiem niezłe pojęcie o tym, co robi. Domyślnie rasteryzuje piksel, jeśli środek piksela znajduje się w wielokącie. Jeśli w środku nie ma nic, nie zostanie rasteryzowane, nawet jeśli w granicach pikseli znajdują się części wielokąta. Aby umożliwić rasteryzację działanie w zamierzony sposób, wypróbuj opcję „ALL_TOUCHED”:

gdal.RasterizeLayer(raster_dataset, [1], shape_layer, None, None, [1], ['ALL_TOUCHED=TRUE'])
Mike T.
źródło
TAK! Najwyraźniej ['ALL_TOUCHED=TRUE']niestety naprawiono tylko warstwy wielokątów. Moje punktowe warstwy plików kształtów są nadal bardzo niewyraźne i pokazują jeden piksel ponad miejscem, w którym zostały umieszczone.
Lark,
W efekcie wygląda to tak . Jest w tej samej projekcji co inne i miałem nadzieję, że to też magicznie to naprawi, ale wydaje się, że uparcie wypala jeden piksel poza miejscem, w którym się faktycznie znajduje.
Lark,
Z pewnością wygląda to na błąd, gdzie punkt zapłonu jest kompensowany przez dx / 2 i dy / 2. Zastanawiam się, czy ten błąd nadal występuje w najnowszym pniu.
Mike T
To nie! Działa w wersji 1.9.0. Dzięki wielkie!
Lark,
1
Jest też całkiem niezły przepis: gis.stackexchange.com/a/16916/9942
j08lue