了解GDAL圖像處理(轉)


原文:http://blog.163.com/chenxingfeng_001/blog/static/10627072006115836560/

1. GDAL庫介紹
簡單地說,GDAL是一個操作各種柵格地理數據格式的庫。包括讀取、寫入、轉換、處理各種柵格數據格式(有些特定的格式對一些操作如寫入等不支持)。它使用了一個單一的抽象數據模型就支持了大多數的柵格數據(GIS對柵格,矢量,3D數據模型的抽象能力實在令人嘆服)。當然除了柵格操作,這個庫還同時包括了操作矢量數據的另一個有名的庫ogr(ogr這個庫另外介紹),這樣這個庫就同時具備了操作柵格和矢量數據的能力,買一送一,這么合算的買賣為什么不做。最最最重要的是這個庫是跨平台的,開源的!如今這個庫對各種數據格式的支持強大到令人嘖嘖的地步了。如果你對他的強大有什么懷疑的話,看看這里一大串的GDAL所支持格式清單,嚇到了吧!再看看它的主頁最后那些使用了它作為底層數據處理的軟件列表吧!其中你可以不知道GRASS,你也可以不知道Quantum GIS (QGIS),但是你總該知道Google Earth吧!
你即使不玩GIS,這個庫也是滿有用的。首先,哪個庫支持這么多柵格(圖片)格式,哪個庫在C/C++/python/ruby/VB/java/C#(這個暫時不完全支持)下都能用,而且都一樣用?退一步講,3S軟件又不一定要用在3S下(很多醫學影像就是用PCI軟件來處理的)。再退一步,你的生活即使和3S一點關系都沒有,柵格數據又不單單只有GIS下才用到。你大可用這個庫來讀取jpg,gif,tif,xpm等格式。而且對各種格式支持得不是一般的好,很大一部分非標准格式照樣支持得非常好。我曾經在java下玩過jai,以及一系列jai的擴展庫,一些圖像格式在很多圖片瀏覽器中都可以正確讀取(有的甚至不是非標准格式),用jai死活就讀不出來!
這個庫的python版和其他的python庫結合的很好。最直接、明顯的支持是使用Numeric庫來進行數據讀取和操作。各種矩陣魔術可以發揮得淋漓盡致(圖像其實就是矩陣)。而且按我的觀點,python對矩陣的操作比其他的語言有明顯的優勢。寫出來的東西比其他語言寫出來的短小的多,而且好看得多。並且python的弱類型在處理柵格數據格式類型的時候代碼量比強類型的語言少了數倍(不用double,byte,short等等分開處理,這簡直就是先天上的優勢)。所以我就喜歡用python做圖像的處理。所以就連GIS界的微軟ESRI也直接在ARCGIS9中用python來作柵格數據的導入導出。一句話,真是太方便啦!
2. 安裝
2.1. windows下的安裝
官方安裝文檔在這里。下面是我自己的實踐步驟:
先去http://www.gdal.org/dl/下一個版本,解壓。 打開控制台,輸入: “D:\Program Files\Microsoft Visual Studio .NET 2003\Vc7\bin\vcvars32.bat" 注冊vc的編譯環境。
打 開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的安裝路徑。 下面的參數愛改什么就把前面的#刪除(要看您有沒有那些庫的源碼),注意一下路徑就可以了。我是都沒改。 后面就依次運行
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以上版本不能順利編譯,有一個地方指針轉換出錯。可能是2005的編譯器比以往的嚴厲一點吧。另外,安裝了QGIS,對編譯也有一些影響,主要是proj庫的沖突,導致一個找不到"d:/program.obj"文件的錯誤,如果你有靜態編譯過proj,那么你可以打開nmake.opt修改有關proj的設置,如果搞不定,就卸載QGIS,然后編譯,編譯后再安裝QGIS.呵呵,還好QGIS的體積沒有ArcGIS那么可怕.
2.2. linux下的安裝
linux下的安裝就更簡單了。直接
Toggle line numbers
1 ./configure
2 make
3 su
4 make install
5 ldconfig

就ok(默認就已經支持python)。當然在第一步的時候需要看看是否依賴的庫都安裝了。如果缺少,就去安裝一個。如果對configure的條件不理解,就用./configure --help看看具體情況。
2.3. 安裝其他驅動
這里講一個安裝hdf4的驅動的例子(默認情況下gdal是不安裝hdf4的),其他驅動應該和這個也差不了多少吧,可以作為其他的參考。完整步驟如下:
在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_DIR)\lib\hd421.lib $(HDF4_DIR)\lib\hm421.lib



用HDF4_LIB=/LIBPATH:這種形式似乎可以建立gdal的庫,但是往下編譯會出錯。而且要把$(HDF4_DIR)\dll和$(HDF4_DIR)\lib拷貝到同一個目錄下,不然會提示找不到庫
你也可以試一試在D:\Program Files\Microsoft Visual Studio .NET 2003\Common7\Tools\vsvars32.bat文件中添加HDF4_LIB路徑到“@set LIB=”這行的末尾(不要忘記;的分割符)。
然后找一下"INC="這行,把 -I$(HDF4_DIR)\include加到下一行的末尾(應該也可以在vsvars32.bat中添加路徑,不過要重啟命令行)。
然后編譯吧!祝你好運。
注意:上面的HDF4_DIR 是我本機的路徑,你要根據你自己的路徑進行設置(想起我的一個老師說過的話:“抄人家的作業可以,不要連名字也一起抄走啊 ” ),下面的$(HDF4_DIR)可以不用改,那個是變量,會自動替代HDF4_DIR 路徑。
編譯成功后,要HDF4能運行,還需要兩個庫,一個是zlib,一個是szip,可以到下面兩個鏈接去下載一個
ftp://hdf.ncsa.uiuc.edu/lib-external/zlib/1.2/bin
ftp://hdf.ncsa.uiuc.edu/lib-external/szip/2.0/bin/windows
把這兩個庫下載后解壓,然后設置PATH系統變量,使得它們在默認狀態下也可以被動態鏈接成功 。
在Linux下比在Windows下簡單:
只要用./configure --help察看一下打開HDF4的編譯開關(包括庫路徑,頭文件路徑等,看清楚),然后在./configure 后面加上那個開關以及hdf4的安裝路徑后就可以了。在configure后gdal會提示是否支持HDF4。 編譯后也要把zlib和szip的動態鏈接庫設置好 。
到此你已經可以用C/C++來操作gdal讀寫hdf4的格式了!
最后,為了讓Python能夠支持hdf的讀寫,別忘了把重新生成安裝后的pymod目錄下的內容,還有gdal13,還有那兩個hdf的庫,還有zlib,szip的庫拷貝到Python的Lib\site-packages目錄下 。
2.4. 下載
如果你實在玩不轉,可以在這里下載已經編譯好的gdal1.3.2程序庫 以及其依賴的其他庫,其中包括hdf4,hdf5支持,以及proj,geos插件。注意,這里的geos是靜態鏈接的,注意版權(geos是LGPL的license)。hdf4和hdf5用的是release版本。這里是我的nmake配置文件,你可以對照你的實際情況參考一下。
3. 快速開始
其實在主站的教程里已經有python的示例了。但是我們還是按照自己的思路來開始吧。
第一步就是打開一個數據集。對於“數據集”這個名詞大家可能不會太習慣,但是對於一般的格式來說,一個“數據集”就是一個文件,比如一個gif文件就是一個以gif為擴展名的文件。但是對於眾多RS數據來說,一個數據集包含的絕對不僅僅是一個文件。對於很多RS數據,他們把一張圖像分成數個圖像文件,然后放在一個文件夾中,用一些額外的文件來組織它們之間的關系,形成一個“數據集”。如果你不理解,那么就算了,當成jpg或者gif文件好了。
下面我們打開一個tiff文件(GeoTIFF)。這個文件是我從GRASS的示例數據spearfish中導出的一個同名影像數據。
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 >>>

這樣我們就打開了這個文件。並且我們可以看到可以供我們調用的函數們(更具體的API列表可以看這里)。現在我們不做修改,不做添加,所以只要帶有Set開頭的函數以及有Write開頭的函數我們暫時都不管。因為RS影像必然要和地理上的位置掛上鈎,才能把圖像正確鋪展到一個坐標系中。其中的信息和對應關系有點復雜,不適合在快速開始中介紹,我們暫時也先不管。這里需要注意的就是幾個函數。
GetDescription 獲得柵格的描述信息。
Toggle line numbers
1 >>> dataset.GetDescription()
2 'j:/gisdata/gtif/spot.tif'
3 >>>

看來這里的圖像描述是圖像的路徑名,但是這是和各種不同數據集相關的,不同數據集可能有不同的描述。這要看讀取驅動的實現作者的高興了。
RasterCount 獲得柵格數據集的波段數。
GetRasterBand 獲得柵格數據集的波段。
Toggle line numbers
1 >>> dataset.RasterCount
2 1
3 >>> band = dataset.GetRasterBand(1)
4 >>>

Band這個詞可以翻譯成“波段。
這里我們獲取了第一個波段(紅色值組成的表)。注意!這里的波段獲取和通常的C數組獲取不一樣,開始是1不是0。獲取了波段,我們就可以在下面的操作中讀取這個波段的所有數值。
RasterXSize 圖像的寬度(X方向上的像素個數)
RasterYSize 圖像的高度(Y方向上的像素個數)
Toggle line numbers
1 >>> dataset.RasterXSize
2 950
3 >>> dataset.RasterYSize
4 700
5 >>>

可以看出我們的圖像大小是950*700。還是很小的一張圖。
ReadRaster 讀取圖像數據(以二進制的形式)
ReadAsArray 讀取圖像數據(以數組的形式)
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 >>>

這兩個函數很重要,它們直接讀取圖像的數據,可以看到兩個函數的幫助中有一大溜的參數。解釋一下: xoff,yoff,xsize,ysize 你可能不想讀取整張圖像。只想讀取其中的一部分。那么就用xoff,yoff指定想要讀取的部分原點位置在整張圖像中距離全圖原點的位置。用xsize和ysize指定要讀取部分圖像的矩形大小。
buf_xsize buf_ysize 你可以在讀取出一部分圖像后進行縮放。那么就用這兩個參數來定義縮放后圖像最終的寬和高,gdal將幫你縮放到這個大小。
buf_type 如果你要讀取的圖像的數據類型不是你想要的(比如原圖數據類型是short,你要把它們縮小成byte),就可以設置它。
band_list 這就適應上面多波段的情況。你可以指定讀取的波段序列。要哪幾個波段,不要哪幾個波段,你說了算。
舉個例子吧:
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的數據讀取出來了。
我們看完了數據集的主要函數。似乎已經夠用了。的確,如果只是為了顯示圖像,這些的確已經夠了。但是如果需要更多信息,我們就不得不進入波段操作數據(實際上我們大多數時候都需要進入band獲取信息)。下面我們現在來看看剛才讀取出來的那個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 >>>

挑幾個有用的吧。
Toggle line numbers
1 >>> band.XSize
2 950
3 >>> band.YSize
4 700
5 >>> band.DataType
6 1
7 >>>

不用解釋了吧,波段圖像的寬和高(象元為單位)。DataType,圖像中實際數值的數據類型。具體數據類型定義在gdalconst模塊里。使用的時候用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開頭的就是數值數據類型。
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 是表示在本波段數值中最大的值,Minimum當然就是表示本波段中最小的值啦。我們可以看到在一開始這兩個都沒有值。因為對於文件格式不會有固有的最大最小值。所以我們通過函數ComputeRasterMinMax計算得到了。注意!這里的最大最小值不包括“無意義值”!也就是上面顯示的NoDataValue。需要解釋一下“無意義值”。不要以為0或者255在任何情況下都無意義。在很多情況下0,255需要和其他值一樣表示一個實際意義。雖然可能它最終會被顯示得和黑色一樣。而一些位置上的點要表示的意思是“什么也不是”,它在那個位置上只是為了占一個位置,使得整體圖像看起來像個矩形而已。在做實際應用的時候兩種值的處理將會完全不一樣。所以需要設置無意義值,來和其他的值區別開來。而用ComputeRasterMinMax算出的最大最小值,是排除了無意義值后計算出來的最大最小值。
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>>>

通過GetRasterColorInterpretation,我們知道我們的圖像是一個顏色表索引的圖像而不是純粹的黑白灰度圖像(PaletteIndex,其他的顏色模型,可以察看gdalconst模塊中GCI打頭的枚舉值)。這意味着我們讀出的數據有可能不是真實的數據。這些數據只是一個個實際數據的索引。真實數據存儲在另一個表中。我們通過ReadRaster讀出的數據值只是對應到這個表的一個索引而已。我們需要通過讀出這些數據,並在真實數據表中找出真實數據,重新組織成一個RGB表才能用來繪制。如果我們不經過對應,我們繪制出來的東西可能什么東西都不是。
用GetRasterColorTable獲得了顏色表,通過GetPaletteInterpretation我們知道我們獲得的顏色表是一個RGB顏色表。GDAL支持多種顏色表,具體可以參考gdalconst模塊中GPI打頭的枚舉值。然后我們可以通過GetCount獲得顏色的數量。通過GetColorEntry獲得顏色表中的值。這里的顏色值都是一個4值的元組。里面有意義的只有前三個(如果顏色模型是GPI_RGB, GPI_HLS,都使用前3個,如果采用GPI_CMYK,則4個值都有意義了)。
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參數顯然比dataset里面的要好用,而ReadRaster則差不多。但是ReadAsArray讀出的是數組,可以用Numeric模塊進行矩陣魔法。ReadRaster讀出的是二進制,雖然可以直接繪制,但是對於一些繪圖API來說,對[[RRR...][GGG...][BBB...]]表的處理明顯不如[[RGB][RGB]...],有的甚至不支持。雖然可以用struct.unpack來拆封,可效率上就差很多(而且拆封出來還是數組)。數組在矩陣魔法的控制之下則會顯得十分方便快捷,最后用tostring直接轉化稱為二進制繪制,速度也相當快。
好了,快速開始已經使我們可以初步看清楚了gdal中圖像的組織。下面用一句話總結一下:波段組成圖像,波段指揮顏色。


免責聲明!

本站轉載的文章為個人學習借鑒使用,本站對版權不負任何法律責任。如果侵犯了您的隱私權益,請聯系本站郵箱yoyou2525@163.com刪除。



 
粵ICP備18138465號   © 2018-2025 CODEPRJ.COM