[Gdal-dev] gdal_translate and Dynamic Range

[Gregory, Matthew]

Frank Warmerdam wrote:
> Scott Elko wrote:
> > Can anybody share any ideas in how to get gdla_translate to produce 
> > the same result as that of OpenEV? Or is it an impossible feat that 
> > I'm trying to accomplish?
> 
> gdal_translate isn't going to have the same logic as OpenEV.  
> So you either have to compute the scaling min/max by some 
> other means (perhaps inspecting it in OpenEV or a script) and 
> use those with gdal_translate, or you are going to have to 
> emulate the logic of OpenEV in a script of some sort.  In 
> fact, depending on the circumstances you might be able to 
> write a Python script that imports the OpenEV core, and use 
> it to do the operation, but does not actually create any user 
> interface components.  This could be tricky though.

Scott,

I think I was addressing a similar desire to get a stretched image on a
subwindow in my recent post
(http://lists.maptools.org/pipermail/gdal-dev/2007-February/011872.html)
.  Frank made the excellent suggestion to write a Python script to
calculate the min/max on that subwindow in code and then stretch using
that remapping.  Here are the steps that my Python script performs:

1. Read in the subwindow's coordinates and calculate its offsets and
size from the source image.

2. Create a VRT XML string which stores the specification of your
subwindow.

3. Open this string with gdal.Open() and calculate the "local" min/max,
e.g. 

  vrt_ds = gdal.Open(vrtXmlStr)
  vrtBand = vrt_ds.GetRasterBand(band)
  (zMinLocal, zMaxLocal) = vrtBand.ComputeRasterMinMax()

4. Create the subsetted image using this local min/max.  Using the
raster_copy function from gdal_merge.py as a template and importing
numpy, I modified the code to this:

  srcBand = srcFile.GetRasterBand( srcBandNum )
  dstBand = dstFile.GetRasterBand( dstBandNum )

  max = zMaxLocal
  min = zMinLocal
  nodata = noDataValue
  x = srcBand.ReadAsArray(src_xoff, src_yoff, src_xsize, src_ysize,
                            dst_xsize, dst_ysize, dstBand.DataType)
  m = (255.0 - 1.0) / (max - min)
  b = 255.0 - m * max
  y = numpy.select([x == nodata], [0], default = numpy.round(m * x + b))
  z = y.astype(numpy.uint8)
  dstBand.WriteArray(z[0])

This creates an unsigned 8-bit image with nodata mapped to 0 and a
linear scaling of your source local min/max to the range 1-255.  Of
course, with no error checking! 

I confess that I haven't looked at what OpenEV gives you, but I think
this gets at the functionality for which you are looking.

Hope this helps,

Matt Gregory
Faculty Research Assistant
Department of Forest Science
Oregon State University
Phone : (541) 758-7778
Email : matt.gregory at oregonstate.edu