本文的文字及圖片來(lái)源于網(wǎng)絡(luò),僅供學(xué)習(xí)、交流使用,不具有任何商業(yè)用途,版權(quán)歸原作者所有,如有問(wèn)題請(qǐng)及時(shí)聯(lián)系我們以作處理

以下文章來(lái)源于騰訊云，作者：bugsuse。

0.前言

因?yàn)闆](méi)有喝上“秋天的第一份奶茶”，準(zhǔn)備來(lái)更新一篇推送。

在上一篇推文中，我展示了如何使用Python結(jié)合Landsat制作遙感影像圖(Python干貨 | 制作遙感影像圖)。

對(duì)于Landsat數(shù)據(jù)來(lái)說(shuō)，對(duì)某個(gè)區(qū)域的重訪周期為16天，每個(gè)位置使用全球參考系(WRS)進(jìn)行索引，即每一個(gè)位置都會(huì)對(duì)應(yīng)一個(gè)Path和Row，相鄰的影像之間會(huì)有部分區(qū)域是重疊的。

Fig.1 World Reference System

在某些遙感影像的應(yīng)用場(chǎng)景中，如果我們關(guān)注的區(qū)域正好處于兩景影像的交界處，如下圖中的象山港，那我們就需要將影像拼接起來(lái)才可以使用。

單張影像是這樣。

本文合并后是這樣。

1.準(zhǔn)備工作

相較于上一篇推送，我們這次為了實(shí)現(xiàn)遙感影像的鑲嵌拼接，我們使用到了兩個(gè)庫(kù)， rasterio和gdal。

import rasterio as rio

import gdal

先介紹一下我們實(shí)現(xiàn)兩組遙感影像拼接的思路，首先選取兩景相鄰的影像，分別得到他們的空間范圍，再得到兩景組合到一起之后的空間范圍，使用gdal新建一個(gè)tif文件(數(shù)據(jù)中轉(zhuǎn)用)，分別得到原來(lái)兩景影像在新建的tif文件中的起始位置，將對(duì)應(yīng)的數(shù)據(jù)寫(xiě)入新的tif文件中，即實(shí)現(xiàn)鑲嵌拼接。

上面說(shuō)的是兩景影像的拼接，如果是更多影像拼接同樣適用，但是現(xiàn)階段的方法如果拼接多的影像的話，需要的內(nèi)存空間很大，容易導(dǎo)致內(nèi)存溢出，感興趣的朋友可以思考一下如何高效實(shí)現(xiàn)多景影像的拼接。

其中還有兩處關(guān)鍵處理，一是如何去除重疊區(qū)域的無(wú)效信息，二是重疊區(qū)域的數(shù)據(jù)如何選擇。希望各位看官能從代碼里面找到答案。

2.動(dòng)起手來(lái)

得到輸入影像的四個(gè)角點(diǎn)。

def tiffileList2filename(tiffileList):

filename = []

prefix = []

for ifile in tiffileList:

file0 = ifile.split("\\")[-1]

prefix.append(os.path.join(ifile, file0))

filename.append(os.path.join(ifile, file0) + "_B1.TIF")

return filename, prefix

def get_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，如果是其他遙感數(shù)據(jù)，需對(duì)應(yīng)進(jìn)行修改。

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)

# 打開(kāi)一個(gè)tif文件

in_ds = gdal.Open(filename[0])

for i in range(len(bands)):

ibands = bands[i]

# 新建一個(gè)tif文件

driver = gdal.GetDriverByName('gtiff')

out_ds = driver.Create(ibands + 'mosaic.tif', cols, rows)

# 設(shè)置tif文件的投影

out_ds.SetProjection(in_ds.GetProjection())

out_band = out_ds.GetRasterBand(1)

# 設(shè)置新tif文件的地理變換

gt = list(in_ds.GetGeoTransform())

gt[0], gt[3] = left, top

out_ds.SetGeoTransform(gt)

# 對(duì)要拼接的影像進(jìn)行循環(huán)讀取

for ifile in prefix:

in_ds = gdal.Open(ifile + '_' + ibands + '.tif')

# 計(jì)算新建的tif文件及本次打開(kāi)的tif文件之間的坐標(biāo)漂移

trans = gdal.Transformer(in_ds, out_ds, [])

# 得到偏移起始點(diǎn)

success, xyz = trans.TransformPoint(False, 0, 0)

x, y, z = map(int, xyz)

# 讀取波段信息

fnBand = in_ds.GetRasterBand(1)

data = fnBand.ReadAsArray()

# 寫(xiě)入tif文件之前，最大值設(shè)置為255，這一步很關(guān)鍵

data = data / 65535 * 255

data[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)

del out_band, out_ds

file2read = ibands + 'mosaic.tif'

arr[i, :, :] = tiff.imread(file2read)

os.remove(file2read)

plot_rgb(arr, rgb=(0, 1, 2))

3.小結(jié) 大功告成！

總結(jié)

以上是生活随笔為你收集整理的python 拼接 遥感影像_如何用Python| 制作遥感影像拼接的全部?jī)?nèi)容，希望文章能夠幫你解決所遇到的問(wèn)題。

如果覺(jué)得生活随笔網(wǎng)站內(nèi)容還不錯(cuò)，歡迎將生活随笔推薦給好友。