Блог

У меня есть шейп-файл береговой линии пяти великих озер и привязанная к географической привязке фотография озера Онтарио. Я хочу наложить шейп-файл береговой линии на фотографию tiff в качестве нового растрового слоя. В частности, я хочу сопоставить пиксели шейп-файла и фотографии tiff с одинаковыми координатами. Результатом будет растровый слой m x n (размер одного канала фотографии tiff), где значения равны 1, если координата совпадает с файлом tiff, и 0, если не совпадает. Что мне делать после прочтения этих файлов, как показано ниже? Я могу использовать базовые gdal, растерио и геопанды в Python. Я искал существующие ответы в течение недели, но не смог найти ни одного.

 # read tiff file
ds = gdal.Open(path)
band = ds.GetRasterBand(1)
array = band.ReadAsArray()
# read shapefile
shoreline_delineation = gpd.read_file(shoreline_path)

 

Ниже приведены два изображения, о которых я упоминал. Большое спасибо!
Файлы Geotiff шейп
-файлов

Добавление растеризованной версии шейп-файла можно выполнить с помощью кода, приведенного ниже.

Предположим, у вас есть два входных файла, растр и вектор, и предполагается, что ваш вектор уже имеет тип linestring (не полигоны).

 from osgeo import gdal

shp_file = 'D:/Temp/input_vector.shp'
ras_file = 'D:/Temp/input_raster.tif'
ras_file_overlay = 'D:/Temp/output_raster_overlay.tif'
 

Прочитайте свойства входного растра:

 ds = gdal.Open(ras_file)
gt = ds.GetGeoTransform()
proj = ds.GetProjection()
n_bands = ds.RasterCount
xsize = ds.RasterXSize
ysize = ds.RasterYSize
ds = None
 

Вычислите экстент растра:

 ulx, xres, _, uly, _, yres = gt

lrx = ulx   xres * xsize
lry = uly   yres * ysize
 

Растрируйте вектор в геотиф в памяти и считайте его как массив Numpy:

 opts = gdal.RasterizeOptions(
    outputType=gdal.GDT_Byte,
    outputBounds=[ulx, lry, lrx, uly],
    outputSRS=proj,
    xRes=xres,
    yRes=yres,
    allTouched=False,
    initValues=0,
    burnValues=1,
)

tmp_file = '/vsimem/tmp'
ds_overlay = gdal.Rasterize(tmp_file, shp_file, options=opts)
overlay = ds_overlay.ReadAsArray()
ds_overlay = None
gdal.Unlink(tmp_file)
 

Вы можете изменить tmp_file местоположение на диске, чтобы записать его в качестве промежуточного вывода (однослойный). Это может быть удобно для отладки.

При необходимости вы можете добавить буфер к береговой линии, чтобы создать большую/более широкую границу на выходе, добавив что-то вроде этого к gdal.RasterizeOptions вышесказанному:

 SQLStatement="SELECT ST_Buffer(geometry, 0.1) FROM <layer_name>",
SQLDialect="sqlite",
 

Драйвер Geotiff не поддерживает метод «AddBand», поэтому вы не можете напрямую добавить его во входной растр. Но добавить его в копию можно с помощью:

 ds = gdal.Open(ras_file)

tmp_ds = gdal.GetDriverByName('MEM').CreateCopy('', ds, 0)
tmp_ds.AddBand()
tmp_ds.GetRasterBand(tmp_ds.RasterCount).WriteArray(overlay)

dst_ds = gdal.GetDriverByName('GTIFF').CreateCopy(ras_file_overlay, tmp_ds, 0)
dst_ds = None
ds = None
 

В качестве альтернативы вы можете записать растрированный слой непосредственно на диск в формате tiff, а затем использовать что-то вроде gdal.BuildVRT их объединения.