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