如何将海岸线shape文件叠加到geo tiff文件中？

Question

我有五个大湖的海岸线形状文件和安大略湖的地理参考湖岸 tiff 照片。我想将海岸线形状文件作为新的栅格图层叠加到 tiff 照片中。具体来说，我想将 shapefile 和 tiff 照片的像素与相同的坐标进行匹配。结果将是一个 m x n（tiff 照片一个通道的尺寸）栅格图层，其中如果坐标与 tiff 文件匹配，则值为 1，如果不匹配，则值为 0。阅读以下这些文件后我应该做什么？我可以在 Python 中使用基本的 gdal、rasterio 和 geopandas。我搜索现有答案一周但找不到。

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

下面是我提到的两张图片。非常感谢！

Answer 1

添加 Shapefile 的光栅化版本可以使用下面的代码来完成。

假设您有两个输入文件，一个栅格和一个向量，并且假设您的向量已经是线串类型（不是多边形）。

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

将向量光栅化为内存中的 Geotiff，并将其读取为 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

的内容将它们堆叠在一起。

Answer 2

我还有一个问题。如果我不想添加额外的band，而是想替换向量（shp）对应位置的像素值（tiff），我应该如何实现