天下无云第二步，逐像素图片合成

在做出第一步除云之后，停滞了很长时间，一方面是因为一直在尝试摸索其余软件比如envi，snap，arcmap等等，但是最终发现都不尽人意，且不说前两个软件根本没有这样功能，第三个位移能实现图片逻辑向加减的也只能对类似于颜色均匀图片实现，根们不可能按照我的想法——按照像素合成
所以，靠天靠地不如靠自己啊，在这里先放所有代码，虽然无云的图片很多网络地图公司都能做，但是没有人公开过技术细节，但是我希望技术一直是开放的，不避让后人重复造车轮，共同进步：

import os
import os.path
from osgeo import gdal
import numpy as np
def readpath(bigpath):
    # 读取文件列表
    # 这里通过筛选掉了enp文件，留下的只有光tif文件，把含有enp的单独储存然后一块删掉
    files = os.listdir(bigpath)
    M=[]
    for file in files:
        if 'enp' in file:
            M.append(file)
        elif 'hdr' in file:
            M.append(file)
    for m in M:
        files.remove(m)
    return files


def read_tiff(inpath):
    ds = gdal.Open(inpath)
    row = ds.RasterYSize
    col = ds.RasterXSize
    band = ds.RasterCount
    geoTransform = ds.GetGeoTransform()
    proj = ds.GetProjection()
    data = np.zeros([row, col, band],np.int16) #这里注意务必要改成uint16不然没法显示图像
    for i in range(1, band+1):#因为波段是从1开始的所以这里要从1开始，其实所能读取的文件远不止tiff，很多遥感格式都没问题
        #而且range（1，n）是从1到n-1的list所以要band+1
        dt = ds.GetRasterBand(i)
        data[:, :, i-1] = dt.ReadAsArray(0, 0, col, row)
    return data, geoTransform, proj
def write_tiff(outpath,data,geoTransform,proj):
    if 'int8' in data.dtype.name:
        datatype = gdal.GDT_Byte
    elif 'int16' in data.dtype.name:
        datatype = gdal.GDT_UInt16
    else:
        datatype = gdal.GDT_Float32
    #先搞清楚读进来的数据的格式,决定输出也使用相同的格式,实际验证刚好不损失精度的大小的类型就是int16了
    if len(data.shape) == 3:
        im_height, im_width,im_bands = data.shape
    elif len(data.shape) == 2:
        im_height, im_width = data.shape
    else:
        im_bands, (im_height, im_width) = 1,data.shape
    #搞清楚矩阵的波段以及高度宽度
    driver = gdal.GetDriverByName('Gtiff')
    dataset = driver.Create(outpath, im_width, im_height, im_bands, datatype)
    if(dataset!= None):
        dataset.SetGeoTransform(geoTransform) #写入仿射变换参数
        dataset.SetProjection(proj) #写入投影
    for i in range(im_bands):
        dataset.GetRasterBand(i+1).WriteArray(data[:,:,i])#这里对应的是readtiff里面的从0开始存储的矩阵数据
    del dataset


def combine(bigpath, fname):
    # 这里得考虑进去每个波段
    cpname = [0]*len(fname)
    for i in range(len(fname)):
        cpname[i] = bigpath +'\\' + fname