您好,登錄后才能下訂單哦!
這篇文章主要介紹了Python空間數(shù)據(jù)處理之GDAL如何讀寫遙感圖像,具有一定借鑒價(jià)值,感興趣的朋友可以參考下,希望大家閱讀完這篇文章之后大有收獲,下面讓小編帶著大家一起了解一下。
GDAL是空間數(shù)據(jù)處理的開(kāi)源包,支持多種數(shù)據(jù)格式的讀寫。遙感圖像是一種帶大地坐標(biāo)的柵格數(shù)據(jù),遙感圖像的柵格模型包含以下兩部分的內(nèi)容:
柵格矩陣:由正方形或者矩形柵格點(diǎn)組成,每個(gè)柵格點(diǎn)所對(duì)應(yīng)的數(shù)值為該點(diǎn)的像元值,在遙感圖像中用于表示地物屬性值;遙感圖像有單波段與多波段,波段表示地物屬性的種類,每個(gè)波段表示地物一種屬性。
大地坐標(biāo):空間數(shù)據(jù)參考表示地圖的投影信息;仿射矩陣能將行列坐標(biāo)映射到面坐標(biāo)上。
GDAL讀寫遙感數(shù)據(jù)的代碼:
from osgeo import gdal import os class GRID: #讀圖像文件 def read_img(self,filename): dataset=gdal.Open(filename) #打開(kāi)文件 im_width = dataset.RasterXSize #柵格矩陣的列數(shù) im_height = dataset.RasterYSize #柵格矩陣的行數(shù) im_geotrans = dataset.GetGeoTransform() #仿射矩陣 im_proj = dataset.GetProjection() #地圖投影信息 im_data = dataset.ReadAsArray(0,0,im_width,im_height) #將數(shù)據(jù)寫成數(shù)組,對(duì)應(yīng)柵格矩陣 del dataset return im_proj,im_geotrans,im_data #寫文件,以寫成tif為例 def write_img(self,filename,im_proj,im_geotrans,im_data): #gdal數(shù)據(jù)類型包括 #gdal.GDT_Byte, #gdal .GDT_UInt16, gdal.GDT_Int16, gdal.GDT_UInt32, gdal.GDT_Int32, #gdal.GDT_Float32, gdal.GDT_Float64 #判斷柵格數(shù)據(jù)的數(shù)據(jù)類型 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 #判讀數(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 #創(chuàng)建文件 driver = gdal.GetDriverByName("GTiff") #數(shù)據(jù)類型必須有,因?yàn)橐?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 if __name__ == "__main__": os.chdir(r'D:\Python_Practice') #切換路徑到待處理圖像所在文件夾 run = GRID() proj,geotrans,data = run.read_img('LC81230402013164LGN00.tif') #讀數(shù)據(jù) print proj print geotrans print data print data.shape run.write_img('LC81230402013164LGN00_Rewrite.tif',proj,geotrans,data) #寫數(shù)據(jù)
在GDAL遙感影像讀寫的基礎(chǔ)上,我們可以進(jìn)行遙感圖像的各種公式計(jì)算和統(tǒng)計(jì)分析。
例如我們所熟知的計(jì)算NDVI(歸一化植被指數(shù)),只要在以上代碼倒數(shù)第二行中插入代碼:
import numpy as np data = data.astype(np.float) ndvi = (data[3]-data[2])/(data[3]+data[2]) #3為近紅外波段;2為紅波段 run.write_img('LC81230402013164LGN00_ndvi.tif',proj,geotrans,ndvi) #寫為ndvi圖像
當(dāng)然,這是理想的NDVI,實(shí)際處理NDVI還會(huì)遇到一些其他要處理的問(wèn)題。例如NDVI值應(yīng)該在區(qū)間[-1,1]內(nèi),但實(shí)際中會(huì)出現(xiàn)大于1或小于-1的情況,或者某些像點(diǎn)是壞點(diǎn),出現(xiàn)空值nan,需要進(jìn)一步的配套處理。
感謝你能夠認(rèn)真閱讀完這篇文章,希望小編分享的“Python空間數(shù)據(jù)處理之GDAL如何讀寫遙感圖像”這篇文章對(duì)大家有幫助,同時(shí)也希望大家多多支持億速云,關(guān)注億速云行業(yè)資訊頻道,更多相關(guān)知識(shí)等著你來(lái)學(xué)習(xí)!
免責(zé)聲明:本站發(fā)布的內(nèi)容(圖片、視頻和文字)以原創(chuàng)、轉(zhuǎn)載和分享為主,文章觀點(diǎn)不代表本網(wǎng)站立場(chǎng),如果涉及侵權(quán)請(qǐng)聯(lián)系站長(zhǎng)郵箱:is@yisu.com進(jìn)行舉報(bào),并提供相關(guān)證據(jù),一經(jīng)查實(shí),將立刻刪除涉嫌侵權(quán)內(nèi)容。