[Gdal-dev] GDAL vs ArcGIS raster calculator

[Frank Warmerdam]

joosthoek wrote:
> Hi,
> 
> 
> I was wondering if someone knows if there is a GDAL alternative for the
> following ArcGIS Raster Calculator code:
> 
> select ( merge ( setnull ( [raster] > 32767 , [raster] ) , ( [raster] -
> 65536 ) ) , 'value <> -32768' )
> 
> The input raster is 16-bit DTM data and as far as I know this code sets the
> values to meters with correct datum.
> 
> I was trying the gdal_translate -scale option but that didn't work somehow
> (I tried -scale 0 65535 -32768 32767). Someone has a better idea?

Jelmer,

The calculation appears to be correcting for 16 bit signed integer data
that was incorrectly treated as unsigned.  It also marks post-adjusted
values of -32768 as nodata.

You should be able to do something similar with a numpy python script like
the following (which works with FWTools):

#!/usr/bin/env python

import gdal
import gdalnumeric

x = gdalnumeric.LoadFile('utm.tif')

x_adjusted = gdalnumeric.choose( \
     gdalnumeric.greater( x, 32767 ),  (x, x-65536) )

x_adjusted = x_adjusted.astype(gdalnumeric.Int16)

prototype_ds = gdal.Open( 'utm.tif' )
gdalnumeric.SaveArray( x_adjusted, 'utm2.tif', prototype = prototype_ds )


Note, this script does not set the nodata value.  That could be done as
an additional step with gdal_translate, or is left as an exercise for the
reader. :-)

Best regards,
-- 
---------------------------------------+--------------------------------------
I set the clouds in motion - turn up   | Frank Warmerdam, warmerdam at pobox.com
light and sound - activate the windows | http://pobox.com/~warmerdam
and watch the world go round - Rush    | President OSGeo, http://osgeo.org