GDAL RasterizeLayer não grava todos os polígonos em raster?

12

Estou tentando gravar um shapefile em uma varredura usando o RasterizeLayer da GDAL. Eu pré-crio uma área de interesse de varredura a partir de um shapefile diferente, considerando um tamanho de pixel específico. Essa AOI serve como base para todas as seguintes rasterizações (mesmo número de colunas e linhas, mesma projeção e geotransformação).

O problema ocorre, no entanto, quando vou gravar as formas em sua própria varredura, com base no mesmo tamanho de pixel e projeções. O link abaixo (não possui representante suficiente para postar a imagem) mostra o shapefile original em marrom e o rosa escuro em que o RasterizeLayer gravou dados. Os valores rosa claro são nodata para os dados raster rosa escuro. O cinza é a AOI com base na qual a gravação do shapefile foi concluída.

Dadas as extensões dos polígonos do shapefile, esperaria ver os valores de gravação nos dois cantos inferiores, bem como os dois pixels abaixo dos dados mostrados. Obviamente, no entanto, este não é o caso.

Imagem de Problema - Queimadas de varredura concluídas

A seguir, é o código que eu usei para gerá-los. Todas as formas foram criadas usando o QGIS e todas criadas na mesma projeção. (Observe que a grade na imagem mostrada era apenas para dar uma idéia do tamanho do pixel que eu estava 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])

Isso é um bug no GDAL ou o RasterizeLayer está gravando dados com base em algo que não seja apenas a presença ou a falta de um polígono dentro de uma área de pixel especificada?

Os arquivos que eu estava usando podem ser encontrados aqui .

Cotovia
fonte
Você pode fornecer um link para 'activity_3.shp' e 'AOI_Raster.tif'? Quero ver se consigo recriar do meu lado.
Ricos

Respostas:

10

Estive jogando com GDALRasterizeLayers esta semana e tenho uma boa ideia do que está fazendo. Por padrão, ele rasterizará um pixel se o centro do pixel estiver dentro do polígono. Se não houver nada no centro, ele não será rasterizado, mesmo se houver partes de um polígono dentro dos limites de pixel. Para permitir que a rasterização funcione da maneira que você pretende, tente a opção "ALL_TOUCHED":

gdal.RasterizeLayer(raster_dataset, [1], shape_layer, None, None, [1], ['ALL_TOUCHED=TRUE'])
Mike T
fonte
SIM! Aparentemente ['ALL_TOUCHED=TRUE'], embora infelizmente, apenas as camadas poligonais fixas. Minhas camadas de shapefile de ponto ainda estão sendo super instáveis ​​e aparecem um pixel acima de onde elas estão colocadas.
Cotovia
Ele acaba parecendo esta . Está na mesma projeção que as outras, e eu esperava que isso também pudesse corrigi-lo magicamente, mas parece queimar teimosamente um pixel de onde ele está realmente localizado.
Cotovia
Isso certamente parece digno de bug, onde o ponto de gravação é compensado por dx / 2 e dy / 2. Gostaria de saber se esse bug ainda persiste com o último tronco.
Mike T
Isso não! Funciona em 1.9.0. Muito obrigado!
Cotovia
1
Há também é bastante uma boa receita aqui: gis.stackexchange.com/a/16916/9942
j08lue