Tuesday, 12 January 2016

Showing only raster attribute table using GDAL?


I have a binary ArcGRID raster (AIG/Arc/Info Binary Grid) with an attribute table. I would like to get only this attribute table but gdalinfo outputs also projection information and bounds and then wanted GDALRasterAttributeTable XML part:


$ gdalinfo us_120fbfm13 -nogcp -nomd -noct filename
Driver: AIG/Arc/Info Binary Grid
...
Coordinate System is:
PROJCS["unnamed",
...
Center ( 1540680.000, 1553400.000) ( 78d42'42.88"W, 35d46' 8.50"N)
Band 1 Block=256x4 Type=Byte, ColorInterp=Undefined

Min=1.000 Max=99.000
NoData Value=255

...

How to show only raster attribut table using GDAL command line tools gdalinfo or perhaps gdal_translate?


Additional question which might be separate but might be also answered together with this is whether there is a possibility to get the attribute table in format such as CSV and not XML.


I'm on Linux, so I can actually do something like:


gdalinfo us_120fbfm13 -nogcp -nomd -noct -nofl | grep "<.*>"


but I would rather use some more general and possibly cross-platform solution.



Answer



Unless you use os/shell specific operations i.e. grep etc... there's no way that I know of to do this from the command line with gdalinfo.


A python alternative is to use the GetDefaultRAT() method and output to XML (which is what gdalinfo does):


import sys
from osgeo import gdal
import xml.etree.ElementTree as et #I much prefer the lxml library
#but it's 3rd party and must be installed

def toxml(rat):

#partial port of Serialize() from svn.osgeo.org/gdal/trunk/gdal/gcore/gdal_rat.cpp

rattree=et.Element('GDALRasterAttributeTable')

coltypes=[]

#Define each column.
icolcount=rat.GetColumnCount()
for icol in range(icolcount):
field=et.SubElement(rattree,'FieldDefn', {'index':str(icol)})

fname=et.SubElement(field,'Name')
ftype=et.SubElement(field,'Type')
fusage=et.SubElement(field,'Usage')

itype=rat.GetTypeOfCol(icol)

fname.text=rat.GetNameOfCol(icol)
ftype.text=str(rat.GetTypeOfCol(icol))
fusage.text=str(rat.GetUsageOfCol(icol))


if itype==gdal.GFT_Integer:
coltypes.append(['%s', rat.GetValueAsInt])
elif itype==gdal.GFT_Real:
coltypes.append(['%.16g', rat.GetValueAsDouble])
else:
coltypes.append(['%s', rat.GetValueAsString])

#Write out each row.
irowcount = rat.GetRowCount()
for irow in range(irowcount):

row=et.SubElement(rattree,'Row', {'index':str(irow)})
for icol in range(icolcount):
rowcol=et.SubElement(row,'F')
itype=rat.GetTypeOfCol(icol)
if itype==gdal.GFT_Integer:
rowcol.text='%s'%rat.GetValueAsInt(irow,icol)
elif itype==gdal.GFT_Real:
rowcol.text='%.16g'%rat.GetValueAsDouble(irow,icol)
else:
rowcol.text='%s'%rat.GetValueAsString(irow,icol)


return et.tostring(rattree)

if __name__ == '__main__':
ds=gdal.Open(sys.argv[1])
rat=ds.GetRasterBand(1).GetDefaultRAT()
sxml=serialize(rat)

print sxml


#Or if you want to pretty print
import xml.dom.minidom as md
print md.parseString(sxml).toprettyxml()

and here's a function to write out a CSV file.


import csv
import sys

from osgeo import gdal


def tocsv(rat,filepath):
with open(filepath, 'wb') as csvfile:
csvwriter = csv.writer(csvfile)

#Write out column headers
icolcount=rat.GetColumnCount()
cols=[]
for icol in range(icolcount):
cols.append(rat.GetNameOfCol(icol))
csvwriter.writerow(cols)


#Write out each row.
irowcount = rat.GetRowCount()
for irow in range(irowcount):
cols=[]
for icol in range(icolcount):
itype=rat.GetTypeOfCol(icol)
if itype==gdal.GFT_Integer:
value='%s'%rat.GetValueAsInt(irow,icol)
elif itype==gdal.GFT_Real:

value='%.16g'%rat.GetValueAsDouble(irow,icol)
else:
value='%s'%rat.GetValueAsString(irow,icol)
cols.append(value)
csvwriter.writerow(cols)

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...