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

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

$ gdalinfo lu75i1.tif

......

251: 80,249,240,255

252: 114,239,249,255

253: 65,254,253,255

......

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

In [11]:
from osgeo import gdal
from osgeo import gdalconst
dataset = gdal.Open('/gdata/lu75i1.tif')
band = dataset.GetRasterBand(1)
band.GetRasterColorInterpretation()
Out[11]:
2
In [12]:
import gdalconst
gdalconst.GCI_PaletteIndex
Out[12]:
2
In [13]:
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 [14]:
from osgeo import gdal
dataset = gdal.Open("/gdata/lu75i1.tif")
band = dataset.GetRasterBand(1)
colormap = band.GetRasterColorTable()
colormap.GetPaletteInterpretation()
Out[14]:
1
In [15]:
colormap.GetCount()
Out[15]:
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 [16]:
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, 但是实际的数值不一样。