使用Python GDAL小工具实现栅格拼接（替代arcgis）

栅格拼接最初使用的是arcgis，但是arcgis拼接每次只能拼接单张，拼接之后还要自动建立投影金字塔，所以如果拼很多图像并且最终的图像很大，就会非常慢（并不知道如何拼接可以不建立投影金字塔），所以就在网上搜了搜Python gdal的解决方法，大概分成两种吧。第一种是比较手动的，需要自己创建tif，指明行列分辨率和投影信息，最后将矩阵塞进去，这种方法实在是太原始了…第二种就是自动拼接，使用gdal的小工具：gdal_merge.py，下面是我学习后自己写的函数，基本上可以实现简单的arcgis中Toolbox的Mosaic To New Raster的功能了：

def m2newraster(pathout,pathin,file=0,options='COMPRESS=LZW'):
    '''
    Equal to Mosaic To New Raster in Arcgis.Using gdal_merge.py to mosaics a set of images.
    This utility will automatically mosaic a set of images. 
    All the images must be in the same coordinate system and have a matching number of bands, but they may be overlapping, and at different resolutions. 
    In areas of overlap, the last image will be copied over earlier ones.
    Ref:https://gdal.org/programs/gdal_merge.html
    
    Parameters
    ----------
    pathout : path+file name
        Output directory with file name (end with '.tif')
    pathin : path
        Path of input tif files
    file : List, optional
        The list of input tif flies. If file=0, the input file will use all tif files in 'pathin'. The default is 0.
    options : str, optional
        The options is the 'create_options' in gdal.GetDriverByName.Create(out_file, xsize, ysize, bands,band_type, create_options) 
        The default is 'COMPRESS=LZW' (Compress the output tif file using LZW method). 
    Returns
    -------
    None.
    '''
    import os
    if os.path.exists(pathin):
        if file==0:
            file_all = os.listdir(pathin)
            file = np.array([])
            for p in file_all:
                if p[-4:] == '.tif' or p[-4:] == '.TIF' :
                    file = np.append(file, p)
    else:
        raise ValueError('Path not exists')
    filein = ''
    for i in range(len(file)):
        filein = filein + file[i]+' '
    os.chdir(pathin)#change working directory so that we can find the input tif files
    path_utility = 'D:/anaconda/Scripts/gdal_merge.py'
    os.system('python {0} -co {3} -o {1} {2}'.format(path_utility,pathout,filein,options))

ps：使用前记得检查gdal_merge.py的路径，我的电脑上是在：D:/anaconda/Scripts/gdal_merge.py。
我没有找到实时打印os.system输出的方法，似乎只有在命令行才可以打印，所以这里没办法监视运行的程度了。
gdal_merge.py还有更多options，有需要去官方文档中查询。
#--------------------------------------------------方法二----------------------------------------------
写完上面这些才突然在网上发现还有一种更方便的方法，如果大家只是简单拼接的话还是直接用下面这句话好了（Warp真是yyds）：

gdal.Warp(outputfilePath,[inputrasfile1,inputrasfile2],options=options)