>>> from env_helper import info; info()
页面更新时间: 2023-06-26 15:11:00
运行环境:
Linux发行版本: Debian GNU/Linux 12 (bookworm)
操作系统内核: Linux-6.1.0-9-amd64-x86_64-with-glibc2.36
Python版本: 3.11.2
2.5. 访问索引图像的处理¶
图像 lu75i1.tif
是索引图像,使用gdalinfo来查看的话,可以看见有如下的信息:
>>> !gdalinfo /gdata/lu75i1.tif
Driver: GTiff/GeoTIFF
Files: /gdata/lu75i1.tif
/gdata/lu75i1.tif.ovr
/gdata/lu75i1.tif.aux.xml
Size is 6122, 4669
Coordinate System is:
PROJCRS["Albers_Conic_Equal_Area",
BASEGEOGCRS["GCS_X",
DATUM["unknown",
ELLIPSOID["unretrievable_using_WGS84",6378137,298.257223563,
LENGTHUNIT["metre",1,
ID["EPSG",9001]]]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433,
ID["EPSG",9122]]]],
CONVERSION["Albers Equal Area",
METHOD["Albers Equal Area",
ID["EPSG",9822]],
PARAMETER["Latitude of false origin",0,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8821]],
PARAMETER["Longitude of false origin",105,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8822]],
PARAMETER["Latitude of 1st standard parallel",25,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8823]],
PARAMETER["Latitude of 2nd standard parallel",47,
ANGLEUNIT["degree",0.0174532925199433],
ID["EPSG",8824]],
PARAMETER["Easting at false origin",0,
LENGTHUNIT["metre",1],
ID["EPSG",8826]],
PARAMETER["Northing at false origin",0,
LENGTHUNIT["metre",1],
ID["EPSG",8827]]],
CS[Cartesian,2],
AXIS["easting",east,
ORDER[1],
LENGTHUNIT["metre",1,
ID["EPSG",9001]]],
AXIS["northing",north,
ORDER[2],
LENGTHUNIT["metre",1,
ID["EPSG",9001]]]]
Data axis to CRS axis mapping: 1,2
Origin = (1852951.760316815227270,5309350.360150607302785)
Pixel Size = (30.000000000000000,-30.000000000000000)
Metadata:
AREA_OR_POINT=Area
DataType=Thematic
Image Structure Metadata:
COMPRESSION=LZW
INTERLEAVE=BAND
Corner Coordinates:
Upper Left ( 1852951.760, 5309350.360) (129d35'30.45"E, 46d56'20.01"N)
Lower Left ( 1852951.760, 5169280.360) (129d 8'55.08"E, 45d43'10.61"N)
Upper Right ( 2036611.760, 5309350.360) (131d54'59.43"E, 46d30'55.88"N)
Lower Right ( 2036611.760, 5169280.360) (131d26' 7.21"E, 45d18'19.30"N)
Center ( 1944781.760, 5239315.360) (130d31'28.35"E, 46d 7'26.29"N)
Band 1 Block=128x128 Type=Byte, ColorInterp=Palette
Min=0.000 Max=18.000
Minimum=0.000, Maximum=18.000, Mean=3.791, StdDev=3.385
NoData Value=255
Overviews: 3061x2335, 1531x1168, 766x584, 383x292, 192x146
Metadata:
RepresentationType=THEMATIC
STATISTICS_COVARIANCES=11.45841938446506
STATISTICS_MAXIMUM=18
STATISTICS_MEAN=3.7912636531876
STATISTICS_MINIMUM=0
STATISTICS_SKIPFACTORX=1
STATISTICS_SKIPFACTORY=1
STATISTICS_STDDEV=3.3850287125023
Color Table (RGB with 256 entries)
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
显示的是索引图像的颜色查找表。
>>> from osgeo import gdal
>>> from osgeo import gdalconst
>>> dataset = gdal.Open('/gdata/lu75i1.tif')
>>> band = dataset.GetRasterBand(1)
>>> band.GetRasterColorInterpretation()
2
>>> from osgeo import gdalconst
>>> gdalconst.GCI_PaletteIndex
2
>>> colormap = band.GetRasterColorTable()
>>> dir(colormap)
>>> colormap.GetCount()
>>> colormap.GetPaletteInterpretation()
>>> gdalconst.GPI_RGB
>>> for i in range(colormap.GetCount() - 10, colormap.GetCount()):
>>> print("%i:%s" % (i, colormap.GetColorEntry(i)))
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, 我们知道我们的图像是一个颜色表索引的图像, 而不是纯粹的黑白灰度图像。这意味着我们读出的数据有可能不是真实的数据。 这些数据只是一个个实际数据的索引。真实数据存储在另一个表中。
2.5.1. 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表示AlphaGPI_CYMP
: c1表示青,c2表示品,c3表示黄,c4表示黑GPI_HLS
: c1表示色调,c2表示亮度,c3表示饱和度
虽然有颜色表示数的区别,但是用GetColorEntry读出的都是4个值, 根据PaletteInterp枚举值看截取其中的几个值形成的颜色。
>>> from osgeo import gdal
>>> dataset = gdal.Open("/gdata/lu75i1.tif")
>>> band = dataset.GetRasterBand(1)
>>> colormap = band.GetRasterColorTable()
>>> colormap.GetPaletteInterpretation()
1
>>> colormap.GetCount()
256
通过用GetRasterColorTable获得了颜色表, 通过GetPaletteInterpretation知道我们获得的颜色表是一个RGB颜色表。 GDAL支持多种颜色表,具体可以参考gdalconst模块中GPI打头的枚举值。 然后我们可以通过GetCount来获得颜色数量, 通过GetColorEntry获得颜色表中的值。这里的颜色值都是一个4值的元组。 里面有意义的只有前三个(如果颜色模型是GPI_RGB, GPI_HLS, 都使用前3个,如果采用GPI_CMYK,则4个值都有意义了)。
2.5.2. 颜色空间¶
gdalconst
与整型的对应值
类型 | 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 |
这里要注意,尽管在Python中GDAL的数据类型(表[tab:gdalconst])与GDAL的color interpretation都是在gdalconst字典中, 但是在C语言中是在不同的枚举类型中定义的。
2.5.3. 访问数据¶
我们通过ReadRaster读出的数据值只是对应到这个表的一个索引而已。 我们需要通过读出这些数据,并在真实数据表中找出真实数据, 重新组织成一个RGB表才能用来绘制。如果不经过对应, 绘制出来的东西可能没有任何意义。
使用自省函数help()
:
>>> from osgeo import gdal
>>> dataset = gdal.Open("/gdata/lu75i1.tif")
>>> band = dataset.GetRasterBand(1)
>>> band.DataType
>>> help(band.ReadAsArray)
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(band.ReadRaster)
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, buf_obj=None) method of osgeo.gdal.Band instance
可以看到这两个函数的原型:
显然,band里的 ReadAsArray()
参数比 dataset
里面的要实用,
而ReadRaster则差不多。
但是ReadAsArray读出的是数组,可以用Numeric模块进行矩阵魔法。
ReadRaster读出的是二进制,虽然可以直接绘制, 但是对于一些绘图API来说,
对[[RRR…][GGG…][BBB…]]表的处理明显不如[[RGB][RGB]…], 有的甚至不支持。
虽然可以用struct.unpack来拆封,可效率上就差很多(而且拆封出来还是数组)。
数组在矩阵魔法的控制之下则会显得十分方便快捷, 最后用 tostring
直接转化成为二进制绘制,速度也相当快。
好了,快速开始已经使我们可以初步看清楚了GDAL中图像的组织。 下面用一句话总结:波段组成图像,波段指挥颜色。
实际这个库主要是用来读取遥感数据的。 真正在图像处理方面还是不如PIL,两个其实是互用的。
2.5.4. 读取索引图像中数据的问题¶
需要注意的是在gdal 使用ColorMap的时候, 对原始的ColorMap已经做过变动了。
比如geotiff的ColorMap的数据类型是 short
, 默认的范围是在
0~65535, 但是在gdal读取出来以后,已经经过了范围压缩。
压缩范围是0~255。虽然都是short类型,但是值已经变化。
参考快速开始那张颜色表,可以看到颜色表中的数据是经过 (原始数据/65535.0*255)调整过的。这里在使用的时候可能比较方便, 但是如果你是有读取过geotiff原始数据背景的,则需要小心。 可能会有习惯性思维的陷阱。因为两者的类型都是short, 但是实际的数值不一样。