1. PythonGIS
  2. GDAL
  3. 对真彩色图像与索引图像的处理

对真彩色图像与索引图像的处理

图像 lu75i1.tif 是索引图像，使用gdalinfo来查看的话，可以看见有如下的信息：

$ gdalinfo lu75i1.tif

......

251: 80,249,240,255

252: 114,239,249,255

253: 65,254,253,255

......

前后都省略了。显示的是索引图像的颜色查找表。

In [1]:
from osgeo import gdal
from osgeo import gdalconst
dataset = gdal.Open('/gdata/lu75i1.tif')
band = dataset.GetRasterBand(1)
band.GetRasterColorInterpretation()
Out[1]:
2
In [2]:
import gdalconst
gdalconst.GCI_PaletteIndex
Out[2]:
2
In [3]:
colormap = band.GetRasterColorTable()
dir(colormap)
colormap.GetCount()
colormap.GetPaletteInterpretation()
gdalconst.GPI_RGB
for i in range(colormap.GetCount()):
    print("%i:%s" % (i, colormap.GetColorEntry(i)))

0:(196, 88, 130, 255)
1:(77, 25, 29, 255)
2:(5, 22, 252, 255)
3:(201, 170, 250, 255)
4:(26, 161, 35, 255)
5:(59, 116, 237, 255)
6:(13, 18, 89, 255)
7:(111, 238, 252, 255)
8:(247, 20, 221, 255)
9:(120, 30, 117, 255)
10:(40, 71, 82, 255)
11:(207, 245, 110, 255)
12:(64, 250, 35, 255)
13:(217, 195, 176, 255)
14:(28, 138, 99, 255)
15:(115, 68, 235, 255)
16:(237, 210, 2, 255)
17:(209, 10, 10, 255)
18:(128, 105, 32, 255)
19:(0, 0, 0, 255)
20:(0, 0, 0, 255)
21:(0, 0, 0, 255)
22:(0, 0, 0, 255)
23:(0, 0, 0, 255)
24:(0, 0, 0, 255)
25:(0, 0, 0, 255)
26:(0, 0, 0, 255)
27:(0, 0, 0, 255)
28:(0, 0, 0, 255)
29:(0, 0, 0, 255)
30:(0, 0, 0, 255)
31:(0, 0, 0, 255)
32:(0, 0, 0, 255)
33:(0, 0, 0, 255)
34:(0, 0, 0, 255)
35:(0, 0, 0, 255)
36:(0, 0, 0, 255)
37:(0, 0, 0, 255)
38:(0, 0, 0, 255)
39:(0, 0, 0, 255)
40:(0, 0, 0, 255)
41:(0, 0, 0, 255)
42:(0, 0, 0, 255)
43:(0, 0, 0, 255)
44:(0, 0, 0, 255)
45:(0, 0, 0, 255)
46:(0, 0, 0, 255)
47:(0, 0, 0, 255)
48:(0, 0, 0, 255)
49:(0, 0, 0, 255)
50:(0, 0, 0, 255)
51:(0, 0, 0, 255)
52:(0, 0, 0, 255)
53:(0, 0, 0, 255)
54:(0, 0, 0, 255)
55:(0, 0, 0, 255)
56:(0, 0, 0, 255)
57:(0, 0, 0, 255)
58:(0, 0, 0, 255)
59:(0, 0, 0, 255)
60:(0, 0, 0, 255)
61:(0, 0, 0, 255)
62:(0, 0, 0, 255)
63:(0, 0, 0, 255)
64:(0, 0, 0, 255)
65:(0, 0, 0, 255)
66:(0, 0, 0, 255)
67:(0, 0, 0, 255)
68:(0, 0, 0, 255)
69:(0, 0, 0, 255)
70:(0, 0, 0, 255)
71:(0, 0, 0, 255)
72:(0, 0, 0, 255)
73:(0, 0, 0, 255)
74:(0, 0, 0, 255)
75:(0, 0, 0, 255)
76:(0, 0, 0, 255)
77:(0, 0, 0, 255)
78:(0, 0, 0, 255)
79:(0, 0, 0, 255)
80:(0, 0, 0, 255)
81:(0, 0, 0, 255)
82:(0, 0, 0, 255)
83:(0, 0, 0, 255)
84:(0, 0, 0, 255)
85:(0, 0, 0, 255)
86:(0, 0, 0, 255)
87:(0, 0, 0, 255)
88:(0, 0, 0, 255)
89:(0, 0, 0, 255)
90:(0, 0, 0, 255)
91:(0, 0, 0, 255)
92:(0, 0, 0, 255)
93:(0, 0, 0, 255)
94:(0, 0, 0, 255)
95:(0, 0, 0, 255)
96:(0, 0, 0, 255)
97:(0, 0, 0, 255)
98:(0, 0, 0, 255)
99:(0, 0, 0, 255)
100:(0, 0, 0, 255)
101:(0, 0, 0, 255)
102:(0, 0, 0, 255)
103:(0, 0, 0, 255)
104:(0, 0, 0, 255)
105:(0, 0, 0, 255)
106:(0, 0, 0, 255)
107:(0, 0, 0, 255)
108:(0, 0, 0, 255)
109:(0, 0, 0, 255)
110:(0, 0, 0, 255)
111:(0, 0, 0, 255)
112:(0, 0, 0, 255)
113:(0, 0, 0, 255)
114:(0, 0, 0, 255)
115:(0, 0, 0, 255)
116:(0, 0, 0, 255)
117:(0, 0, 0, 255)
118:(0, 0, 0, 255)
119:(0, 0, 0, 255)
120:(0, 0, 0, 255)
121:(0, 0, 0, 255)
122:(0, 0, 0, 255)
123:(0, 0, 0, 255)
124:(0, 0, 0, 255)
125:(0, 0, 0, 255)
126:(0, 0, 0, 255)
127:(0, 0, 0, 255)
128:(0, 0, 0, 255)
129:(0, 0, 0, 255)
130:(0, 0, 0, 255)
131:(0, 0, 0, 255)
132:(0, 0, 0, 255)
133:(0, 0, 0, 255)
134:(0, 0, 0, 255)
135:(0, 0, 0, 255)
136:(0, 0, 0, 255)
137:(0, 0, 0, 255)
138:(0, 0, 0, 255)
139:(0, 0, 0, 255)
140:(0, 0, 0, 255)
141:(0, 0, 0, 255)
142:(0, 0, 0, 255)
143:(0, 0, 0, 255)
144:(0, 0, 0, 255)
145:(0, 0, 0, 255)
146:(0, 0, 0, 255)
147:(0, 0, 0, 255)
148:(0, 0, 0, 255)
149:(0, 0, 0, 255)
150:(0, 0, 0, 255)
151:(0, 0, 0, 255)
152:(0, 0, 0, 255)
153:(0, 0, 0, 255)
154:(0, 0, 0, 255)
155:(0, 0, 0, 255)
156:(0, 0, 0, 255)
157:(0, 0, 0, 255)
158:(0, 0, 0, 255)
159:(0, 0, 0, 255)
160:(0, 0, 0, 255)
161:(0, 0, 0, 255)
162:(0, 0, 0, 255)
163:(0, 0, 0, 255)
164:(0, 0, 0, 255)
165:(0, 0, 0, 255)
166:(0, 0, 0, 255)
167:(0, 0, 0, 255)
168:(0, 0, 0, 255)
169:(0, 0, 0, 255)
170:(0, 0, 0, 255)
171:(0, 0, 0, 255)
172:(0, 0, 0, 255)
173:(0, 0, 0, 255)
174:(0, 0, 0, 255)
175:(0, 0, 0, 255)
176:(0, 0, 0, 255)
177:(0, 0, 0, 255)
178:(0, 0, 0, 255)
179:(0, 0, 0, 255)
180:(0, 0, 0, 255)
181:(0, 0, 0, 255)
182:(0, 0, 0, 255)
183:(0, 0, 0, 255)
184:(0, 0, 0, 255)
185:(0, 0, 0, 255)
186:(0, 0, 0, 255)
187:(0, 0, 0, 255)
188:(0, 0, 0, 255)
189:(0, 0, 0, 255)
190:(0, 0, 0, 255)
191:(0, 0, 0, 255)
192:(0, 0, 0, 255)
193:(0, 0, 0, 255)
194:(0, 0, 0, 255)
195:(0, 0, 0, 255)
196:(0, 0, 0, 255)
197:(0, 0, 0, 255)
198:(0, 0, 0, 255)
199:(0, 0, 0, 255)
200:(0, 0, 0, 255)
201:(0, 0, 0, 255)
202:(0, 0, 0, 255)
203:(0, 0, 0, 255)
204:(0, 0, 0, 255)
205:(0, 0, 0, 255)
206:(0, 0, 0, 255)
207:(0, 0, 0, 255)
208:(0, 0, 0, 255)
209:(0, 0, 0, 255)
210:(0, 0, 0, 255)
211:(0, 0, 0, 255)
212:(0, 0, 0, 255)
213:(0, 0, 0, 255)
214:(0, 0, 0, 255)
215:(0, 0, 0, 255)
216:(0, 0, 0, 255)
217:(0, 0, 0, 255)
218:(0, 0, 0, 255)
219:(0, 0, 0, 255)
220:(0, 0, 0, 255)
221:(0, 0, 0, 255)
222:(0, 0, 0, 255)
223:(0, 0, 0, 255)
224:(0, 0, 0, 255)
225:(0, 0, 0, 255)
226:(0, 0, 0, 255)
227:(0, 0, 0, 255)
228:(0, 0, 0, 255)
229:(0, 0, 0, 255)
230:(0, 0, 0, 255)
231:(0, 0, 0, 255)
232:(0, 0, 0, 255)
233:(0, 0, 0, 255)
234:(0, 0, 0, 255)
235:(0, 0, 0, 255)
236:(0, 0, 0, 255)
237:(0, 0, 0, 255)
238:(0, 0, 0, 255)
239:(0, 0, 0, 255)
240:(0, 0, 0, 255)
241:(0, 0, 0, 255)
242:(0, 0, 0, 255)
243:(0, 0, 0, 255)
244:(0, 0, 0, 255)
245:(0, 0, 0, 255)
246:(0, 0, 0, 255)
247:(0, 0, 0, 255)
248:(0, 0, 0, 255)
249:(0, 0, 0, 255)
250:(0, 0, 0, 255)
251:(0, 0, 0, 255)
252:(0, 0, 0, 255)
253:(0, 0, 0, 255)
254:(0, 0, 0, 255)
255:(0, 0, 0, 0)

可以看到最后输出的几行，同样得到了这幅索引图像的颜色查找表。

通过GetRasterColorInterpretation()的返回值2， 我们知道我们的图像是一个颜色表索引的图像， 而不是纯粹的黑白灰度图像。这意味着我们读出的数据有可能不是真实的数据。 这些数据只是一个个实际数据的索引。真实数据存储在另一个表中。

ColorMap颜色定义

ColorMap的定义在下面有详细的解释 :

颜色表： 一个颜色表可能包含一个或者更多的如下面这种C结构描述的颜色元素。

typedef struct
{
/- gray, red, cyan or hue -/
short      c1;
/- green, magenta, or lightness -/
short      c2;
/- blue, yellow, or saturation -/
short      c3;
/- alpha or blackband -/
short      c4;
} GDALColorEntry;

颜色表也有一个色彩解析值。（GDALPaletteInterp）这个值有可能是下面值的一种。

并且描述了c1,c2,c3,c4该如何解释。

  • GPI_GRAY: c1表示灰度值
  • GPI_RGB: c1表示红，c2表示绿，c3表示蓝，c4表示Alpha
  • GPI_CYMP: c1表示青，c2表示品，c3表示黄，c4表示黑
  • GPI_HLS: c1表示色调，c2表示亮度，c3表示饱和度

虽然有颜色表示数的区别，但是用GetColorEntry读出的都是4个值， 根据PaletteInterp枚举值看截取其中的几个值形成的颜色。

In [4]:
from osgeo import gdal
dataset = gdal.Open("/gdata/lu75i1.tif")
band = dataset.GetRasterBand(1)
colormap = band.GetRasterColorTable()
colormap.GetPaletteInterpretation()
Out[4]:
1
In [5]:
colormap.GetCount()
Out[5]:
256

通过用GetRasterColorTable获得了颜色表， 通过GetPaletteInterpretation知道我们获得的颜色表是一个RGB颜色表。 GDAL支持多种颜色表，具体可以参考gdalconst模块中GPI打头的枚举值。 然后我们可以通过GetCount来获得颜色数量， 通过GetColorEntry获得颜色表中的值。这里的颜色值都是一个4值的元组。 里面有意义的只有前三个(如果颜色模型是GPI_RGB, GPI_HLS， 都使用前3个，如果采用GPI_CMYK，则4个值都有意义了)。

访问索引图像中的数据

颜色空间

类型 gdalconst属性 整型值
未定义 GCI_Undefined 0
Greyscale GCI_GrayIndex 1
Paletted (see associated color table) GCI_PaletteIndex 2
Red band of RGBA image GCI_RedBand 3
Green band of RGBA image GCI_GreenBand 4
Blue band of RGBA image GCI_BlueBand 5
Alpha (0 = transparent; 255 = opaque) GCI_AlphaBand 6
Hue band of HLS image GCI_HueBand 7
Saturation band of HLS image GCI_SaturationBand 8
Lightness band of HLS image GCI_LightnessBand 9
Cyan band of CMYK image GCI_CyanBand 10
Magenta band of CMYK image GCI_MagentaBand 11
Yellow band of CMYK image GCI_YellowBand 12
Black band of CMLY image GCI_BlackBand 13
Y Luminance GCI_YCbCr_YBand 14
Cb Chroma GCI_YCbCr_CbBand 15
Cr Chroma GCI_YCbCr_CrBand 16

Table: gdalconst与整型的对应值

这里要注意，尽管在Python中GDAL的数据类型（表[tab:gdalconst]）与GDAL的color interpretation都是在gdalconst字典中， 但是在C语言中是在不同的枚举类型中定义的。

访问数据

我们通过ReadRaster读出的数据值只是对应到这个表的一个索引而已。 我们需要通过读出这些数据，并在真实数据表中找出真实数据， 重新组织成一个RGB表才能用来绘制。如果不经过对应， 绘制出来的东西可能没有任何意义。

使用自省函数help()：

In [6]:
from osgeo import gdal
dataset = gdal.Open("/gdata/lu75i1.tif")
band = dataset.GetRasterBand(1)
band.DataType
help(band.ReadAsArray)
help(band.ReadRaster)  

Help on method ReadAsArray in module osgeo.gdal:

ReadAsArray(xoff=0, yoff=0, win_xsize=None, win_ysize=None, buf_xsize=None, buf_ysize=None, buf_type=None, buf_obj=None, resample_alg=0, callback=None, callback_data=None) method of osgeo.gdal.Band instance
    Reading a chunk of a GDAL band into a numpy array. The optional (buf_xsize,buf_ysize,buf_type)
    parameters should generally not be specified if buf_obj is specified. The array is returned

Help on method ReadRaster in module osgeo.gdal:

ReadRaster(xoff=0, yoff=0, xsize=None, ysize=None, buf_xsize=None, buf_ysize=None, buf_type=None, buf_pixel_space=None, buf_line_space=None, resample_alg=0, callback=None, callback_data=None) method of osgeo.gdal.Band instance

可以看到这两个函数的原型：

ReadAsArray(self, xoff=0, yoff=0, win_xsize=None, win_ysize=None, buf_xsize=None, buf_ysize=None, buf_obj=None)

ReadRaster(self, int xoff, int yoff, int xsize, int ysize, int buf_xsize = None, 
int buf_ysize = None, int buf_type = None)

显然，band里的 ReadAsArray() 参数比 dataset 里面的要实用， 而ReadRaster则差不多。 但是ReadAsArray读出的是数组，可以用Numeric模块进行矩阵魔法。 ReadRaster读出的是二进制，虽然可以直接绘制， 但是对于一些绘图API来说， 对[[RRR...][GGG...][BBB...]]表的处理明显不如[[RGB][RGB]...]， 有的甚至不支持。 虽然可以用struct.unpack来拆封，可效率上就差很多(而且拆封出来还是数组)。 数组在矩阵魔法的控制之下则会显得十分方便快捷， 最后用 tostring 直接转化成为二进制绘制，速度也相当快。

好了，快速开始已经使我们可以初步看清楚了GDAL中图像的组织。 下面用一句话总结：波段组成图像，波段指挥颜色。

实际这个库主要是用来读取遥感数据的。 真正在图像处理方面还是不如PIL，两个其实是互用的。

读取索引图像中数据的问题

需要注意的是在gdal 使用ColorMap的时候， 对原始的ColorMap已经做过变动了。 比如geotiff的ColorMap的数据类型是 short， 默认的范围是在 0～65535， 但是在gdal读取出来以后，已经经过了范围压缩。 压缩范围是0～255。虽然都是short类型，但是值已经变化。

参考快速开始那张颜色表，可以看到颜色表中的数据是经过 （原始数据/65535.0*255）调整过的。这里在使用的时候可能比较方便， 但是如果你是有读取过geotiff原始数据背景的，则需要小心。 可能会有习惯性思维的陷阱。因为两者的类型都是short， 但是实际的数值不一样。

In [1]:
import helper; helper.info()

页面更新时间： 2019-04-14 22:05:20
操作系统/OS: Linux-4.9.0-8-amd64-x86_64-with-debian-9.8
Python: 3.5.3
版权归 OSGeo中国中心（http://www.osgeo.cn/ ） 所有
邮箱： bukun@iga.ac.cn