Использование Python для преобразования шейп-файла в Geotiff

У меня есть файл .shp, и я хотел бы преобразовать его в GEOTIFF. Мой шейп-файл состоит из большого многоугольника с множеством многоугольников внутри. Я использую следующий код, но выходной TIF состоит только из большого многоугольника.

from osgeo import ogr, gdal
import subprocess

InputVector = Shapefile
OutputImage = OutputfileName

gdalformat = 'GTiff'
datatype = gdal.GDT_Byte
burnVal = 1 
Shapefile = ogr.Open(InputVector)
Shapefile_layer = Shapefile.GetLayer()


Output = gdal.GetDriverByName(gdalformat).Create(OutputImage, RasterXSize, RasterYSize, 1, datatype, options=['COMPRESS=DEFLATE'])
Output.SetProjection(Projection)
Output.SetGeoTransform(GeoTransform) 

Band = Output.GetRasterBand(1)
Band.SetNoDataValue(0)
gdal.RasterizeLayer(Output, [1], Shapefile_layer, burn_values=[burnVal])

subprocess.call("gdaladdo --config COMPRESS_OVERVIEW DEFLATE "+OutputImage+" 2 4 8 16 32 64", shell=True)

Я не уверен, что делаю здесь неправильно.

Спасибо


person Remote Senser    schedule 04.03.2019    source источник


Ответы (1)


Мое первое предположение - убедиться, что у большего полигона есть «дыры» / NoData там, где встречаются меньшие полигоны. Я постараюсь пояснить пример из инструмента стирания ArcMap. То, что вы хотите, чтобы ваш большой многоугольник выглядел так, это то, что с синей галочкой (с очевидным наложением функции стирания), как я думаю, ваш большой многоугольник в настоящее время выглядит так, как будто ваш большой многоугольник имеет красную галочку, что приведет к тому, что ваш растр будет состоять только из большой многоугольник. Приношу свои извинения, если я полностью ошибаюсь, тогда мы найдем другое решение.

введите здесь описание изображения

В этом случае (ваш большой многоугольник похож на синюю галочку), сжигание всех многоугольников до 1 также может создать путаницу. Вы можете создать новое поле в своем шейп-файле (например, UID) и присвоить ему уникальный числовой идентификатор (он должен быть числовым). Вы можете растрировать на основе нового поля, например:

gdal.RasterizeLayer(Output, [1], Shapefile_layer, options = ['ATTRIBUTE=UID'])
person Jascha Muller    schedule 12.03.2019