为栅格数据添加属性表

学习《Python地理数据处理》的笔记

import os
import gdal

os.chdir(r'E:\pycode\DATA\osgeopy-data\Switzerland')
ds = gdal.Open('dem_class2.tif')
band = ds.GetRasterBand(1)
band.SetNoDataValue(-1)

rat = gdal.RasterAttributeTable()
rat.CreateColumn('Value', gdal.GFT_Integer, gdal.GFU_Name)
rat.CreateColumn('Count', gdal.GFT_Integer, gdal.GFU_PixelCount)
rat.CreateColumn('Elevation', gdal.GFT_String, gdal.GFU_Generic)
rat.SetRowCount(6)

rat.WriteArray(range(6), 0)
rat.WriteArray(band.GetHistogram(-0.5, 5.5, 6, False, False), 1)
rat.SetValueAsString(1, 2, '0 - 800')
rat.SetValueAsString(2, 2, '800 - 1300')
rat.SetValueAsString(3, 2, '1300 - 2000')
rat.SetValueAsString(4, 2, '2000 - 2600')
rat.SetValueAsString(5, 2, '2600 +')

band.SetDefaultRAT(rat)
band.SetNoDataValue(0)
del band, ds

生成属性表后，在同名文件.tif.aux中会以XML形式记录（部分）：

<GDALRasterAttributeTable>
      <FieldDefn index="0">
        <Name>Value</Name>
        <Type>0</Type>
        <Usage>2</Usage>
      </FieldDefn>
      <FieldDefn index="1">
        <Name>Count</Name>
        <Type>0</Type>
        <Usage>1</Usage>
      </FieldDefn>
      <FieldDefn index="2">
        <Name>Elevation</Name>
        <Type>2</Type>
        <Usage>0</Usage>
      </FieldDefn>
      <Row index="0">
        <F>0</F>
        <F>55467</F>
        <F>NoData</F>
      </Row>
      <Row index="1">
        <F>1</F>
        <F>27213</F>
        <F>0 - 800</F>
      </Row>
      <Row index="2">
        <F>2</F>
        <F>12986</F>
        <F>800 - 1300</F>
      </Row>
      <Row index="3">
        <F>3</F>
        <F>13642</F>
        <F>1300 - 2000</F>
      </Row>
      <Row index="4">
        <F>4</F>
        <F>10632</F>
        <F>2000 - 2600</F>
      </Row>
      <Row index="5">
        <F>5</F>
        <F>5414</F>
        <F>2600 +</F>
      </Row>
    </GDALRasterAttributeTable>