|
1. GDAL庫(kù)介紹可能你不玩GIS，不懂這個(gè)庫(kù)到底有什么用，或者和python有什么關(guān)系。但是你要玩GIS，RS，你就應(yīng)當(dāng)知道這個(gè)庫(kù)的價(jià)值。就算你不玩GIS，我想這個(gè)庫(kù)對(duì)你也應(yīng)該有致命的吸引力。為什么？看下面的介紹吧！ 先看看這段GDAL主頁(yè)上的英文介紹吧！ 有人說(shuō)我又不玩GIS。不錯(cuò)，但是，你即使不玩GIS，這個(gè)庫(kù)也是滿有用的。首先，哪個(gè)庫(kù)支持這么多柵格（圖片）格式，哪個(gè)庫(kù)在C/C++/python/ruby/VB/java/C#(這個(gè)暫時(shí)不完全支持)下都能用，而且都一樣用？退一步講，3S軟件又不一定要用在3S下（很多醫(yī)學(xué)影像就是用PCI軟件來(lái)處理的）。再退一步，你的生活即使和3S一點(diǎn)關(guān)系都沒(méi)有，柵格數(shù)據(jù)又不單單只有GIS下才用到。你大可用這個(gè)庫(kù)來(lái)讀取jpg，gif，tif，xpm等格式。而且對(duì)各種格式支持得不是一般的好，很大一部分非標(biāo)準(zhǔn)格式照樣支持得非常好。我曾經(jīng)在java下玩過(guò)jai，以及一系列jai的擴(kuò)展庫(kù)，一些圖像格式在很多圖片瀏覽器中都可以正確讀取（有的甚至不是非標(biāo)準(zhǔn)格式），用jai死活就讀不出來(lái)！ 這個(gè)庫(kù)的python版和其他的python庫(kù)結(jié)合的很好。最直接、明顯的支持是使用Numeric庫(kù)來(lái)進(jìn)行數(shù)據(jù)讀取和操作。各種矩陣魔術(shù)可以發(fā)揮得淋漓盡致（圖像其實(shí)就是矩陣）。而且按我的觀點(diǎn)，python對(duì)矩陣的操作比其他的語(yǔ)言有明顯的優(yōu)勢(shì)。寫出來(lái)的東西比其他語(yǔ)言寫出來(lái)的短小的多，而且好看得多。并且python的弱類型在處理柵格數(shù)據(jù)格式類型的時(shí)候代碼量比強(qiáng)類型的語(yǔ)言少了數(shù)倍（不用double，byte，short等等分開(kāi)處理，這簡(jiǎn)直就是先天上的優(yōu)勢(shì)）。所以我就喜歡用python做圖像的處理。所以就連GIS界的微軟ESRI也直接在ARCGIS9中用python來(lái)作柵格數(shù)據(jù)的導(dǎo)入導(dǎo)出。一句話，真是太方便啦！
2. 安裝
2.1. windows下的安裝官方安裝文檔在這里。下面是我自己的實(shí)踐步驟： 先去http://www./dl/下一個(gè)版本，解壓。 打開(kāi)控制臺(tái)，輸入： “D:\Program Files\Microsoft Visual Studio .NET 2003\Vc7\bin\vcvars32.bat" 注冊(cè)vc的編譯環(huán)境。 打 開(kāi)gdal文件夾下的nmake.opt修改GDAL_HOME = "C:\warmerda\bld"把路徑改到需要把gdal安裝的地方。不改也可以。這里需要添加python支持，所以修改PY_INST_DIR = $(GDAL_HOME)\pymod把路徑改成python下的Lib\site-packages文件夾下。PYDIR = "C:\Software\Python24" 改成python的安裝路徑。 下面的參數(shù)愛(ài)改什么就把前面的#刪除（要看您有沒(méi)有那些庫(kù)的源碼），注意一下路徑就可以了。我是都沒(méi)改。 后面就依次運(yùn)行
Toggle line numbers
1 nmake /f makefile.vc
2 nmake /f makefile.vc install
3 nmake /f makefile.vc devinstall
最后最后，還要去GDAL_HOME目錄下的bin文件夾下把gdal13.dll （也有可能是gdal12.dll）copy到PY_INST_DIR路徑下 到此處就完成安裝gdal（python）的工作。 最后需要注意一下，gdal在.net2005下只能順利編譯1.2，1.3以上版本不能順利編譯，有一個(gè)地方指針轉(zhuǎn)換出錯(cuò)。可能是2005的編譯器比以往的嚴(yán)厲一點(diǎn)吧。另外,安裝了QGIS,對(duì)編譯也有一些影響,主要是proj庫(kù)的沖突,導(dǎo)致一個(gè)找不到"d:/program.obj"文件的錯(cuò)誤,如果你有靜態(tài)編譯過(guò)proj,那么你可以打開(kāi)nmake.opt修改有關(guān)proj的設(shè)置,如果搞不定,就卸載QGIS,然后編譯,編譯后再安裝QGIS.呵呵,還好QGIS的體積沒(méi)有ArcGIS那么可怕.
2.2. linux下的安裝linux下的安裝就更簡(jiǎn)單了。直接
Toggle line numbers
1 ./configure
2 make
3 su
4 make install
5 ldconfig
就ok(默認(rèn)就已經(jīng)支持python)。當(dāng)然在第一步的時(shí)候需要看看是否依賴的庫(kù)都安裝了。如果缺少，就去安裝一個(gè)。如果對(duì)configure的條件不理解，就用./configure --help看看具體情況。
2.3. 安裝其他驅(qū)動(dòng)這里講一個(gè)安裝hdf4的驅(qū)動(dòng)的例子（默認(rèn)情況下gdal是不安裝hdf4的），其他驅(qū)動(dòng)應(yīng)該和這個(gè)也差不了多少吧，可以作為其他的參考。完整步驟如下： 在windows下的安裝： 從ftp://ftp.ncsa.uiuc.edu/HDF/HDF/HDF%5FCurrent/bin/windows/下載42r1-win.ZIP，解壓。 編輯gdal根目錄下的nmake.opt，找到“# Uncomment the following and update to enable NCSA HDF Release 4 support.”這一行 把下面兩行前面的#去掉，然后改成：
HDF4_DIR = D:\warmerda\42r1-win\release #HDF4_LIB = /LIBPATH:$(HDF4_DIR)\lib hd421m.lib HDF4_LIB = $(HDF4_DIR)\dll\hd421m.lib $(HDF4_DIR)\dll\hm421m.lib \
用HDF4_LIB=/LIBPATH：這種形式似乎可以建立gdal的庫(kù)，但是往下編譯會(huì)出錯(cuò)。而且要把$(HDF4_DIR)\dll和$(HDF4_DIR)\lib拷貝到同一個(gè)目錄下，不然會(huì)提示找不到庫(kù) 你也可以試一試在D:\Program Files\Microsoft Visual Studio .NET 2003\Common7\Tools\vsvars32.bat文件中添加HDF4_LIB路徑到“@set LIB=”這行的末尾（不要忘記;的分割符）。 然后找一下"INC="這行,把 -I$(HDF4_DIR)\include加到下一行的末尾（應(yīng)該也可以在vsvars32.bat中添加路徑，不過(guò)要重啟命令行）。 然后編譯吧！祝你好運(yùn)。 注意：上面的HDF4_DIR 是我本機(jī)的路徑，你要根據(jù)你自己的路徑進(jìn)行設(shè)置（想起我的一個(gè)老師說(shuō)過(guò)的話：“抄人家的作業(yè)可以，不要連名字也一起抄走啊 ” 編譯成功后，要HDF4能運(yùn)行，還需要兩個(gè)庫(kù)，一個(gè)是zlib,一個(gè)是szip，可以到下面兩個(gè)鏈接去下載一個(gè) ftp://hdf.ncsa.uiuc.edu/lib-external/zlib/1.2/bin ftp://hdf.ncsa.uiuc.edu/lib-external/szip/2.0/bin/windows 把這兩個(gè)庫(kù)下載后解壓，然后設(shè)置PATH系統(tǒng)變量，使得它們?cè)谀J(rèn)狀態(tài)下也可以被動(dòng)態(tài)鏈接成功 。 在Linux下比在Windows下簡(jiǎn)單： 只要用./configure --help察看一下打開(kāi)HDF4的編譯開(kāi)關(guān)（包括庫(kù)路徑，頭文件路徑等，看清楚），然后在./configure 后面加上那個(gè)開(kāi)關(guān)以及hdf4的安裝路徑后就可以了。在configure后gdal會(huì)提示是否支持HDF4。 編譯后也要把zlib和szip的動(dòng)態(tài)鏈接庫(kù)設(shè)置好 。 到此你已經(jīng)可以用C/C++來(lái)操作gdal讀寫hdf4的格式了！ 最后，為了讓Python能夠支持hdf的讀寫，別忘了把重新生成安裝后的pymod目錄下的內(nèi)容，還有g(shù)dal13，還有那兩個(gè)hdf的庫(kù)，還有zlib，szip的庫(kù)拷貝到Python的Lib\site-packages目錄下 。
2.4. 下載如果你實(shí)在玩不轉(zhuǎn)，可以在這里下載已經(jīng)編譯好的gdal1.3.2程序庫(kù) 以及其依賴的其他庫(kù)，其中包括hdf4，hdf5支持，以及proj，geos插件。注意，這里的geos是靜態(tài)鏈接的，注意版權(quán)（geos是LGPL的license）。hdf4和hdf5用的是release版本。這里是我的nmake配置文件，你可以對(duì)照你的實(shí)際情況參考一下。
3. 快速開(kāi)始其實(shí)在主站的教程里已經(jīng)有python的示例了。但是我們還是按照自己的思路來(lái)開(kāi)始吧。 第一步就是打開(kāi)一個(gè)數(shù)據(jù)集。對(duì)于“數(shù)據(jù)集”這個(gè)名詞大家可能不會(huì)太習(xí)慣，但是對(duì)于一般的格式來(lái)說(shuō)，一個(gè)“數(shù)據(jù)集”就是一個(gè)文件，比如一個(gè)gif文件就是一個(gè)以gif為擴(kuò)展名的文件。但是對(duì)于眾多RS數(shù)據(jù)來(lái)說(shuō)，一個(gè)數(shù)據(jù)集包含的絕對(duì)不僅僅是一個(gè)文件。對(duì)于很多RS數(shù)據(jù)，他們把一張圖像分成數(shù)個(gè)圖像文件，然后放在一個(gè)文件夾中，用一些額外的文件來(lái)組織它們之間的關(guān)系，形成一個(gè)“數(shù)據(jù)集”。如果你不理解，那么就算了，當(dāng)成jpg或者gif文件好了。 下面我們打開(kāi)一個(gè)tiff文件(GeoTIFF)。這個(gè)文件是我從GRASS的示例數(shù)據(jù)spearfish中導(dǎo)出的一個(gè)同名影像數(shù)據(jù)。
Toggle line numbers
1 >>> import gdal
2 >>> dataset = gdal.Open("j:/gisdata/gtif/spot.tif")
3 >>> dir(dataset)
4 [‘AddBand‘, ‘AdviseRead‘, ‘BuildOverviews‘, ‘FlushCache‘, ‘GetDescription‘, ‘Get
5 Driver‘, ‘GetGCPCount‘, ‘GetGCPProjection‘, ‘GetGCPs‘, ‘GetGeoTransform‘, ‘GetMe
6 tadata‘, ‘GetProjection‘, ‘GetProjectionRef‘, ‘GetRasterBand‘, ‘GetSubDatasets‘,
7 ‘RasterCount‘, ‘RasterXSize‘, ‘RasterYSize‘, ‘ReadAsArray‘, ‘ReadRaster‘, ‘Refr
8 eshBandInfo‘, ‘SetDescription‘, ‘SetGCPs‘, ‘SetGeoTransform‘, ‘SetMetadata‘, ‘Se
9 tProjection‘, ‘WriteRaster‘, ‘__del__‘, ‘__doc__‘, ‘__init__‘, ‘__module__‘, ‘_b
10 and‘, ‘_o‘]
11 >>>
這樣我們就打開(kāi)了這個(gè)文件。并且我們可以看到可以供我們調(diào)用的函數(shù)們(更具體的API列表可以看這里)?，F(xiàn)在我們不做修改，不做添加，所以只要帶有Set開(kāi)頭的函數(shù)以及有Write開(kāi)頭的函數(shù)我們暫時(shí)都不管。因?yàn)镽S影像必然要和地理上的位置掛上鉤，才能把圖像正確鋪展到一個(gè)坐標(biāo)系中。其中的信息和對(duì)應(yīng)關(guān)系有點(diǎn)復(fù)雜，不適合在快速開(kāi)始中介紹，我們暫時(shí)也先不管。這里需要注意的就是幾個(gè)函數(shù)。 GetDescription 獲得柵格的描述信息。
Toggle line numbers
1 >>> dataset.GetDescription()
2 ‘j:/gisdata/gtif/spot.tif‘
3 >>>
看來(lái)這里的圖像描述是圖像的路徑名，但是這是和各種不同數(shù)據(jù)集相關(guān)的，不同數(shù)據(jù)集可能有不同的描述。這要看讀取驅(qū)動(dòng)的實(shí)現(xiàn)作者的高興了。 RasterCount 獲得柵格數(shù)據(jù)集的波段數(shù)。 GetRasterBand 獲得柵格數(shù)據(jù)集的波段。
Toggle line numbers
1 >>> dataset.RasterCount
2 1
3 >>> band = dataset.GetRasterBand(1)
4 >>>
這里需要解釋的是Band這個(gè)詞。這個(gè)詞可以翻譯成“波段”，“通道”等等。我這里把它統(tǒng)一稱為“波段”。因?yàn)檫b感衛(wèi)星的傳感器有很多個(gè)。一個(gè)傳感器只負(fù)責(zé)接收一個(gè)頻率范圍的地物反射光波，一個(gè)頻率范圍的光波記錄稱為一個(gè)波段。是不是暈了？其實(shí)說(shuō)得簡(jiǎn)單一點(diǎn)。其實(shí)你可以把波段看成紅綠藍(lán)幾種顏色。圖像不是分RGB三色嗎？把R，G，B值都提取出來(lái)成為三個(gè)表。R值表就是波段一，G值表就是波段二，B值表就是波段三。 這里我們看到這張圖只有一個(gè)波段(一種顏色)。就可以把它看成是一個(gè)灰度圖(類似黑白照片)。如果RasterCount是3，就有可能是彩色圖。如果RasterCount是比3大的數(shù)，恭喜你，你看到一張遙感影像。有很多衛(wèi)星的傳感器大于3個(gè)，比如TM就有7個(gè)波段，不僅有可見(jiàn)光，還有紅外等其他非可見(jiàn)光。，所以，波段一般比RGB能表達(dá)的豐富地多。不過(guò)這樣一來(lái)就需要我們從中挑出3個(gè)波段然后組合成RGB，當(dāng)然這樣就有可能使圖像顯示出來(lái)的東西不像平常我們看到的那樣。這樣安排是因?yàn)閷?duì)科學(xué)有幫助(一些波段在科學(xué)家眼里比真實(shí)的彩色照片有價(jià)值)。不理解就跳過(guò)，很正常，我第一次聽(tīng)這種東西也覺(jué)得很玄:) 這里我們獲取了第一個(gè)波段(紅色值組成的表)。注意！這里的波段獲取和通常的C數(shù)組獲取不一樣，開(kāi)始是1不是0。獲取了波段，我們就可以在下面的操作中讀取這個(gè)波段的所有數(shù)值。 RasterXSize 圖像的寬度(X方向上的像素個(gè)數(shù)) RasterYSize 圖像的高度(Y方向上的像素個(gè)數(shù))
Toggle line numbers
1 >>> dataset.RasterXSize
2 950
3 >>> dataset.RasterYSize
4 700
5 >>>
可以看出我們的圖像大小是950*700。還是很小的一張圖。 ReadRaster 讀取圖像數(shù)據(jù)(以二進(jìn)制的形式) ReadAsArray 讀取圖像數(shù)據(jù)(以數(shù)組的形式)
Toggle line numbers
1 >>> help(dataset.ReadRaster)
2 Help on method ReadRaster in module gdal:
3 ReadRaster(self, xoff, yoff, xsize, ysize, buf_xsize=None, buf_ysize=None, buf_t
4 ype=None, band_list=None) method of gdal.Dataset instance
5 >>> help(dataset.ReadAsArray)
6 Help on method ReadAsArray in module gdal:
7 ReadAsArray(self, xoff=0, yoff=0, xsize=None, ysize=None) method of gdal.Dataset
8 instance
9 >>>
這兩個(gè)函數(shù)很重要，它們直接讀取圖像的數(shù)據(jù)，可以看到兩個(gè)函數(shù)的幫助中有一大溜的參數(shù)。解釋一下： xoff,yoff,xsize,ysize 你可能不想讀取整張圖像。只想讀取其中的一部分。那么就用xoff，yoff指定想要讀取的部分原點(diǎn)位置在整張圖像中距離全圖原點(diǎn)的位置。用xsize和ysize指定要讀取部分圖像的矩形大小。 buf_xsize buf_ysize 你可以在讀取出一部分圖像后進(jìn)行縮放。那么就用這兩個(gè)參數(shù)來(lái)定義縮放后圖像最終的寬和高，gdal將幫你縮放到這個(gè)大小。 buf_type 如果你要讀取的圖像的數(shù)據(jù)類型不是你想要的(比如原圖數(shù)據(jù)類型是short，你要把它們縮小成byte)，就可以設(shè)置它。 band_list 這就適應(yīng)上面多波段的情況。你可以指定讀取的波段序列。要哪幾個(gè)波段，不要哪幾個(gè)波段，你說(shuō)了算。 舉個(gè)例子吧：
Toggle line numbers
1 >>> dataset.ReadAsArray(230,270,10,10)
2 array([[255, 255, 255, 232, 232, 255, 255, 255, 255, 222],
3 [255, 255, 255, 255, 255, 255, 210, 110, 11, 122],
4 [255, 255, 255, 255, 255, 255, 210, 255, 11, 243],
5 [201, 255, 255, 255, 255, 200, 200, 110, 122, 243],
6 [111, 211, 255, 201, 255, 255, 100, 11, 132, 243],
7 [255, 100, 100, 100, 110, 100, 110, 111, 122, 243],
8 [255, 255, 255, 255, 255, 255, 122, 222, 255, 255],
9 [255, 255, 255, 255, 255, 255, 243, 243, 255, 255],
10 [255, 255, 255, 255, 255, 255, 255, 255, 255, 255],
11 [255, 255, 255, 255, 255, 255, 255, 255, 255, 255]],‘b‘)
12 >>> dataset.ReadRaster(230,270,10,10)
13 ‘\xff\xff\xff\xe8\xe8\xff\xff\xff\xff\xde\xff\xff\xff\xff\xff\xff\xd2n\x0bz\xff\
14 xff\xff\xff\xff\xff\xd2\xff\x0b\xf3\xc9\xff\xff\xff\xff\xc8\xc8nz\xf3o\xd3\xff\x
15 c9\xff\xffd\x0b\x84\xf3\xffdddndnoz\xf3\xff\xff\xff\xff\xff\xffz\xde\xff\xff\xff
16 \xff\xff\xff\xff\xff\xf3\xf3\xff\xff\xff\xff\xff\xff\xff\xff\xff\xff\xff\xff\xff
17 \xff\xff\xff\xff\xff\xff\xff\xff\xff‘
18 >>>
我們就把圖像中位于230,270，寬度10高度10的數(shù)據(jù)讀取出來(lái)了。 我們看完了數(shù)據(jù)集的主要函數(shù)。似乎已經(jīng)夠用了。的確，如果只是為了顯示圖像，這些的確已經(jīng)夠了。但是如果需要更多信息，我們就不得不進(jìn)入波段操作數(shù)據(jù)(實(shí)際上我們大多數(shù)時(shí)候都需要進(jìn)入band獲取信息)。下面我們現(xiàn)在來(lái)看看剛才讀取出來(lái)的那個(gè)band有些什么東西可以供我們操作的(具體的API列表看這里)。
Toggle line numbers
1 >>> dir(band)
2 [‘AdviseRead‘, ‘Checksum‘, ‘ComputeBandStats‘, ‘ComputeRasterMinMax‘, ‘DataType‘
3 , ‘Fill‘, ‘FlushCache‘, ‘GetDefaultHistogram‘, ‘GetDescription‘, ‘GetHistogram‘,
4 ‘GetMaximum‘, ‘GetMetadata‘, ‘GetMinimum‘, ‘GetNoDataValue‘, ‘GetOffset‘, ‘GetO
5 verview‘, ‘GetOverviewCount‘, ‘GetRasterColorInterpretation‘, ‘GetRasterColorTab
6 le‘, ‘GetScale‘, ‘GetStatistics‘, ‘ReadAsArray‘, ‘ReadRaster‘, ‘SetDefaultHistog
7 ram‘, ‘SetDescription‘, ‘SetMetadata‘, ‘SetNoDataValue‘, ‘SetRasterColorInterpre
8 tation‘, ‘SetRasterColorTable‘, ‘WriteArray‘, ‘WriteRaster‘, ‘XSize‘, ‘YSize‘, ‘
9 __doc__‘, ‘__init__‘, ‘__module__‘, ‘_o‘]
10 >>>
挑幾個(gè)有用的吧。
Toggle line numbers
1 >>> band.XSize
2 950
3 >>> band.YSize
4 700
5 >>> band.DataType
6 1
7 >>>
不用解釋了吧，波段圖像的寬和高(象元為單位)。DataType，圖像中實(shí)際數(shù)值的數(shù)據(jù)類型。具體數(shù)據(jù)類型定義在gdalconst模塊里。使用的時(shí)候用import gdalconst引入。
Toggle line numbers
1 >>> import gdalconst
2 >>> dir(gdalconst)
3 [‘CE_Debug‘, ‘CE_Failure‘, ‘CE_Fatal‘, ‘CE_None‘, ‘CE_Warning‘, ‘CPLES_Backslash
4 Quotable‘, ‘CPLES_CSV‘, ‘CPLES_SQL‘, ‘CPLES_URL‘, ‘CPLES_XML‘, ‘CPLE_AppDefined‘
5 , ‘CPLE_AssertionFailed‘, ‘CPLE_FileIO‘, ‘CPLE_IllegalArg‘, ‘CPLE_NoWriteAccess‘
6 , ‘CPLE_None‘, ‘CPLE_NotSupported‘, ‘CPLE_OpenFailed‘, ‘CPLE_OutOfMemory‘, ‘CPLE
7 _UserInterrupt‘, ‘CXT_Attribute‘, ‘CXT_Comment‘, ‘CXT_Element‘, ‘CXT_Literal‘, ‘
8 CXT_Text‘, ‘DCAP_CREATE‘, ‘DCAP_CREATECOPY‘, ‘DMD_CREATIONDATATYPES‘, ‘DMD_CREAT
9 IONOPTIONLIST‘, ‘DMD_EXTENSION‘, ‘DMD_HELPTOPIC‘, ‘DMD_LONGNAME‘, ‘DMD_MIMETYPE‘
10 , ‘GA_ReadOnly‘, ‘GA_Update‘, ‘GCI_AlphaBand‘, ‘GCI_BlackBand‘, ‘GCI_BlueBand‘,
11 ‘GCI_CyanBand‘, ‘GCI_GrayIndex‘, ‘GCI_GreenBand‘, ‘GCI_HueBand‘, ‘GCI_LightnessB
12 and‘, ‘GCI_MagentaBand‘, ‘GCI_PaletteIndex‘, ‘GCI_RedBand‘, ‘GCI_SaturationBand‘
13 , ‘GCI_Undefined‘, ‘GCI_YellowBand‘, ‘GDT_Byte‘, ‘GDT_CFloat32‘, ‘GDT_CFloat64‘,
14 ‘GDT_CInt16‘, ‘GDT_CInt32‘, ‘GDT_Float32‘, ‘GDT_Float64‘, ‘GDT_Int16‘, ‘GDT_Int
15 32‘, ‘GDT_TypeCount‘, ‘GDT_UInt16‘, ‘GDT_UInt32‘, ‘GDT_Unknown‘, ‘GF_Read‘, ‘GF_
16 Write‘, ‘GPI_CMYK‘, ‘GPI_Gray‘, ‘GPI_HLS‘, ‘GPI_RGB‘, ‘GRA_Bilinear‘, ‘GRA_Cubic
17 ‘, ‘GRA_CubicSpline‘, ‘GRA_NearestNeighbour‘, ‘__builtins__‘, ‘__doc__‘, ‘__file
18 __‘, ‘__name__‘]
19 >>>
那些GDT開(kāi)頭的就是數(shù)值數(shù)據(jù)類型。
Toggle line numbers
1 >>> band.GetNoDataValue()
2 65535.0
3 >>> band.GetMaximum()
4 >>> band.GetMinimum()
5 >>> band.ComputeRasterMinMax()
6 (1.0, 255.0)
7 >>>
Maximum 是表示在本波段數(shù)值中最大的值，Minimum當(dāng)然就是表示本波段中最小的值啦。我們可以看到在一開(kāi)始這兩個(gè)都沒(méi)有值。因?yàn)閷?duì)于文件格式不會(huì)有固有的最大最小值。所以我們通過(guò)函數(shù)ComputeRasterMinMax計(jì)算得到了。注意！這里的最大最小值不包括“無(wú)意義值”！也就是上面顯示的NoDataValue。需要解釋一下“無(wú)意義值”。不要以為0或者255在任何情況下都無(wú)意義。在很多情況下0，255需要和其他值一樣表示一個(gè)實(shí)際意義。雖然可能它最終會(huì)被顯示得和黑色一樣。而一些位置上的點(diǎn)要表示的意思是“什么也不是”，它在那個(gè)位置上只是為了占一個(gè)位置，使得整體圖像看起來(lái)像個(gè)矩形而已。在做實(shí)際應(yīng)用的時(shí)候兩種值的處理將會(huì)完全不一樣。所以需要設(shè)置無(wú)意義值，來(lái)和其他的值區(qū)別開(kāi)來(lái)。而用ComputeRasterMinMax算出的最大最小值，是排除了無(wú)意義值后計(jì)算出來(lái)的最大最小值。
Toggle line numbers
1 >>> band.GetRasterColorInterpretation()
2 2
3 >>> gdalconst.GCI_PaletteIndex
4 2
5 >>> colormap = band.GetRasterColorTable()
6 >>> dir(colormap)
7 [‘Clone‘, ‘GetColorEntry‘, ‘GetColorEntryAsRGB‘, ‘GetCount‘, ‘GetPaletteInterpre
8 tation‘, ‘SetColorEntry‘, ‘__del__‘, ‘__doc__‘, ‘__init__‘, ‘__module__‘, ‘__str
9 __‘, ‘_o‘, ‘own_o‘, ‘serialize‘]
10 >>> colormap.GetCount()
11 256
12 >>> colormap.GetPaletteInterpretation()
13 1
14 >>> gdalconst.GPI_RGB
15 1
16 >>> for i in range(colormap.GetCount()):
17 ... print colormap.GetColorEntry(i),
18 ...
19 (0, 0, 0, 255) (0, 0, 28, 255) (0, 0, 56, 255) (0, 0, 85, 255) (0, 0, 113, 255)
20 (0, 0, 142, 255) (0, 0, 170, 255) (0, 0, 199, 255) (0, 0, 227, 255) (0, 0, 255,
21 255) (0, 28, 0, 255) (0, 28, 28, 255) (0, 28, 56, 255) (0, 28, 85, 255) (0, 28,
22 113, 255) (0, 28, 142, 255) (0, 28, 170, 255) (0, 28, 199, 255) (0, 28, 227, 255
23 ) (0, 28, 255, 255) (0, 56, 0, 255) (0, 56, 28, 255) (0, 56, 56, 255) (0, 56, 85
24 , 255) (0, 56, 113, 255) (0, 56, 142, 255) (0, 56, 170, 255) (0, 56, 199, 255) (
25 0, 56, 227, 255) (0, 56, 255, 255) (0, 85, 0, 255) (0, 85, 28, 255) (0, 85, 56,
26 255) (0, 85, 85, 255) (0, 85, 113, 255) (0, 85, 142, 255) (0, 85, 170, 255) (0,
27 85, 199, 255) (0, 85, 227, 255) (0, 85, 255, 255) (0, 113, 0, 255) (0, 113, 28,
28 255) (0, 113, 56, 255) (0, 113, 85, 255) (0, 113, 113, 255) (0, 113, 142, 255) (
29 0, 113, 170, 255) (0, 113, 199, 255) (0, 113, 227, 255) (0, 113, 255, 255) (0, 1
30 42, 0, 255) (0, 142, 28, 255) (0, 142, 56, 255) (0, 142, 85, 255) (0, 142, 113,
31 255) (0, 142, 142, 255) (0, 142, 170, 255) (0, 142, 199, 255) (0, 142, 227, 255)
32 (0, 142, 255, 255) (0, 170, 0, 255) (0, 170, 28, 255) (0, 170, 56, 255) (0, 170
33 , 85, 255) (0, 170, 113, 255) (0, 170, 142, 255) (0, 170, 170, 255) (0, 170, 199
34 , 255) (0, 170, 227, 255) (0, 170, 255, 255) (0, 199, 0, 255) (0, 199, 28, 255)
35 (0, 199, 56, 255) (0, 199, 85, 255) (0, 199, 113, 255) (0, 199, 142, 255) (0, 19
36 9, 170, 255) (0, 199, 199, 255) (0, 199, 227, 255) (0, 199, 255, 255) (0, 227, 0
37 , 255) (0, 227, 28, 255) (0, 227, 56, 255) (0, 227, 85, 255) (0, 227, 113, 255)
38 (0, 227, 142, 255) (0, 227, 170, 255) (0, 227, 199, 255) (0, 227, 227, 255) (0,
39 227, 255, 255) (0, 255, 0, 255) (0, 255, 28, 255) (0, 255, 56, 255) (0, 255, 85,
40 255) (0, 255, 113, 255) (0, 255, 142, 255) (0, 255, 170, 255) (0, 255, 199, 255
41 ) (0, 255, 227, 255) (0, 255, 255, 255) (28, 0, 0, 255) (28, 0, 28, 255) (28, 0,
42 56, 255) (28, 0, 85, 255) (28, 0, 113, 255) (28, 0, 142, 255) (28, 0, 170, 255)
43 (28, 0, 199, 255) (28, 0, 227, 255) (28, 0, 255, 255) (28, 28, 0, 255) (28, 28,
44 28, 255) (28, 28, 56, 255) (28, 28, 85, 255) (28, 28, 113, 255) (28, 28, 142, 2
45 55) (28, 28, 170, 255) (28, 28, 199, 255) (28, 28, 227, 255) (28, 28, 255, 255)
46 (28, 56, 0, 255) (28, 56, 28, 255) (28, 56, 56, 255) (28, 56, 85, 255) (28, 56,
47 113, 255) (28, 56, 142, 255) (28, 56, 170, 255) (28, 56, 199, 255) (28, 56, 227,
48 255) (28, 56, 255, 255) (28, 85, 0, 255) (28, 85, 28, 255) (28, 85, 56, 255) (2
49 8, 85, 85, 255) (28, 85, 113, 255) (28, 85, 142, 255) (28, 85, 170, 255) (28, 85
50 , 199, 255) (28, 85, 227, 255) (28, 85, 255, 255) (28, 113, 0, 255) (28, 113, 28
51 , 255) (28, 113, 56, 255) (28, 113, 85, 255) (28, 113, 113, 255) (28, 113, 142,
52 255) (28, 113, 170, 255) (28, 113, 199, 255) (28, 113, 227, 255) (28, 113, 255,
53 255) (28, 142, 0, 255) (28, 142, 28, 255) (28, 142, 56, 255) (28, 142, 85, 255)
54 (28, 142, 113, 255) (28, 142, 142, 255) (28, 142, 170, 255) (28, 142, 199, 255)
55 (28, 142, 227, 255) (28, 142, 255, 255) (28, 170, 0, 255) (28, 170, 28, 255) (28
56 , 170, 56, 255) (28, 170, 85, 255) (28, 170, 113, 255) (28, 170, 142, 255) (28,
57 170, 170, 255) (28, 170, 199, 255) (28, 170, 227, 255) (28, 170, 255, 255) (28,
58 199, 0, 255) (28, 199, 28, 255) (28, 199, 56, 255) (28, 199, 85, 255) (28, 199,
59 113, 255) (28, 199, 142, 255) (28, 199, 170, 255) (28, 199, 199, 255) (28, 199,
60 227, 255) (28, 199, 255, 255) (28, 227, 0, 255) (28, 227, 28, 255) (28, 227, 56,
61 255) (28, 227, 85, 255) (28, 227, 113, 255) (28, 227, 142, 255) (28, 227, 170,
62 255) (28, 227, 199, 255) (28, 227, 227, 255) (28, 227, 255, 255) (28, 255, 0, 25
63 5) (28, 255, 28, 255) (28, 255, 56, 255) (28, 255, 85, 255) (28, 255, 113, 255)
64 (28, 255, 142, 255) (28, 255, 170, 255) (28, 255, 199, 255) (28, 255, 227, 255)
65 (28, 255, 255, 255) (56, 0, 0, 255) (56, 0, 28, 255) (56, 0, 56, 255) (56, 0, 85
66 , 255) (56, 0, 113, 255) (56, 0, 142, 255) (56, 0, 170, 255) (56, 0, 199, 255) (
67 56, 0, 227, 255) (56, 0, 255, 255) (56, 28, 0, 255) (56, 28, 28, 255) (56, 28, 5
68 6, 255) (56, 28, 85, 255) (56, 28, 113, 255) (56, 28, 142, 255) (56, 28, 170, 25
69 5) (56, 28, 199, 255) (56, 28, 227, 255) (56, 28, 255, 255) (56, 56, 0, 255) (56
70 , 56, 28, 255) (56, 56, 56, 255) (56, 56, 85, 255) (56, 56, 113, 255) (56, 56, 1
71 42, 255) (56, 56, 170, 255) (56, 56, 199, 255) (56, 56, 227, 255) (56, 56, 255,
72 255) (56, 85, 0, 255) (56, 85, 28, 255) (56, 85, 56, 255) (56, 85, 85, 255) (56,
73 85, 113, 255) (56, 85, 142, 255) (56, 85, 170, 255) (56, 85, 199, 255) (56, 85,
74 227, 255) (56, 85, 255, 255) (56, 113, 0, 255) (56, 113, 28, 255) (56, 113, 56,
75 255) (56, 113, 85, 255) (56, 113, 113, 255) (56, 113, 142, 255) (56, 113, 170,
76 255) (56, 113, 199, 255) (56, 113, 227, 255) (56, 113, 255, 255) (56, 142, 0, 25
77 5) (56, 142, 28, 255) (56, 142, 56, 255) (56, 142, 85, 255) (56, 142, 113, 255)
78 (56, 142, 142, 255)
79 >>>
通過(guò)GetRasterColorInterpretation，我們知道我們的圖像是一個(gè)顏色表索引的圖像而不是純粹的黑白灰度圖像(PaletteIndex，其他的顏色模型，可以察看gdalconst模塊中GCI打頭的枚舉值)。這意味著我們讀出的數(shù)據(jù)有可能不是真實(shí)的數(shù)據(jù)。這些數(shù)據(jù)只是一個(gè)個(gè)實(shí)際數(shù)據(jù)的索引。真實(shí)數(shù)據(jù)存儲(chǔ)在另一個(gè)表中。我們通過(guò)ReadRaster讀出的數(shù)據(jù)值只是對(duì)應(yīng)到這個(gè)表的一個(gè)索引而已。我們需要通過(guò)讀出這些數(shù)據(jù)，并在真實(shí)數(shù)據(jù)表中找出真實(shí)數(shù)據(jù)，重新組織成一個(gè)RGB表才能用來(lái)繪制。如果我們不經(jīng)過(guò)對(duì)應(yīng)，我們繪制出來(lái)的東西可能什么東西都不是。 用GetRasterColorTable獲得了顏色表，通過(guò)GetPaletteInterpretation我們知道我們獲得的顏色表是一個(gè)RGB顏色表。GDAL支持多種顏色表，具體可以參考gdalconst模塊中GPI打頭的枚舉值。然后我們可以通過(guò)GetCount獲得顏色的數(shù)量。通過(guò)GetColorEntry獲得顏色表中的值。這里的顏色值都是一個(gè)4值的元組。里面有意義的只有前三個(gè)(如果顏色模型是GPI_RGB, GPI_HLS，都使用前3個(gè)，如果采用GPI_CMYK，則4個(gè)值都有意義了)。
Toggle line numbers
1 >>> help(band.ReadAsArray)
2 Help on method ReadAsArray in module gdal:
3 ReadAsArray(self, xoff=0, yoff=0, win_xsize=None, win_ysize=None, buf_xsize=None
4 , buf_ysize=None, buf_obj=None) method of gdal.Band instance
5 >>> help(band.ReadRaster)
6 Help on method ReadRaster in module gdal:
7 ReadRaster(self, xoff, yoff, xsize, ysize, buf_xsize=None, buf_ysize=None, buf_t
8 ype=None) method of gdal.Band instance
9 >>>
顯然，band里的ReadAsArray參數(shù)顯然比dataset里面的要好用，而ReadRaster則差不多。但是ReadAsArray讀出的是數(shù)組，可以用Numeric模塊進(jìn)行矩陣魔法。ReadRaster讀出的是二進(jìn)制，雖然可以直接繪制，但是對(duì)于一些繪圖API來(lái)說(shuō)，對(duì)[[RRR...][GGG...][BBB...]]表的處理明顯不如[[RGB][RGB]...]，有的甚至不支持。雖然可以用struct.unpack來(lái)拆封，可效率上就差很多(而且拆封出來(lái)還是數(shù)組)。數(shù)組在矩陣魔法的控制之下則會(huì)顯得十分方便快捷，最后用tostring直接轉(zhuǎn)化稱為二進(jìn)制繪制，速度也相當(dāng)快。 好了，快速開(kāi)始已經(jīng)使我們可以初步看清楚了gdal中圖像的組織。下面用一句話總結(jié)一下：波段組成圖像，波段指揮顏色。
4. 反饋如果您發(fā)現(xiàn)我寫的東西中有問(wèn)題，或者您對(duì)我寫的東西有意見(jiàn)，請(qǐng)一定要發(fā)郵件跟我講，Email( linux_23@163.com ) |
|
|
is a translator library for raster geospatial data formats that is released under an
），下面的$(HDF4_DIR)可以不用改，那個(gè)是變量，會(huì)自動(dòng)替代HDF4_DIR 路徑。
