如何将海岸线 shapefile 叠加到地理 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)

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

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 更改为磁盘上的某个位置，以将其写入为中间输出（单层）。这可以方便调试。

如果需要，您可以向海岸线添加缓冲区以在输出中创建 bigger/wider 边界，方法是向上面的 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 之类的东西将它们堆叠在一起。

如何将海岸线 shapefile 叠加到地理 tiff 文件中？

How to overlay shoreline shapefile into geo tiff file?

python

gis

shapefile

gdal