图像 lu75i1.tif 是索引图像,使用gdalinfo来查看的话,可以看见有如下的信息:
!gdalinfo /gdata/lu75i1.tif
ERROR 4: /gdata/lu75i1.tif: No such file or directory gdalinfo failed - unable to open '/gdata/lu75i1.tif'.
显示的是索引图像的颜色查找表。
from osgeo import gdal
from osgeo import gdalconst
dataset = gdal.Open('/data/gdata/lu75i1.tif')
band = dataset.GetRasterBand(1)
band.GetRasterColorInterpretation()
/opt/conda/lib/python3.12/site-packages/osgeo/gdal.py:311: FutureWarning: Neither gdal.UseExceptions() nor gdal.DontUseExceptions() has been explicitly called. In GDAL 4.0, exceptions will be enabled by default. warnings.warn(
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, 我们知道我们的图像是一个颜色表索引的图像, 而不是纯粹的黑白灰度图像。这意味着我们读出的数据有可能不是真实的数据。 这些数据只是一个个实际数据的索引。真实数据存储在另一个表中。
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("/data/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个值都有意义了)。
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 |
from osgeo import gdal
dataset = gdal.Open("/data/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
Read a window of this raster band into a NumPy array.
Parameters
----------
xoff : float, default=0
The pixel offset to left side of the region of the band to
be read. This would be zero to start from the left side.
yoff : float, default=0
The line offset to top side of the region of the band to
be read. This would be zero to start from the top side.
win_xsize : float, optional
The number of pixels to read in the x direction. By default,
equal to the number of columns in the raster.
win_ysize : float, optional
The number of rows to read in the y direction. By default,
equal to the number of bands in the raster.
buf_xsize : int, optional
The number of columns in the returned array. If not equal
to ``win_xsize``, the returned values will be determined
by ``resample_alg``.
buf_ysize : int, optional
The number of rows in the returned array. If not equal
to ``win_ysize``, the returned values will be determined
by ``resample_alg``.
buf_type : int, optional
The data type of the returned array
buf_obj : np.ndarray, optional
Optional buffer into which values will be read. If ``buf_obj``
is specified, then ``buf_xsize``/``buf_ysize``/``buf_type``
should generally not be specified.
resample_alg : int, default = :py:const:`gdal.GRIORA_NearestNeighbour`.
Specifies the resampling algorithm to use when the size of
the read window and the buffer are not equal.
callback : function, optional
A progress callback function
callback_data: optional
Optional data to be passed to callback function
Returns
-------
np.ndarray
Examples
--------
>>> import numpy as np
>>> ds = gdal.GetDriverByName("GTiff").Create("test.tif", 4, 4, eType=gdal.GDT_Float32)
>>> ds.WriteArray(np.arange(16).reshape(4, 4))
0
>>> band = ds.GetRasterBand(1)
>>> # Reading an entire band
>>> band.ReadAsArray()
array([[ 0., 1., 2., 3.],
[ 4., 5., 6., 7.],
[ 8., 9., 10., 11.],
[12., 13., 14., 15.]], dtype=float32)
>>> # Reading a window of a band
>>> band.ReadAsArray(xoff=2, yoff=2, win_xsize=2, win_ysize=2)
array([[10., 11.],
[14., 15.]], dtype=float32)
>>> # Reading a band into a new buffer at higher resolution
>>> band.ReadAsArray(xoff=0.5, yoff=0.5, win_xsize=2.5, win_ysize=2.5, buf_xsize=5, buf_ysize=5)
array([[ 0., 1., 1., 2., 2.],
[ 4., 5., 5., 6., 6.],
[ 4., 5., 5., 6., 6.],
[ 8., 9., 9., 10., 10.],
[ 8., 9., 9., 10., 10.]], dtype=float32)
>>> # Reading a band into an existing buffer at lower resolution
>>> band.ReadAsArray(buf_xsize=2, buf_ysize=2, buf_type=gdal.GDT_Float64, resample_alg=gdal.GRIORA_Average)
array([[ 2.5, 4.5],
[10.5, 12.5]])
>>> buf = np.zeros((2,2))
>>> band.ReadAsArray(buf_obj=buf)
array([[ 5., 7.],
[13., 15.]])
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,两个其实是互用的。
读取索引图像中数据的问题
需要注意的是在gdal 使用ColorMap的时候, 对原始的ColorMap已经做过变动了。
比如geotiff的ColorMap的数据类型是 short, 默认的范围是在 0~65535,
但是在gdal读取出来以后,已经经过了范围压缩。
压缩范围是0~255。虽然都是short类型,但是值已经变化。
参考快速开始那张颜色表,可以看到颜色表中的数据是经过 (原始数据/65535.0*255)调整过的。这里在使用的时候可能比较方便, 但是如果你是有读取过geotiff原始数据背景的,则需要小心。 可能会有习惯性思维的陷阱。因为两者的类型都是short, 但是实际的数值不一样。