python gdal 栅格影像拼接工具

最新推荐文章于 2024-07-09 09:42:52 发布

Chaoying.

最新推荐文章于 2024-07-09 09:42:52 发布

阅读量3k

点赞数 4

分类专栏： GDAL 文章标签： GIS 遥感图像处理 GDAL

本文链接：https://blog.csdn.net/TJLCY/article/details/124262686

版权

GDAL 专栏收录该内容

14 篇文章 2 订阅

订阅专栏

先上代码

from osgeo import gdal, gdalconst, osr
import argparse
import os


if __name__ == "__main__":
    parser = argparse.ArgumentParser(
        prog='raster files mosaic', description='对输入目录中的栅格影像进行拼接，要求这些影像具有相同的空间参考')
    parser.add_argument('inputDir', type=str, help='待拼接影像所在的栅格目录')
    parser.add_argument('destFile', type=str, help='拼接影像输出路径')
    parser.add_argument('--resampleAlg', type=int, default=0,
                        help='重采样方法：\n{}\n{}\n{}'.format('0:NearestNeighbour', '1:Bilinear', '2:Cubic'))
    parser.add_argument('--cutlineSHP', type=str, default='',
                        required=False, help='一个SHP文件用于裁剪拼接后的栅格影像')
    args = parser.parse_args()
    inputDir = args.inputDir
    destFile = args.destFile
    print(destFile)
    resampleAlg = args.resampleAlg
    cutlineSHP = args.cutlineSHP
    tifs = os.listdir(inputDir)
    tifs = list(filter(lambda fileName: fileName.endswith('.tif'), tifs))
    if len(tifs) == 0:
        print('栅格目录为空！')
        exit(0)
    tifs = [os.path.join(inputDir, tif) for tif in tifs]
    print(tifs)
    if cutlineSHP and os.path.exists(cutlineSHP) == False:
        print("裁剪SHP文件{}不存在！".format(cutlineSHP))
        exit()
    # 检查这些待拼接影像是否句有相同的空间参考
    gdal.AllRegister()
    osrs = []
    for tif in tifs:
        ds = gdal.Open(tif, gdalconst.GA_ReadOnly)
        osr_ = gdal.Dataset.GetSpatialRef(ds)
        osrs.append(osr_)
    osr_ = osrs[0]
    for osri in osrs:
        flag = osr.SpatialReference.IsSame(osr_, osri)
        if not(flag):
            print('待拼接的栅格影像必须有相同的空间参考！')
            exit(0)
    if resampleAlg == 0:
        resampleType = gdalconst.GRA_NearestNeighbour
    elif resampleAlg == 1:
        resampleType = gdalconst.GRA_Bilinear
    else:
        resampleType = gdalconst.GRA_Cubic
    if cutlineSHP:
        options = gdal.WarpOptions(
            srcSRS=osr_, dstSRS=osr_, format='GTiff', resampleAlg=resampleType, creationOptions=["COMPRESS=LZW"], cutlineDSName=cutlineSHP, cropToCutline=True)
    else:
        options = gdal.WarpOptions(
            srcSRS=osr_, dstSRS=osr_, format='GTiff', resampleAlg=resampleType, creationOptions=["COMPRESS=LZW"])
    gdal.Warp(destFile, tifs, options=options)
    print('处理完毕！')