2

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.

Improve this question
4
  • 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, 2015 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, 2015 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, 2015 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, 2015 at 23:26

1 Answer 1

Reset to default
3

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)
Improve this answer

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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