Я пытаюсь записать шейп-файл в растр, используя RasterizeLayer из GDAL. Я предварительно создаю растр области интереса из другого шейп-файла, учитывая определенный размер пикселя. Этот AOI затем служит базой для всех последующих растеризаций (одинаковое количество столбцов и строк, одинаковая проекция и геотрансформация).
Проблема возникает, однако, когда я иду, чтобы записать формы в их собственный растр, основываясь на том же размере пикселей и проекций. Ссылка ниже (не хватает представителя, чтобы опубликовать изображение), показывает оригинальный шейп-файл в tan и темно-розовый, где RasterizeLayer записал данные. Светло-розовый - это значения узлов для темно-розовых растровых данных. Серый цвет - это AOI, на основе которого был завершен процесс записи шейп-файла.
Учитывая размеры полигонов шейп-файлов, я ожидаю увидеть значения записи в двух нижних углах, а также в двух пикселях под данными, которые отображаются. Очевидно, однако, это не тот случай.
Ниже приведен код, который я использовал для их создания. Все фигуры были созданы с использованием QGIS, и все они были созданы в одной проекции. (Следует отметить, что сетка на показанном рисунке была просто для того, чтобы дать представление о размере пикселя, который я использовал.)
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])
Это ошибка в GDAL или RasterizeLayer записывает данные на основе чего-то другого, кроме наличия или отсутствия многоугольника в пределах определенной области пикселей?
Файлы, которые я использовал, можно найти здесь .
источник
Ответы:
Я играл с GDALRasterizeLayers на этой неделе и довольно хорошо представляю, что он делает. По умолчанию он растеризует пиксель, если центр пикселя находится внутри многоугольника. Если в центре ничего нет, оно не будет растеризовано, даже если в пределах пикселя есть части многоугольника. Чтобы растеризация работала так, как вы хотите, попробуйте опцию «ALL_TOUCHED»:
источник
['ALL_TOUCHED=TRUE']
, хотя, к сожалению, это только фиксированные полигональные слои. Слои моего точечного шейп-файла все еще супер-шаткие, и отображаются на один пиксель от того места, где они размещены.