本文共 8649 字，大约阅读时间需要 28 分钟。

在做出第一步除云之后，停滞了很长时间，一方面是因为一直在尝试摸索其余软件比如envi，snap，arcmap等等，但是最终发现都不尽人意，且不说前两个软件根本没有这样功能，第三个位移能实现图片逻辑向加减的也只能对类似于颜色均匀图片实现，根们不可能按照我的想法——按照像素合成

所以，靠天靠地不如靠自己啊，在这里先放所有代码，虽然无云的图片很多网络地图公司都能做，但是没有人公开过技术细节，但是我希望技术一直是开放的，不避让后人重复造车轮，共同进步：

import osimport os.pathfrom osgeo import gdalimport numpy as npdef 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 filesdef 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, projdef 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 datasetdef combine(bigpath, fname):    # 这里得考虑进去每个波段    cpname = [0]*len(fname)    for i in range(len(fname)):        cpname[i] = bigpath +'\\' + fname[i]  # 给出了完整的路径    S = [0] * len(cpname) * 3  # 三给为一组，分别可以代表数据，地理转换，投影坐标    for i in range(len(cpname)):        S[3 * i:3 * i + 3] = read_tiff(cpname[i])        # 读取出了所有的数据，3i表示第i个数据，3i+1是他的转换，+2是投影信息    Cgeo = [0] * 6    for i in range(len(cpname)):        if i == 0:            Cgeo = S[3 * i + 1]        else:            if Cgeo[0] > S[3 * i + 1][0]:                Cgeo[0] = S[3 * i + 1][0]            if Cgeo[3] < S[3 * i + 1][3]:                Cgeo[3] = S[3 * i + 1][3]    # 这样一来就得到整个大图的最左上角坐标，也就相当于输出图像转换坐标，下面来求右下角    Rloco = [[0, 0]] * len(cpname)  # 这个用来存储每个图片的右下角坐标    for i in range(len(cpname)):        d_y, d_x, _ = S[3 * i].shape#返回格式是行数，列数，维度数        X = S[3 * i + 1][0] + d_x * S[3 * i + 1][1]        Y = S[3 * i + 1][3] + d_y * S[3 * i + 1][5]        Rloco[i] = [X, Y]    # 接下来需要得到图片最右下角的坐标    RXY = [0, 0]  # 存储最右下角坐标    for i in range(len(cpname)):        if i == 0:            RXY = Rloco[i]        if RXY[0] < Rloco[i][0]:            RXY[0] = Rloco[i][0]        if RXY[1] > Rloco[i][1]:            RXY[1] = Rloco[i][1]    # 至此得到右下角坐标    _, _, band = S[0].shape    Xsize = int(np.abs((Cgeo[0] - RXY[0]) / Cgeo[1]))    Ysize = int(np.abs((Cgeo[3] - RXY[1]) / Cgeo[5]))    Cd = np.zeros((Ysize, Xsize,  band), np.int16)    Cproj = S[2]    for i in range(len(cpname)):        a, b, c = S[3 * i].shape        deltx = int(np.abs((Cgeo[0] - S[3 * i + 1][0]) / Cgeo[1]))        delty = int(np.abs((Cgeo[3] - S[3 * i + 1][3]) / Cgeo[5]))        if i == 0:            Cd[delty :delty  + a, deltx :deltx + b, :] = S[3 * i]        else:            for j in range(c):#波段                for k in range(a):#行                    for m in range(b):#列                        pix = Cd[k + delty , m + deltx , j]                        if pix == 0:                            pix = S[3 * i][k, m, j]                        Cd[delty  + k, deltx  + m, j] = pix    return Cd, Cgeo, Cprojif __name__ == '__main__':    path='E:\yi_data\Processed_shanghai\hecheng'    filename=readpath(path)    A,B,C=combine(path,filename)    write_tiff('E:\yi_data\Processed_shanghai\hecheng\\big.tif',A,B,C)

接下来开始讲解原理：

首先，在理解像素镶嵌之前，必须要了解怎么读取信息，怎么把读取到的信息提取出来，怎么处理，怎么写出。 好，问题就来了，第一步，怎么读数据。关于数据源，我使用的是哨兵二号的数据，相关下载方式网络有很多就不赘述了。我们先看看这里面数据是什么样的吧 打开这些裸数据之后，得到的是这些 打开E:\yi_data\shanghai\S2A_MSIL1C_20181018T022701_N0206_R046_T51RUQ_20181018T053000.SAFE\GRANULE\L1C_T51RUQ_A017346_20181018T023228 这个链接以后，你会看到这么个 最这个数据的软件当然就是snap了，同样可以在官网下，直接打开MTD_TL.xml就可以了，但是其实这个软件功能不是很强大所以我需要envi处理呀，当然了如果你使用的是envi5.5这些新版本可以无视，所以首先我们需要先把数据重采样输出为envi格式，这些具体操作在上一篇已经讲过了，不再赘述。 当你得到一个遥感图像时，云层的出现位置是不确定的，他可大可小，也可以出现在任何一个区域，因此常规的大范围拼接是绝对无法满足需求的。 需要一个逐个像素的高精度拼接，而且必须对准所有的坐标，与此同时还要有波段的匹配。 那么对于图像，处理速度最快的方式就是使用数组了，我们第一步是要把图像放在一个数组里面，维度根据长宽以及波段决定。 可是，最重要的一个问题就出现了，你怎么把图像上的某个点映射到数组里面，换言之如果仅仅读取出来不清楚位置是无法拼接的。因此重点之一就是确定每个像素点的坐标。坐标确定好之后，再把图像按照位置拼接好，然后再输出为一个大图像。 好了既然大方向搞清楚了，那么就要关注一些具体细节了，第一点你怎么知道哪个点有云没云的？这个在第一篇文章里已经把所有的判断为有云的点赋值为0了，所以只需要检测有没有0，就知道是不是云点，就可以知道把不把另一幅图可能的无云点加上去。

我们可以把功能分解为3个部分，读，处理，写三大部分。

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

这是用来读数数据的，通过一个大路径读取出所有内部文件名，然后，因为envi格式有头文件，同时还可能有你使用envi软件处理时候产生的金字塔文件，所以我就检测一遍文件名，凡是含有这两个后缀的文件都被提出来，最后统一去除

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

这是读取文件的函数，里面的变量都名字比较清晰就不做多的解释了，需要强调的一点是，里面的X是向右为正方向，Y向上为正方向，投影变换坐标是至关重要的，仿射矩阵信息有六个参数，描述的是栅格行列号和地理坐标之间的关系：

‘’’ 0：图像左上角的X坐标； 1：图像东西方向分辨率； 2：旋转角度，如果图像北方朝上，该值为0； 3：图像左上角的Y坐标； 4：旋转角度，如果图像北方朝上，该值为0； 5：图像南北方向分辨率； ‘’’ 所以根据这些数据就可以得到每个点的绝对像素坐标 得到这六个参数之后就可以进行图像像素坐标（即行列号）和地理坐标之间的变换：

XGeo = GeoTransform[0]+GeoTransform[1]*Xpixel+Ypixel*GeoTransform[2];YGeo = GeoTransform[3]+GeoTransform[5]*Ypixel+Ypixel*GeoTransform[5];

不过现在大家所下载到的绝大多数遥感影像，都是旋转角度为0，面向正北方的，其他地方用0填补了而已。所以这就给了我们很大的方便，不用在单独存取每个点的坐标，否则可能还需要用字典来存储这些信息。

下面讲解处理函数：

cpname = [0]*len(fname)for i in range(len(fname)):    cpname[i] = bigpath +'\\' + fname[i]  # 给出了完整的路径S = [0] * len(cpname) * 3  # 三给为一组，分别可以代表数据，地理转换，投影坐标

用cpname存储了完整的路径，同时预备好list存每张图的三份数据。

Cgeo = [0] * 6    for i in range(len(cpname)):        if i == 0:            Cgeo = S[3 * i + 1]        else:            if Cgeo[0] > S[3 * i + 1][0]:                Cgeo[0] = S[3 * i + 1][0]            if Cgeo[3] < S[3 * i + 1][3]:                Cgeo[3] = S[3 * i + 1][3]

Cgeo是存储输出的图像的转换坐标的，这里检测所有图像里最偏坐上方的点，以它为我们的输出坐标，而旋角之类是0，还有一点要注意的是geotransform的0,5位置存的是X，Y方向的分辨率，注意他是带符号的，正代表正方向，所以你可以看到Y方向是负数就是说明图像是向下发展的

Rloco = [[0, 0]] * len(cpname)  # 这个用来存储每个图片的右下角坐标    for i in range(len(cpname)):        d_y, d_x, _ = S[3 * i].shape#返回格式是行数，列数，维度数        X = S[3 * i + 1][0] + d_x * S[3 * i + 1][1]        Y = S[3 * i + 1][3] + d_y * S[3 * i + 1][5]        Rloco[i] = [X, Y]    # 接下来需要得到图片最右下角的坐标

因为要一片一片的填补数据，所以需要把每个图的右下角坐标得到

RXY = [0, 0]  # 存储最右下角坐标    for i in range(len(cpname)):        if i == 0:            RXY = Rloco[i]        if RXY[0] < Rloco[i][0]:            RXY[0] = Rloco[i][0]        if RXY[1] > Rloco[i][1]:            RXY[1] = Rloco[i][1]    # 至此得到右下角坐标

这里是存放整个图片右下角坐标，用来确定图像大小的

_, _, band = S[0].shape    Xsize = int(np.abs((Cgeo[0] - RXY[0]) / Cgeo[1]))    Ysize = int(np.abs((Cgeo[3] - RXY[1]) / Cgeo[5]))    Cd = np.zeros((Ysize, Xsize,  band), np.int16)    Cproj = S[2]

这里用band得到波段数，同时决定xy方向大小

for i in range(len(cpname)):        a, b, c = S[3 * i].shape        deltx = int(np.abs((Cgeo[0] - S[3 * i + 1][0]) / Cgeo[1]))        delty = int(np.abs((Cgeo[3] - S[3 * i + 1][3]) / Cgeo[5]))        if i == 0:            Cd[delty :delty  + a, deltx :deltx + b, :] = S[3 * i]        else:            for j in range(c):#波段                for k in range(a):#行                    for m in range(b):#列                        pix = Cd[k + delty , m + deltx , j]                        if pix == 0:                            pix = S[3 * i][k, m, j]                        Cd[delty  + k, deltx  + m, j] = pix    del S    Cmdata=gdal.Open(cpname[0]).GetMetadata()

这是代码最花时间的部分，三个维度的检测降决定是否填补，最后的一行是得到头文件信息，最后用来输出文件补上头文件的

好了说完了接下来自己运行一下就可以得到结果了，速度比较慢是很正常的

转载地址：http://kmhrn.baihongyu.com/