[gdal-dev] Creating aux.xml files for GTiff using Python

Jason Roberts jason.roberts at duke.edu
Thu Dec 9 10:21:23 EST 2010


Matt,

Since nobody else more knowledgeable replied and I was in the original
discussion, I'll take a shot.

To review: for some situations and formats, the following code would work:

statistics = band.GetStatistics(False, True)     # Force calculation of
statistics
band.SetStatistics(*statistics)                  # Write statistics to the
file

But that was not working for GTiff.

And now a short digression. I had a further problem which I did not mention
on gdal-dev, best described by this comment in my code:

        # Do not use GDAL to calculate statistics or build pyramids.
        # Even though GDAL can do it faster than the geoprocessor, we
        # prefer to use the geoprocessor to maximize the compatibility
        # of the output raster with ArcGIS. For example, I do not know
        # how to get GDAL to build an ArcGIS-compatible histogram.
        # Without the histogram ArcGIS defaults to displaying integer
        # rasters using a black-to-white stretch rather than a unique
        # value classifier with random colors. ArcGIS users will
        # expect colors; to get the necessary histogram, we have to
        # calculate statistics with the geoprocessor.

By "histogram", I meant that when I loaded the integer raster into ArcMap
and tried to set the symbology to unique values, ArcMap prompted with
"Histogram does not exist. Create one?" Try as I might, I could not figure
out how to create the raster with GDAL in a way that did not result in that
prompt. I discovered that calling ArcGIS's own Calculate Statistics
geoprocessing tool solved the problem. Because that was ok in my situation,
I gave up on GDAL and just did that.

As a result, I never ended up implementing Frank's suggestion in my own
situation. But it looked unusual to me. In English, that call to
gdal_translate means "copy xxx.tif to xxx.aux, using HFA as the output
format (ERDAS .img format), and for that output, calculate statistics and
use an aux file". I'm not fully certain, but I think xxx.aux would then be a
fully copy of xxx.tif, but in HFA format. That seemed at first like a
bizarre suggestion, because I did not want a duplicate raster. But I
concluded that Frank was saying that a side effect of this procedure was to
create a .aux file that was compatible with the original GTiff. Presumably,
I would then delete the HFA file that I did not need.

If my understanding is correct, it should be possible to reproduce it with
the Python bindings by just calling CreateCopy with the equivalent
parameters. Have you tried that?

See also this thread
http://lists.osgeo.org/pipermail/gdal-dev/2010-November/026790.html. I'm not
sure whether or not that is relevant.

Hope that helps,

Jason


-----Original Message-----
From: gdal-dev-bounces at lists.osgeo.org
[mailto:gdal-dev-bounces at lists.osgeo.org] On Behalf Of Gregory, Matthew
Sent: Tuesday, December 07, 2010 6:22 PM
To: gdal-dev at lists.osgeo.org
Subject: [gdal-dev] Creating aux.xml files for GTiff using Python

Hi all,

There have been a few threads on this listserv about ArcGIS not recognizing
statistics built into GeoTiffs.  For example, see:

  http://lists.osgeo.org/pipermail/gdal-dev/2009-September/022044.html

Because GeoTiffs don't build .aux.xml files by default (even with
GDAL_PAM_ENABLED set to YES or at least I think this is the case), I'm
running into the same issues.  When a GeoTiff is brought into ArcMap, it
complains about statistics not being built and then builds them for you.  

Frank's fix, however, in the above thread:

  'gdal_translate -of HFA -co AUX=YES -co STATISTICS=YES xxx.tif xxx.aux'

looks like it first creates the accompanying .aux.xml file for the GeoTiff
(e.g. xxx.tif.aux.xml) and then creates two other files xxx.aux and
xxx.aux.aux.xml.  When I bring xxx.tif into ArcMap, it correctly reads the
statistics and everything is displayed correctly.

My question is whether or not it's possible to force creation of the
.aux.xml files through the Python bindings.  I've tried both
ComputeBandStats() and ComputeStatistics(False) with no success.

This falls under the category of annoyances, so not a huge issue.

Thanks for pointers.

matt


_______________________________________________
gdal-dev mailing list
gdal-dev at lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev



More information about the gdal-dev mailing list