python gdal教程之：用gdal讀取柵格數據

GDAL原生支援超過100種柵格資料類型，涵蓋所有主流GIS與RS資料格式，包括

 ArcInfo grids, ArcSDE raster, Imagine, Idrisi, ENVI, GRASS, GeoTIFF 

 HDF4, HDF5oo DEM 

 ECW, MrSID 

 TIFF, JPEG, JPEG2000, PNG, GIF, BMP 

完整的支援清單可以參考http://www.gdal.org/formats_list.htmlcom （1.5以前）：import gdal, gdalconst

新版本（1.6以後）：from osgeo import gdal, gdalconst

gdal和gdalconst最好都要導入，其中modgdalconst中的常數都加了前綴，力圖與其他的uleule衝突最小。所以對gdalconst你可以直接這樣導入：from osgeo.gdalconst import *

GDAL資料驅動，與OGR資料驅動類似，需要先建立某一類型的資料驅動，再建立回應的柵格資料集。

一次註冊所有的資料驅動，但是只能讀不能寫：gdal.AllRegister()

單獨註冊某一類型的資料驅動，這樣的話可以讀也可以寫，可以新建資料集：

driver = gdal.GetDriverByName('HFA') 

driver.Register()

開啟現有的柵格資料集：

   fn = 'aster.img' 

 ds is None: 

      print 'Could not open ' + fn 

       sys.exit(1)

RasterXSize 

   rows = ds.RasterYSize 

   bands = ds.RasterCount

注意後面沒有括號，因為他們是屬性(properties)不是方法(methods))一個li​​st，儲存著柵格資料集的地理座標資訊

   adfGeoTransform[0] /* top left x 左上角x座標*/ 

   adfGeoTransform[1] /* w--e pixel resolution 東西方向上的像素率*/ 

   adfGeoTransform[2] /* rotation, 0 if image is "north up" 如果北邊朝上，地圖的旋轉角*/ 

   adfGeoTransform[3] / top

   adfGeoTransform[3] / top

   adfGeoTransform[3] / top

   adfGeoTransform[4] /* rotation, 0 if image is "north up" 如果北邊朝上，地圖的旋轉角度*/ 

   adfGeoTransform[5] /* n-s pixel resolution 分辨率/注意柵格資料集的座標一般都是以左上角為基準的。

下面的例子是從一個柵格資料集中取出Geotransform作為一個list，然後讀取其中的資料

   geotransform = ds.GetGeoTransform() 

   originX = geotransgin[0] geotransgin = geotransform[3] 

   pixelWidth = geotransform[1] 

   pixelHeight = geotransform[5]

計算某一座標像素的相對位置(pixel 相對位置(pixel 的相對位置)，也就是該座標對應的像素對應以像素數計算，計算公式如下：

xOffset = int((x – originX) / pixelWidth) 

yOffset = int((y – originY) / pixelHeight)

兩步驟

先讀取一個波段(band)：GetRasterBand()，其參數為波段的索引號碼

然後用ReadAsArray(, , , )，讀出從(xoff,yoff)開始，大小為(xsize,ysize)的矩陣。如果將矩陣大小設為1X1，就是讀取一個像素了。但這方法只能將讀出的資料放到矩陣中，就算只讀取一個像素也是一樣。例如：

band = ds.GetRasterBand(1)

data = band.ReadAsArray(xOffset, yOffset, 1, 1)

如果想一次讀取一整張圖，那麼將offset都設定為0，size則設定為整個圖幅的size，例如：

data = band.ReadAsArray(0, 0, cols, rows)

但是要注意，從data讀取某一像素的值，必須要用data[ yoff, xoff]。注意不要搞反了。數學中的矩陣是[row,col]，而這裡恰恰相反！這裡面row對應y軸，col對應x軸。

注意在適當的時候釋放內存，例如band = None 或dataset = None。尤其圖很大的時候

如何更有效率的讀取柵格資料？顯然一個一個的讀取效率非常低，將整個柵格資料集都塞進二維數組也不是個好辦法，因為這樣佔的記憶體還是很多。更好的方法是按塊(block)來存取數據，只把要使用的那一塊放進記憶體。本週的範例程式碼中有一個utils模組，可以讀取block大小。

例如：

   import utils 

   blockSize = utils.GetBlockSize(band)

   xBlockSize = blockSize[0]

平鋪(tiled)，即柵格資料以block儲存。有的格式，例如GeoTiff沒有平鋪，一行是block。 Erdas imagine格式則以64x64像素平鋪。

如果一行是一個block，那麼按行讀取是比較節省資源的。

如果是平舖的資料結構，那麼設定ReadAsArray()的參數值，讓它一次只讀入一個block，就是效率最高的方法了。例如：

rows = 13, cols = 11, xBSize = 5, yBSize = 5 

for i in range(0, rows, yBSize): 

if i + yBSize

else : 

        numRows = rows – i 

    for j in range(0, cols, xBSize): 

       numCols = xBSize 

        else: 

            numCols = colsnumCols = cols – j)

        data = band.ReadAsArray(j, i, numCols, numRows)

這段程式碼具有一般性，可以時常拿來用的。

下面介紹一點二維陣列的處理技巧

這裡要用到兩個函式庫，Numeric和numpy。 Numeric比較老了，FWTools用它。自己安裝設定的話還是配功能更強的numpy。

資料型別轉換：

data = band.ReadAsArray(j, i, nCols, nRows) 

data = data.astype(Numeric.Float) # Numeric 

data = data.as Numeric 

data = .as

或簡單點只寫一句

data = band.ReadAsArray(j, i, nCols, nRows).astype(Numeric.Float)

膜mask

這是Numeric和numpy掩膜的功能，輸入一個數組和條件，輸出一個二值數組。例如

mask = Numeric.greater(data, 0)mask = Numeric.greater(data, 0) 

>>> a = Numeric.array([0, 4, 6, 0, 2]) 

> > print a 

[0 4 6 0 2] 

>>> mask = Numeric.greater(a, 0) 

>>> print mask 

[0 1 1 0 1]

和求數字>>> a = Numeric.array([0, 4, 6, 0, 2]) 

>>> print a>>> print a 

>

[0 4 6 0 2] 

>

[0 4 6 0 2] 

>>>> print Numeric. sum(a) 

12

如果是二維數組，那麼sum就會回傳一個一維數組

>>> b = Numeric.array([a, [5, 10, 0, 3, 0]] ) 

>>> print b 

[[ 0      4  6  0  2] 

[ 5 10  0

[ 5 14  6  3  2]

所以，二維陣列的求和要這樣

>>>> print Numeric.sum(Numeric.sum(b)) 

30

統計大於0的像素個數，可以聯合運用mask和sum兩個函數

>>> print a 

[0 4 6 0 2] 

>>> mask = Numeric.greater(a, 0) 

>>> print mask

[0 1 1 0 1] 

>>> print Numeric.sum(mask) 

3

 以上是python gdal教程相關內容請關注PHP中文網（www.php.cn）！