[轉(zhuǎn)自：]（https://blog.csdn.net/tanlangqie/article/details/79486144）
GDAL是空間數(shù)據(jù)處理的開源包，支持多種數(shù)據(jù)格式的讀寫。遙感圖像是一種帶大地坐標的柵格數(shù)據(jù)，遙感圖像的柵格模型包含以下兩部分的內(nèi)容：

柵格矩陣：由正方形或者矩形柵格點組成，每個柵格點所對應的數(shù)值為該點的像元值，在遙感圖像中用于表示地物屬性值；遙感圖像有單波段與多波段，波段表示地物屬性的種類，每個波段表示地物一種屬性。

大地坐標：空間數(shù)據(jù)參考表示地圖的投影信息；仿射矩陣能將行列坐標映射到面坐標上。

GDAL讀寫遙感數(shù)據(jù)的代碼：

from osgeo import gdal
import os

class GRID:

    #讀圖像文件
    def read_img(self,filename):
        dataset=gdal.Open(filename)       #打開文件

        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ù)組，對應柵格矩陣

        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ù)類型必須有，因為要計算需要多大內(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遙感影像讀寫的基礎上，我們可以進行遙感圖像的各種公式計算和統(tǒng)計分析。
例如我們所熟知的計算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圖像

當然，這是理想的NDVI，實際處理NDVI還會遇到一些其他要處理的問題。例如NDVI值應該在區(qū)間[-1，1]內(nèi)，但實際中會出現(xiàn)大于1或小于-1的情況，或者某些像點是壞點，出現(xiàn)空值nan，需要進一步的配套處理。