>>> 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表示Alpha

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