floating point grass rasters and mapserver - a resume

Sun Oct 16 03:35:14 PDT 2005

Some time ago I was dealing with mapserver getting maps from a GRASS location. 
The maps are all floating point maps and as you (whoever interested) probably 
followed in a few mails, the result of most of my maps simply resulted in a 
single color layer (in my case it was yellow).
Frank Warmerdam helped me to solve what was solvable and workaround what was 
"workaroundable".
The following is an extract and resume of the mails exchange.
Any reply and correction is appreciated.

Alright, the problem was not only one, but a bunch concurring:
1) GRASS novalues are not supported and give problems

Running gdalinfo -mm on the map, we find that part of the problem is that 
nodata values are being treated as "nan" (not a number).  GDAL and mapserver 
do not know how to excluse "nan" nodata values from the min/max calculation 
so the autoscaling gets all screwed up.

Here we can see it in the gdalinfo output:
<SNIPPET>
Band 1 Block=385x1 Type=Float64, ColorInterp=Palette
  Min=0.000 Max=229172.410   Computed Min/Max=nan,nan
  NoData Value=nan
<SNIPPET>

Therefore the min and max have to ge defined and supplied to the map file with 
the following:
PROCESSING SCALE=min max

2) There is a problem with the creation of color ramps in the case of a big 
number of not homogeneous distributed values:

<SNIPPET>
  Color Table (RGB with 100001 entries)
    0: 0,0,0,0
    1: 255,255,0,255
    2: 255,255,0,255
    3: 255,255,0,255
    ...
<SNIPPET>

This is why I had a yellow map.  The first 256 values of the color table were 
used by MapServer for the after-scaling values.  Effectively the first 256 
entries of the color table are all shades of yellow.

3) MapServer does not support the GRASS color rules metadata

<SNIPPET>
Metadata:
  COLOR_TABLE_RULES_COUNT=5
  COLOR_TABLE_RULE_RGB_0=1.000000e+00 4.583448e+04 255 255 0 0 255 0
  COLOR_TABLE_RULE_RGB_1=4.583448e+04 9.166896e+04 0 255 0 0 255 255
  COLOR_TABLE_RULE_RGB_2=9.166896e+04 1.375034e+05 0 255 255 0 0 255
  COLOR_TABLE_RULE_RGB_3=1.375034e+05 1.833379e+05 0 0 255 255 0 255
  COLOR_TABLE_RULE_RGB_4=1.833379e+05 2.291724e+05 255 0 255 255 0 0
<SNIPPET/>

The metadata contains the coloring rules that should be applied. However, 
mapserver does not currently know how to use this sort of color rule 
metadata. For now the only solution is to recreate classes that are value 
similar. One problem though is that MapServer classes don't allow applying a 
range of color to a range of values like GRASS (or QGIS).  Applying own 
CLASSes in the LAYER makes mapserver ignore the color table.

The CLASSes option has two main drawbacks:

a) The CLASS does not allow color ranges. There is no big solution (only 
workarounds), if not to wait on the new range based class coloring that Bill 
Binko has been working on. An implementation of that is available in 4.6.x 
but it is under revision.

b) If we decide to create a set of CLASS rules to get the map drawn "as if" it 
was interpolated, you will need to create a lot of rules, which gives big big 
performance problems in drawing the map. This can however be solved with a 
PROCESSING command. Here the explanation:
The problem is that the class lookup is done once for each "bucket" in the 
lookup table and that by default there are 65536 buckets. It helps a lot to 
add the following line in the LAYER definition:

  PROCESSING "SCALE_BUCKETS=200"

This basically alters the resolution of the lookup table but substantially 
accelerates things.

After some testing on what would be the best to do for now, I'm now letting 
the following bash script snippet for automated creation of the raster part 
mapfile. It contains the solutions and workarounds reported in this document:

    # the path to the grass cellheader file
    absolutepath=$mapsetpath/cellhd/$i

    echo "LAYER" > rastertmpfile
    echo "  NAME $i" >> rastertmpfile
    echo "  TYPE RASTER" >> rastertmpfile
    echo "  STATUS ON" >> rastertmpfile
    echo "  DATA \"$absolutepath\"" >> rastertmpfile

    # extract the min and max
    minmax=`gdalinfo -mm $absolutepath | grep Min | awk -F "=" '{print $2, 
$3}' | awk '{print $1, $3}'`
    # extract the gdal computed min and max
    gdalrange=`gdalinfo -mm $absolutepath | grep "Computed Min" |awk -F 
"Computed Min/Max=" '{print $2}' | awk -F "," '{print $1}'`

    # if gdal has problems with colortables and novalues, it is not able to 
compute the range,
    # therefore we have to supply it and create a set of ad hoc colorrules
    if [ "$gdalrange" == "nan" ]
    then    
        # minmax has the format "123 234", so just extract them
        min=`echo $minmax | awk '{print $1}' | awk -F "." '{print $1}'`
        max=`echo $minmax | awk '{print $2}' | awk -F "." '{print $1}'`

        # define a default colortable, rainbow of grass, but with "bin + 1" 
number of rules,
        # to get more colors and give the sensation to be interpolated
        bins=19
        delta=$((($max-$min)/$bins))

        if [ $delta -ne 0 ]
        then
            colors[1]="255 255 0"

            colors[2]="195 255 0"
            colors[3]="130 255 0"
            colors[4]="65  255 0"

            colors[5]="0   255 0"

            colors[6]="0   255 65"
            colors[7]="0   255 130"
            colors[8]="0   255 195"

            colors[9]="0   255 255"

            colors[10]="0   195 255"
            colors[11]="0   130 255"
            colors[12]="0   65  255"

            colors[13]="0   0   255"

            colors[14]="65  0   255"
            colors[15]="130 0   255"
            colors[16]="195 0   255"

            colors[17]="255 0   255"
            colors[18]="255 0   195"
            colors[19]="255 0   130"
            colors[20]="255 0   65"

            # add the processing commands
            echo "  PROCESSING \"SCALE=$min $max\"" >> rastertmpfile
            echo "  PROCESSING \"SCALE_BUCKETS=100\"" >> rastertmpfile

            # create the CLASSes, labeling them with the pixel value
            first=$min
            second=$(($min+$delta))
            for j in `seq 1 $bins`
            do    
                echo "  CLASS " >> rastertmpfile
                echo "      NAME \"$first\"" >> rastertmpfile
                echo "      EXPRESSION ([pixel] >= $first and[pixel] < 
$second)" >> rastertmpfile
                echo "      COLOR ${colors[$j]}" >> rastertmpfile
                echo "  END" >> rastertmpfile
                first=$(($first+$delta))
                second=$(($second+$delta))
            done   

            # add a last CLASSes, due to bash maths limitations
            echo "  CLASS " >> rastertmpfile
            echo "      NAME \"$max\"" >> rastertmpfile
            echo "      EXPRESSION ([pixel] >= $first and[pixel] < $max)" >> 
rastertmpfile
            echo "      COLOR ${colors[20]}" >> rastertmpfile
            echo "  END" >> rastertmpfile
        fi
    fi    

Hope this helps someone, the script is not very stylistic, but at least it 
does the job for me. :)

Ciao,
Andrea

-- 
____________________________________________________________________________
HydroloGIS - Environmental Open Source Solutions
www.hydrologis.com

Andrea Antonello
Environmental Engineer
mobile:  +393288497722

"Let it be as much a great honour to take as to give learning,
if you want to be called wise."
Skuggsja' - The King's mirror - 1240 Reykjavik
____________________________________________________________________________