您好,登錄后才能下訂單哦!
這篇文章主要介紹python+gdal+遙感圖像拼接的示例分析,文中介紹的非常詳細(xì),具有一定的參考價值,感興趣的小伙伴們一定要看完!
關(guān)于遙感圖像的鑲嵌,主要分為6大步驟:
step1:
1)對于每一幅圖像,計算其行與列;
2)獲取左上角X,Y
3)獲取像素寬和像素高
4)計算max X 和 min Y,切記像素高是負(fù)值
maxX1 = minX1 + (cols1 * pixelWidth)
minY1 = maxY1 + (rows1 * pixelHeight)
step2 :計算輸出圖像的min X ,max X,min Y,max Y
minX = min(minX1, minX2, …)
maxX = max(maxX1, maxX2, …)
y坐標(biāo)同理
step3:計算輸出圖像的行與列
cols = int((maxX – minX) / pixelWidth)
rows = int((maxY – minY) / abs(pixelHeight)
step 4:創(chuàng)建一個輸出圖像
driver.create()
step 5:
1)計算每幅圖像左上角坐標(biāo)在新圖像的偏移值
2)依次讀入每幅圖像的數(shù)據(jù)并利用1)計算的偏移值將其寫入新圖像中
step6 :對于輸出圖像
1)刷新磁盤并計算統(tǒng)計值
2)設(shè)置輸出圖像的幾何和投影信息
3)建立金字塔
下面附上筆者的代碼:
#mosica 兩張圖像 import os, sys, gdal from gdalconst import * os.chdir('c:/temp/****')#改變文件夾路徑 # 注冊gdal(required) gdal.AllRegister() # 讀入第一幅圖像 ds1 = gdal.Open('**.img') band1 = ds1.GetRasterBand(1) rows1 = ds1.RasterYSize cols1 = ds1.RasterXSize # 獲取圖像角點(diǎn)坐標(biāo) transform1 = ds1.GetGeoTransform() minX1 = transform1[0] maxY1 = transform1[3] pixelWidth2 = transform1[1] pixelHeight1 = transform1[5]#是負(fù)值(important) maxX1 = minX1 + (cols1 * pixelWidth2) minY1 = maxY1 + (rows1 * pixelHeight1) # 讀入第二幅圖像 ds2 = gdal.Open('**.img') band2 = ds2.GetRasterBand(1) rows2 = ds2.RasterYSize cols2 = ds2.RasterXSize # 獲取圖像角點(diǎn)坐標(biāo) transform2 = ds2.GetGeoTransform() minX2 = transform2[0] maxY2 = transform2[3] pixelWidth3 = transform2[1] pixelHeight2 = transform2[5] maxX2 = minX2 + (cols2 * pixelWidth3) minY2 = maxY2 + (rows2 * pixelHeight2) # 獲取輸出圖像坐標(biāo) minX = min(minX1, minX2) maxX = max(maxX1, maxX2) minY = min(minY1, minY2) maxY = max(maxY1, maxY2) #獲取輸出圖像的行與列 cols = int((maxX - minX) / pixelWidth2) rows = int((maxY - minY) / abs(pixelHeight1)) # 計算圖1左上角的偏移值(在輸出圖像中) xOffset1 = int((minX1 - minX) / pixelWidth2) yOffset1 = int((maxY1 - maxY) / pixelHeight1) # 計算圖2左上角的偏移值(在輸出圖像中) xOffset2 = int((minX2 - minX) / pixelWidth2) yOffset2 = int((maxY2 - maxY) / pixelHeight1) # 創(chuàng)建一個輸出圖像 driver = ds1.GetDriver() dsOut = driver.Create('mosiac.img', cols, rows, 1, band1.DataType)#1是bands,默認(rèn) bandOut = dsOut.GetRasterBand(1) # 讀圖1的數(shù)據(jù)并將其寫到輸出圖像中 data1 = band1.ReadAsArray(0, 0, cols1, rows1) bandOut.WriteArray(data1, xOffset1, yOffset1) #讀圖2的數(shù)據(jù)并將其寫到輸出圖像中 data2 = band2.ReadAsArray(0, 0, cols2, rows2) bandOut.WriteArray(data2, xOffset2, yOffset2) ''' 寫圖像步驟''' # 統(tǒng)計數(shù)據(jù) bandOut.FlushCache()#刷新磁盤 stats = bandOut.GetStatistics(0, 1)#第一個參數(shù)是1的話,是基于金字塔統(tǒng)計,第二個 #第二個參數(shù)是1的話:整幅圖像重度,不需要統(tǒng)計 # 設(shè)置輸出圖像的幾何信息和投影信息 geotransform = [minX, pixelWidth2, 0, maxY, 0, pixelHeight1] dsOut.SetGeoTransform(geotransform) dsOut.SetProjection(ds1.GetProjection()) # 建立輸出圖像的金字塔 gdal.SetConfigOption('HFA_USE_RRD', 'YES') dsOut.BuildOverviews(overviewlist=[2,4,8,16])#4層
補(bǔ)充知識:運(yùn)用Python的第三方庫:GDAL進(jìn)行遙感數(shù)據(jù)的讀寫
0 背景及配置環(huán)境
0.1 背景
GDAL(Geospatial Data Abstraction Library)是一個在X/MIT許可協(xié)議下的開源柵格空間數(shù)據(jù)轉(zhuǎn)換庫。它利用抽象數(shù)據(jù)模型來表達(dá)所支持的各種文件格式。它還有一系列命令行工具來進(jìn)行數(shù)據(jù)轉(zhuǎn)換和處理。
這個開源柵格空間數(shù)據(jù)轉(zhuǎn)換庫擁有許多和其他語言的接口,對于python,他有對應(yīng)的第三方包GDAL,下載安裝已在上篇文章中提到。
目的: 可以使用Python的第三方包:GDAL進(jìn)行遙感數(shù)據(jù)的讀寫,方便批處理。
0.2 配置環(huán)境
電腦系統(tǒng): win7x64
Python版本: 3.6.4
GDAL版本: 2.3.2
1 讀
1.1 TIFF格式
標(biāo)簽圖像文件格式(Tag Image File Format,簡寫為TIFF)是一種靈活的位圖格式,主要用來存儲包括照片和藝術(shù)圖在內(nèi)的圖像。它最初由Aldus公司與微軟公司一起為PostScript打印開發(fā)。TIFF與JPEG和PNG一起成為流行的高位彩色圖像格式。
TIFF文件以.tif為擴(kuò)展名。
def tif_read(tifpath, bandnum): """ Use GDAL to read data and transform them into arrays. :param tifpath:tif文件的路徑 :param bandnum:需要讀取的波段 :return:該波段的數(shù)據(jù),narray格式。len(narray)是行數(shù),len(narray[0])列數(shù) """ image = gdal.Open(tifpath) # 打開該圖像 if image == None: print(tifpath + "該tif不能打開!") return lie = image.RasterXSize # 柵格矩陣的列數(shù) hang = image.RasterYSize # 柵格矩陣的行數(shù) im_bands = image.RasterCount # 波段數(shù) im_proj = image.GetProjection() # 獲取投影信息 im_geotrans = image.GetGeoTransform() # 仿射矩陣 print('該tif:{}個行,{}個列,{}層波段, 取出第{}層.'.format(hang, lie, im_bands, bandnum)) band = image.GetRasterBand(bandnum) # Get the information of band num. band_array = band.ReadAsArray(0,0,lie,hang) # Getting data from zeroth rows and 0 columns # band_df = pd.DataFrame(band_array) del image # 減少冗余 return band_array, im_proj, im_geotrans
2 寫
2.1 TIFF格式
TIFF格式的數(shù)據(jù)格式有:Byete、int16、uint16、int32、uint32、float32、float64等7余種。
首先,要判斷數(shù)據(jù)的格式,才能按需求寫出。
def tif_write(self, filename, im_data, im_proj, im_geotrans): """ gdal數(shù)據(jù)類型包括 gdal.GDT_Byte, gdal.GDT_UInt16, gdal.GDT_Int16, gdal.GDT_UInt32, gdal.GDT_Int32, gdal.GDT_Float32, gdal.GDT_Float64 :param filename: 存出文件名 :param im_data: 輸入數(shù)據(jù) :param im_proj: 投影信息 :param im_geotrans: 放射變換信息 :return: 0 """ if 'int8' in im_data.dtype.name: # 判斷柵格數(shù)據(jù)的數(shù)據(jù)類型 datatype = gdal.GDT_Byte elif 'int16' in im_data.dtype.name: datatype = gdal.GDT_UInt16 else: datatype = gdal.GDT_Float32 # 判讀數(shù)組維數(shù) if len(im_data.shape) == 3: im_bands, im_height, im_width = im_data.shape else: im_bands, (im_height, im_width) = 1,im_data.shape # 多維或1.2維 #創(chuàng)建文件 driver = gdal.GetDriverByName("GTiff") #數(shù)據(jù)類型必須有,因為要計算需要多大內(nèi)存空間 dataset = driver.Create(filename, im_width, im_height, im_bands, datatype) dataset.SetGeoTransform(im_geotrans) #寫入仿射變換參數(shù) dataset.SetProjection(im_proj) #寫入投影 if im_bands == 1: dataset.GetRasterBand(1).WriteArray(im_data) #寫入數(shù)組數(shù)據(jù) else: for i in range(im_bands): dataset.GetRasterBand(i+1).WriteArray(im_data[i]) del dataset
3 展示
3.1 TIFF格式
# 這個展示的效果并不是太好,當(dāng)做示意圖用 def tif_display(self,im_data): """ :param im_data: 影像數(shù)據(jù),narray :return: 展出影像 """ # plt.imshow(im_data,'gray') # 必須規(guī)定為顯示的為什么圖像 plt.imshow(im_data) # 必須規(guī)定為顯示的為什么圖像 plt.xticks([]), plt.yticks([]) # 隱藏坐標(biāo)線 plt.show() # 顯示出來,不要也可以,但是一般都要了
以上是“python+gdal+遙感圖像拼接的示例分析”這篇文章的所有內(nèi)容,感謝各位的閱讀!希望分享的內(nèi)容對大家有幫助,更多相關(guān)知識,歡迎關(guān)注億速云行業(yè)資訊頻道!
免責(zé)聲明:本站發(fā)布的內(nèi)容(圖片、視頻和文字)以原創(chuàng)、轉(zhuǎn)載和分享為主,文章觀點(diǎn)不代表本網(wǎng)站立場,如果涉及侵權(quán)請聯(lián)系站長郵箱:is@yisu.com進(jìn)行舉報,并提供相關(guān)證據(jù),一經(jīng)查實(shí),將立刻刪除涉嫌侵權(quán)內(nèi)容。