Saturday 22 October 2016

Python, GDAL and building raster attribute tables


I have an integer raster for which I would like to build a raster attribute table using Python and GDAL. I can create a GDAL raster attribute table in Python as follows:


>>> rat = gdal.RasterAttributeTable()

This works fine, as we can see:


>>> rat
>

The table thus created has no rows or columns:


>>> rat.GetRowCount()

0
>>> rat.GetColumnCount()
0

I create a column called "Value" to store the unique values in the raster:


>>> rat.CreateColumn("Value", gdalconst.GFT_Integer, gdalconst.GFU_MinMax)
0

This is fine, and the column count is updated:


>>> rat.GetColumnCount()

1

Now I have to add values (records) to the column for it to be of any use. I can get a list of unique values from the raster's band like so:


>>> data = band.ReadAsArray(0, 0, dataset.RasterXSize, dataset.RasterYSize)
>>> vals = list(numpy.unique(data))
>>> vals
[3, 7, 8, 10, 11, 12, 13, 14, 17, 18, 20, 22, 23, 25, 27, 28, 41, 45, 52, 56]

What I would like to do is create a for loop to loop through vals and populate the column in the attribute table. I thought I could do something like this:


>>> for i in range(len(vals)):

rat.SetValueAsInt(i, 0, vals[i])

...where i is the row (record), 0 is the field index and vals[i] is the integer value I want to insert. But it causes an error:


Traceback (most recent call last):
File "", line 2, in
rat.SetValueAsInt(i, 0, vals[i])
File "C:\Python27\lib\site-packages\osgeo\gdal.py", line 1139, in SetValueAsInt
return _gdal.RasterAttributeTable_SetValueAsInt(self, *args)
TypeError: in method 'RasterAttributeTable_SetValueAsInt', argument 4 of type 'int'


The error is caused because I use vals[i] in the call to SetValueAsInt() rather than using an integer directly. For example, rat.SetValueAsInt(0, 0, 0) works fine, but is useless for populating the column if I just want to loop over the list of unique values.


Is this a known issue? Google has so far not been very useful. What can I do to get around this problem?



Answer



The SetValueAsInt method is expecting a python int type, not a numpy uint16 type.


>>> print type(vals[0])


The following works:


rat.SetValueAsInt(i, 0, int(vals[i]))

No comments:

Post a Comment

arcpy - Changing output name when exporting data driven pages to JPG?

Is there a way to save the output JPG, changing the output file name to the page name, instead of page number? I mean changing the script fo...