图像 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, 但是实际的数值不一样。