作為攝影測量與遙感的從業者，筆者最近開始深入研究gdal，為工作打基礎！個人覺得gdal也是沒有什么技術含量，調用別人的api。但是想想這也是算法應用的一個技能，多學無害！

關于遙感圖像的鑲嵌，主要分為6大步驟：

step1：

1）對于每一幅圖像，計算其行與列；

2）獲取左上角X,Y

3）獲取像素寬和像素高

4）計算max X 和 min Y，切記像素高是負值

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坐標同理

step3：計算輸出圖像的行與列

cols = int((maxX ? minX) / pixelWidth)rows = int((maxY ? minY) / abs(pixelHeight)

step 4：創建一個輸出圖像

driver.create()

step 5：

1）計算每幅圖像左上角坐標在新圖像的偏移值

2）依次讀入每幅圖像的數據并利用1)計算的偏移值將其寫入新圖像中

step6 ：對于輸出圖像

1）刷新磁盤并計算統計值

2）設置輸出圖像的幾何和投影信息

3）建立金字塔

下面附上筆者的代碼：

#mosica 兩張圖像import os, sys, gdalfrom gdalconst import *os.chdir(’c:/temp/****’)#改變文件夾路徑# 注冊gdal(required)gdal.AllRegister()# 讀入第一幅圖像ds1 = gdal.Open(’**.img’)band1 = ds1.GetRasterBand(1)rows1 = ds1.RasterYSizecols1 = ds1.RasterXSize# 獲取圖像角點坐標transform1 = ds1.GetGeoTransform()minX1 = transform1[0]maxY1 = transform1[3]pixelWidth1 = transform1[1]pixelHeight1 = transform1[5]#是負值（important）maxX1 = minX1 + (cols1 * pixelWidth1)minY1 = maxY1 + (rows1 * pixelHeight1)# 讀入第二幅圖像ds2 = gdal.Open(’**.img’)band2 = ds2.GetRasterBand(1)rows2 = ds2.RasterYSizecols2 = ds2.RasterXSize# 獲取圖像角點坐標transform2 = ds2.GetGeoTransform()minX2 = transform2[0]maxY2 = transform2[3]pixelWidth2 = transform2[1]pixelHeight2 = transform2[5]maxX2 = minX2 + (cols2 * pixelWidth2)minY2 = maxY2 + (rows2 * pixelHeight2)# 獲取輸出圖像坐標minX = min(minX1, minX2)maxX = max(maxX1, maxX2)minY = min(minY1, minY2)maxY = max(maxY1, maxY2)#獲取輸出圖像的行與列cols = int((maxX - minX) / pixelWidth1)rows = int((maxY - minY) / abs(pixelHeight1))# 計算圖1左上角的偏移值（在輸出圖像中）xOffset1 = int((minX1 - minX) / pixelWidth1)yOffset1 = int((maxY1 - maxY) / pixelHeight1)# 計算圖2左上角的偏移值（在輸出圖像中）xOffset2 = int((minX2 - minX) / pixelWidth1)yOffset2 = int((maxY2 - maxY) / pixelHeight1)# 創建一個輸出圖像driver = ds1.GetDriver()dsOut = driver.Create(’mosiac.img’, cols, rows, 1, band1.DataType)#1是bands，默認bandOut = dsOut.GetRasterBand(1)# 讀圖1的數據并將其寫到輸出圖像中data1 = band1.ReadAsArray(0, 0, cols1, rows1)bandOut.WriteArray(data1, xOffset1, yOffset1)#讀圖2的數據并將其寫到輸出圖像中data2 = band2.ReadAsArray(0, 0, cols2, rows2)bandOut.WriteArray(data2, xOffset2, yOffset2)’’’ 寫圖像步驟’’’# 統計數據bandOut.FlushCache()#刷新磁盤stats = bandOut.GetStatistics(0, 1)#第一個參數是1的話，是基于金字塔統計，第二個#第二個參數是1的話：整幅圖像重度，不需要統計# 設置輸出圖像的幾何信息和投影信息geotransform = [minX, pixelWidth1, 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層

補充知識：運用Python的第三方庫：GDAL進行遙感數據的讀寫

0 背景及配置環境

0.1 背景

GDAL(Geospatial Data Abstraction Library)是一個在X/MIT許可協議下的開源柵格空間數據轉換庫。它利用抽象數據模型來表達所支持的各種文件格式。它還有一系列命令行工具來進行數據轉換和處理。

這個開源柵格空間數據轉換庫擁有許多和其他語言的接口，對于python,他有對應的第三方包GDAL，下載安裝已在上篇文章中提到。

目的: 可以使用Python的第三方包：GDAL進行遙感數據的讀寫，方便批處理。

0.2 配置環境

電腦系統： win7x64Python版本: 3.6.4GDAL版本: 2.3.2

1 讀

1.1 TIFF格式

標簽圖像文件格式（Tag Image File Format，簡寫為TIFF）是一種靈活的位圖格式，主要用來存儲包括照片和藝術圖在內的圖像。它最初由Aldus公司與微軟公司一起為PostScript打印開發。TIFF與JPEG和PNG一起成為流行的高位彩色圖像格式。

TIFF文件以.tif為擴展名。

def tif_read(tifpath, bandnum): ''' Use GDAL to read data and transform them into arrays. :param tifpath:tif文件的路徑 :param bandnum:需要讀取的波段 :return:該波段的數據，narray格式。len(narray)是行數，len(narray[0])列數 ''' image = gdal.Open(tifpath) # 打開該圖像 if image == None: print(tifpath + '該tif不能打開!') return lie = image.RasterXSize # 柵格矩陣的列數 hang = image.RasterYSize # 柵格矩陣的行數 im_bands = image.RasterCount # 波段數 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格式的數據格式有：Byete、int16、uint16、int32、uint32、float32、float64等7余種。

首先，要判斷數據的格式，才能按需求寫出。

def tif_write(self, filename, im_data, im_proj, im_geotrans): ''' gdal數據類型包括 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: 輸入數據 :param im_proj: 投影信息 :param im_geotrans: 放射變換信息 :return: 0 ''' if ’int8’ in im_data.dtype.name: # 判斷柵格數據的數據類型 datatype = gdal.GDT_Byte elif ’int16’ in im_data.dtype.name: datatype = gdal.GDT_UInt16 else: datatype = gdal.GDT_Float32 # 判讀數組維數 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維 #創建文件 driver = gdal.GetDriverByName('GTiff') #數據類型必須有，因為要計算需要多大內存空間 dataset = driver.Create(filename, im_width, im_height, im_bands, datatype) dataset.SetGeoTransform(im_geotrans) #寫入仿射變換參數 dataset.SetProjection(im_proj) #寫入投影 if im_bands == 1: dataset.GetRasterBand(1).WriteArray(im_data) #寫入數組數據 else: for i in range(im_bands): dataset.GetRasterBand(i+1).WriteArray(im_data[i]) del dataset

3 展示

3.1 TIFF格式

# 這個展示的效果并不是太好，當做示意圖用 def tif_display(self,im_data): ''' :param im_data: 影像數據,narray :return: 展出影像 ''' # plt.imshow(im_data,’gray’) # 必須規定為顯示的為什么圖像 plt.imshow(im_data) # 必須規定為顯示的為什么圖像 plt.xticks([]), plt.yticks([]) # 隱藏坐標線 plt.show() # 顯示出來，不要也可以，但是一般都要了

以上這篇python+gdal+遙感圖像拼接(mosaic)的實例就是小編分享給大家的全部內容了，希望能給大家一個參考，也希望大家多多支持好吧啦網。