【GDAL】tif影像拼接和目标截取

原文作者：我辈李想
版权声明：文章原创，转载时请务必加上原文超链接、作者信息和本声明。

---

一、gdal.Warp拼接tif

from osgeo import gdal, gdalconstminX,minY,maxX,maxY = latmin_input,lonmin_input, latmax_input,lonmax_inputprint(minX,minY,maxX,maxY)
tifPaths = [path1,path2]
inputFiles = []
outputfilePath = './wrap.tif'for path in tifPaths:inputrasfile = gdal.Open(path, gdal.GA_ReadOnly) # 读取影像inputProj = inputrasfile.GetProjection() # 获取坐标系inputFiles.append(inputrasfile) # 推入列表options = gdal.WarpOptions(srcSRS=inputProj, # 输入坐标系dstSRS=inputProj,  # 输出坐标系format='GTiff',  # 图像格式resampleAlg=gdalconst.GRA_Bilinear, # 重采样算法，这里是双线性内插dstNodata=-1, # 缺省值outputBounds =(minX,minY,maxX,maxY),# cutlineLayer=outline, # 输出范围，这里可以是一个外轮廓shp数据cropToCutline=True, # 是否对输出边界使用剪切线范围，即outputBounds或cutlineLayeroutputType=gdalconst.GDT_Int32) # 数据类型，这里是有符号32位整型gdal.Warp(outputfilePath,inputFiles,options=options) # 图像镶嵌out_my = gdal.Open(outputfilePath, gdal.GA_ReadOnly).ReadAsArray()
plt.figure(figsize=(7,7))
plt.imshow(out_my,'gray')

二、截取

outputBounds 和cutlineLayer二选一，cropToCutline参数确定是否按照outputBounds 和cutlineLayer输出tif。

1.通过经纬范围截取拼接的影像

outputBounds 输入经纬范围的元组，经纬应该与拼接的inputFiles坐标系一致

2.通过shp范围截取凭借后影像

cutlineLayer输入shp文件的路径，经纬应该与拼接的inputFiles坐标系一致

三、WarpOptions其他参数

创建一个可以传递给gdal.Warp（）的WarpOptions（）对象关键字参数为：options---可以是一个字符串数组、一个字符串，也可以由其他关键字清空和填充。format---输出格式（“GTiff”等）outputBounds——目标SRS中的输出边界为（minX，minY，maxX，maxY）outputBoundsSRS——如果输出边界未在dstSRS中表示，则表示输出边界的SRSxRes，yRes——目标SRS中的输出分辨率targetAlignedPixels——是否强制输出边界为输出分辨率的倍数width——输出光栅的宽度（以像素为单位）height—输出光栅的高度（以像素为单位）srcSRS---源SRSdstSRS--输出SRScoordinateOperation—作为PROJ字符串或WKT字符串的坐标操作srcAlpha——是否强制将输入数据集的最后一个波段视为阿尔法波段dstAlpha——是否强制创建输出alpha波段outputType—输出类型（gdalcons.GDT_Byte等）workingType—工作类型（gdalcons.GDT_Byte等）warpOptions--扭曲选项列表errorThreshold—近似变换器的错误阈值（以像素为单位）warpMemoryLimit—工作缓冲区的大小（MB）resampleAlg-重新采样模式creationOptions--创建选项列表srcNodata--源nodata值dstNodata—输出nodata值多线程—是否执行多线程计算和I/O操作tps--是否使用薄板花键GCP变压器rpc——是否使用rpc转换器geoloc——是否使用GeoLocation阵列转换器多项式阶数——多项式GCP插值的阶数transformerOptions-transformer选项列表cutlineDSName--剪切线数据集名称cutlineLayer—剪切线图层名称cutlineWhere--cutlineWhere子句cutlineSQL—cutline SQL语句cutlineBlend—以像素为单位的剪切线混合距离cropToCutline—是否对输出边界使用剪切线范围copyMetadata—是否复制源元数据metadataConflictValue--元数据数据冲突值setColorInterpretation—是否将输入波段的颜色解释强制到输出波段overviewLevel—指定必须使用的源文件的概述级别callback—回调方法callback_data——回调的用户数据

四、其他方式裁剪

1.通过shapely的Polygon裁剪

from shapely.geometry import Point, Polygonpoints = [(latmin_input,lonmin_input), (latmax_input,lonmin_input), (latmax_input,lonmax_input), (latmin_input,lonmax_input)]
polygon = Polygon(points)import rasterio as rio
from rasterio.mask import mask
with rio.open(rasterPath1) as rasterdata1:out_crs1 = rasterdata1.crsfeature1 = [polygon] # 1.3中得到的polygonout_image1, out_transform1 = mask(rasterdata1, feature1, all_touched=True, crop=True, nodata=0)