本文的文字及图片来源于网络,仅供学习、交流使用,不具有任何商业用途,版权归原作者所有,如有问题请及时联系我们以作处理
以下文章来源于腾讯云,作者:bugsuse。
0.前言
因为没有喝上“秋天的第一份奶茶”,准备来更新一篇推送。
在上一篇推文中,我展示了如何使用Python结合Landsat制作遥感影像图(Python干货 | 制作遥感影像图)。
对于Landsat数据来说,对某个区域的重访周期为16天,每个位置使用全球参考系(WRS)进行索引,即每一个位置都会对应一个Path和Row,相邻的影像之间会有部分区域是重叠的。
Fig.1 World Reference System
在某些遥感影像的应用场景中,如果我们关注的区域正好处于两景影像的交界处,如下图中的象山港,那我们就需要将影像拼接起来才可以使用。
单张影像是这样。
本文合并后是这样。
1.准备工作
相较于上一篇推送,我们这次为了实现遥感影像的镶嵌拼接,我们使用到了两个库, rasterio和gdal。
importrasterio as rioimport gdal
先介绍一下我们实现两组遥感影像拼接的思路,首先选取两景相邻的影像,分别得到他们的空间范围,再得到两景组合到一起之后的空间范围,使用gdal新建一个tif文件(数据中转用),分别得到原来两景影像在新建的tif文件中的起始位置,将对应的数据写入新的tif文件中,即实现镶嵌拼接。
上面说的是两景影像的拼接,如果是更多影像拼接同样适用,但是现阶段的方法如果拼接多的影像的话,需要的内存空间很大,容易导致内存溢出,感兴趣的朋友可以思考一下如何高效实现多景影像的拼接。
其中还有两处关键处理,一是如何去除重叠区域的无效信息,二是重叠区域的数据如何选择。希望各位看官能从代码里面找到答案。
2.动起手来
得到输入影像的四个角点。
deftiffileList2filename(tiffileList):
filename=[]
prefix=[]for ifile intiffileList:
file0= ifile.split("\\")[-1]
prefix.append(os.path.join(ifile, file0))
filename.append(os.path.join(ifile, file0)+ "_B1.TIF")returnfilename, prefixdefget_extent(tiffileList):
filename, prefix=tiffileList2filename(tiffileList)
rioData=rio.open(filename[0])
left=rioData.bounds[0]
bottom= rioData.bounds[1]
right= rioData.bounds[2]
top= rioData.bounds[3]for ifile in filename[1:]:
rioData=rio.open(ifile)
left=min(left, rioData.bounds[0])
bottom= min(bottom, rioData.bounds[1])
right= max(right, rioData.bounds[2])
top= max(top, rioData.bounds[3])return left, bottom, right, top, filename, prefix
得到新建tif文件的size,这里已知Landsat空间分辨率为30m,如果是其他遥感数据,需对应进行修改。
```pythondefgetRowCol(left, bottom, right, top):
cols= int((right - left) / 30.0)
rows= int((top - bottom) / 30.0)return cols, rows
主程序,其中plot_rgb为上一篇推送中用到的函数。
if __name__ == '__main__':
tiffileList= [r'PathofLandsat8\LC08_L1TP_118039_20160126_20170330_01_T1',
r'PathofLandsat8\LC08_L1TP_118040_20160126_20170330_01_T1']
left, bottom, right, top, filename, prefix=get_extent(tiffileList)
cols, rows=getRowCol(left, bottom, right, top)
bands= ['B7', 'B5', 'B3']
n_bands=len(bands)
arr= np.zeros((n_bands, rows, cols), dtype=np.float)#打开一个tif文件
in_ds =gdal.Open(filename[0])for i inrange(len(bands)):
ibands=bands[i]#新建一个tif文件
driver = gdal.GetDriverByName('gtiff')
out_ds= driver.Create(ibands + 'mosaic.tif', cols, rows)#设置tif文件的投影
out_ds.SetProjection(in_ds.GetProjection())
out_band= out_ds.GetRasterBand(1)#设置新tif文件的地理变换
gt =list(in_ds.GetGeoTransform())
gt[0], gt[3] =left, top
out_ds.SetGeoTransform(gt)#对要拼接的影像进行循环读取
for ifile inprefix:
in_ds= gdal.Open(ifile + '_' + ibands + '.tif')#计算新建的tif文件及本次打开的tif文件之间的坐标漂移
trans =gdal.Transformer(in_ds, out_ds, [])#得到偏移起始点
success, xyz =trans.TransformPoint(False, 0, 0)
x, y, z=map(int, xyz)#读取波段信息
fnBand = in_ds.GetRasterBand(1)
data=fnBand.ReadAsArray()#写入tif文件之前,最大值设置为255,这一步很关键
data = data / 65535 * 255data[np.where(data== 255)] =0#影像重合部分处理,重合部分取最大值
xSize =fnBand.XSize
ySize=fnBand.YSize
outData=out_band.ReadAsArray(x, y, xSize, ySize)
data=np.maximum(data, outData)
out_band.WriteArray(data, x, y)delout_band, out_ds
file2read= ibands + 'mosaic.tif'arr[i, :, :]=tiff.imread(file2read)
os.remove(file2read)
plot_rgb(arr, rgb=(0, 1, 2))
3.小结
大功告成