如何将海岸线图形文件叠加成地理图形文件？
如何将海岸线图形文件叠加成地理图形文件？
提问于 2021-09-12 17:09:37
我有一个五大湖的海岸线形状文件和一张安大略湖的地理参考湖岸照片。我想把海岸线的shapefile叠加成一张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
回答 1
Stack Overflow用户
回答已采纳
发布于 2021-09-14 07:11:29
添加一个光栅化版本的Shapefile可以用下面的代码来完成。
假设你有两个输入文件，一个栅格和一个矢量，并且假设你的矢量已经是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将向量栅格化为内存中的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之类的工具将它们堆叠在一起。
页面原文内容由Stack Overflow提供。腾讯云小微IT领域专用引擎提供翻译支持
原文链接：
https://stackoverflow.com/questions/69153624
复制相关文章
相似问题
腾讯云开发者
