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.
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 .