1

I have projected a raster with gdal as part of an R workflow:

gdalwarp -s_srs EPSG:4283 -t_srs EPSG:3577 -r bilinear -multi -tr 1000 1000 -dstnodata -9999 in.tif out.tif

When I add out.tif to ArcMap, the minimum value indicated beside the colour ramp is incorrect (268.048, while the true minimum is 91.2931), and the colour ramp applied to the data is stretched between this incorrect minimum and the (correct) maximum.

I can calculate statistics (data management > raster > raster properties toolbox) to remedy this, but I'm concerned that if I fail to do this for one of the many rasters I'm working with, I (or my colleagues) might misinterpret a map.

The trick of switching to Classified symbology and back to Stretched returns the error: This dataset does not have a valid histogram required for classification method.... Attempting the same with Unique Values returns the error: The number of unique values exceeds the default limitation (>65536)....

Why has ArcMap chosen this seemingly arbitrary minimum of 268.048? Is there anything I can do to prevent it from happening (e.g. to have ArcMap either calculate the correct minimum when the raster is added, or to use default statistics that span the possible range of values given the pixel type and depth)?

NB: I am using a fresh installation of ArcGIS 10.2, and "Calculate statistics" is checked in Geoprocessing > Environments > Raster Storage.

  • 1
    Has R written an aux.xml file? Normally, when using GDAL tools like warp or translate the statistics don't exist and need to be built, which ArcMap can do (and does) in order to use a stretched renderer if it can't find pre-calculated values. If R is tricking ArcMap into using a bad value then you should calculate or remove the existing statistics. Try using ERDAS Imagine format (.img), that format is more likely to be 'bare' as it doesn't support tags and calculated properly. If you like I can provide code to calculate statistics for a folder of rasters.. – Michael Stimson Jan 5 '15 at 22:45
  • @Michael I'm calling gdalwarp from R with system, so the result is identical to if I run gdalwarp from the command line. No aux.xml file is created. Using -of HFA (ERDAS Imagine) produces aux.xml but ArcMap reads the same incorrect minimum. Same story with -of EHdr (.bil). Your offer of code is enticing - is it Python, or R-based? – jbaums Jan 5 '15 at 23:00
  • @Michael Ah - I see that gdalinfo -stats generates the required statistics and the aux.xml file. Once I do this, ArcMap displays the file as expected. A very easy fix - thank you for (maybe inadvertently!) pointing me in this direction. :) – jbaums Jan 5 '15 at 23:09
  • 1
    Code is arcpy, ready to add to a toolbox. It's a really simple task until you need to do it several thousand times. There is the batch option in the toolbox but I find it easier to tell it to find all the rasters rather than trying to supply multiple rasters. If you can do it with gdalinfo then I'd do it that way; note GdalAddO is recognized by Esri so you can do your pyramids that way. – Michael Stimson Jan 5 '15 at 23:26
2

Adding a call to gdalinfo -stats generates an aux.xml file for the raster and solves this problem.

My R code is then:

system('gdalwarp -s_srs EPSG:4283 -t_srs EPSG:3577 -r bilinear -multi -tr 1000 1000 -dstnodata -9999 in.tif out.tif')
system('gdalinfo -stats out.tif', show.output.on.console=FALSE)

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.