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

在做出第一步除云之后，停滞了很长时间，一方面是因为一直在尝试摸索其余软件比如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 filesdef read_tiff(inpath):ds = gdal.Open(inpath)row = ds.RasterYSizecol = ds.RasterXSizeband = ds.RasterCountgeoTransform = 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+1dt = 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_Byteelif 'int16' in data.dtype.name:datatype = gdal.GDT_UInt16else:datatype = gdal.GDT_Float32#先搞清楚读进来的数据的格式,决定输出也使用相同的格式,实际验证刚好不损失精度的大小的类型就是int16了if len(data.shape) == 3:im_height, im_width,im_bands = data.shapeelif len(data.shape) == 2:im_height, im_width = data.shapeelse: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] * 6for 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].shapeXsize = 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].shapedeltx = 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] = pixreturn 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.RasterYSizecol = ds.RasterXSizeband = ds.RasterCountgeoTransform = 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+1dt = 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] * 6for 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].shapeXsize = 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].shapedeltx = 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] = pixdel SCmdata=gdal.Open(cpname[0]).GetMetadata()

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

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

天下无云第二步，逐像素图片合成相关推荐

php合并多张gif图,多张静态图片合成一张动态图-静图合成动图制作
现如今制作gif图片也不是什么难事了,巧用动态图合成软件,即能轻松将多张静态图片合成一张动态图.动态图片是由多张不同的静态图片组合而成的gif格式图片,它会按照一定的顺序和时间进行逐帧播放.做好的动态 ...
图片合成gif_谈谈有哪些好用的制作GIF的方式
今天咱们谈谈计算机上几种制作 GIF 的方式,而且他们都是免费的,但可能需要你会一点计算机操作技能.本文会简单的介绍一些常用的GIF软件,比较一下功能和上手程度.常见的GIF录制方式大概有这么几种方式 ...
MATLAB小白之图片合成
MATLAB之图片合成将两张图片进行合成,先上代码: img = imread('b1.png'); %R通道 R = double(img(:,:,1)); %G通道 G = double(img ...
视频和图片合成软件，简单快速合成视频和图片
怎么把小视频和图片合成起来?有简单好上手的教程吗?今天就教大家简单几步,把视频和图片合成为照片视频.先看看用数码大师合成视频和图片的效果截图: 第一步:把图片一次性导入,为照片配上文字点击" ...
百度飞桨AI抠图+图片合成
考虑到本科学校校庆即将到来,而又刚好学习了百度飞桨AI抠图以及图片合成的相关课程,因而想合成一张自己和本科学校的合照.(由于才疏学浅,略有翻车,请见谅) 使用工具:百度PaddleHub DeepLa ...
使用canvas在前端实现图片合成
看着总结的不错,我也就拿来主义了,做个记录,侵权必删图片合成最常见的需求有验证码图片,亦或者图片加水印等,这种实现一般都是后端实现的. 随着HTML5发展和现代浏览器的占比越来越高,我们其实也可以在 ...
python合成图片_python图片合成的示例
python的PIL库简直好用的不得了,PIL下面的Image库更是封装了很多对图片处理的函数,关于Image库的介绍和使用,看这里:http://effbot.org/imagingbook/ima ...
前端生成二维码与图片合成
首先前端生成二维码使用插件完成,插件为DrawQRCode.unitypackage 合成二维码直接使用DrawQRCode 类里提供的方法即可生成二维码的方法 DrawCode_Color32 ...
c语言.jpg图片转成数组_多张jpg图片合成pdf文件
唐县职称公众号微信视频号评审条件(小程序) 关于如何将多张jpg图片合成pdf文件首先再强调一下"扫描",不要用手机拍照.高拍仪代替扫描.否则责任自负. 如何将多张JPG图片 ...

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

天下无云第二步，逐像素图片合成相关推荐

最新文章

热门文章