[mapserver-users] Continuous Raster Data

[Frank Warmerdam]

David W Bitner wrote:
>  I am looking at using mapserver to serve raster information for several
>  raster datasets.  These datasets are continuos floating points gridded at
>  1km resolution across the coterminous US.  I plan to use tiled geotiffs
> with
>  overlays for display using gdal support for the tiffs.  My questions are
> the
>  following:
> 
>  --- Is it possible to query the values in the raster?

David,

To the best of my knowledge there isn't any support in mapserver or mapscript
for querying values from raster layers.  If you use python mapscript you
could likely use the direct GDAL python bindings to access data values but
it would be going around MapServers abstraction of a raster layer which is
fairly simple.

>  ---  Each of the different rasters that I will be using deal with snow
> cover
>  information.... due to the nature of this, different scaling will be
>  necessary in different areas of the country when zoomed in (ie if one zooms
>  into an area including both the mountains in Southern California and the
>  desert just to the east, the value range will be huge, whereas it would be
>  important to be able to distinguish between subtle differences when zoomed
>  to the snowpack in the Red River Valley of the North).  Is it possible(or,
>  what would it take to make it possible) to be able to create a linear
>  stretch and histogram equalization for just the data values that are
>  displayed on the screen.

This isn't currently possible, but if you are adventurous you could hack the
drawGDAL() code for your local needs.  I would suggest passing extra
information in as metadata from the mapfile, and then keying on this in
drawGDAL().  Currently the greyscale code does this:


   /*
    * Read the entire region of interest in one gulp.  GDAL will take
    * care of downsampling and window logic.
    */
   pabyRaw1 = (unsigned char *) malloc(dst_xsize * dst_ysize);
   if( pabyRaw1 == NULL )
       return -1;

   GDALRasterIO( hBand1, GF_Read, src_xoff, src_yoff, src_xsize, src_ysize,
                 pabyRaw1, dst_xsize, dst_ysize, GDT_Byte, 0, 0 );


   /*
   ** Process single band using the provided, or greyscale colormap
   */
   if( hBand2 == NULL )
   {
       k = 0;
       for( i = dst_yoff; i < dst_yoff + dst_ysize; i++ )
       {
           int	result;

           for( j = dst_xoff; j < dst_xoff + dst_xsize; j++ )
           {
               result = cmap[pabyRaw1[k++]];
               if( result != -1 )
#ifndef USE_GD_1_2
                   img->pixels[i][j] = result;
#else
                   img->pixels[j][i] = result;
#endif
           }
       }
   }

Basically, the GDALRasterIO() call is loading into an 8bit buffer.  Any values
over 255 are just trimmed to 255 ... no scaling is applied.  You would need
to allocate a float buffer, and pass GDT_Float32 as the data type in the
RasterIO() call.  Then you could do more interesting sorts of transformations
on the data in the following loop, including histograming for equalization
and so forth.

I hope this gives you a few hints to get into this.

I would like to incorporate some kind of improved scaling and handling of
non-eight bit data in drawGDAL() eventually, but what you want to do is
fairly sophisticated and likely would have been fully addressed anyways.

Have fun, and drop me a line if you need more information.

-- 
---------------------------------------+--------------------------------------
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    | Geospatial Programmer for Rent