floating point grass rasters and mapserver - a resume
Jachym Cepicky
jachym.cepicky at CENTRUM.CZ
Sun Oct 16 03:43:44 PDT 2005
very well!
I did it by hand till now (well I do not have so many raster in GRASS
and ... how often do one need to change the color table?)
could you post this script on GRASS-Addons site?
(http://grass.gdf-hannover.de/twiki/bin/view/GRASS/GrassAddOns)
Jáchym
On Sun, Oct 16, 2005 at 12:35:14PM +0200, Andrea Antonello wrote:
> 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
> ____________________________________________________________________________
--
Jachym Cepicky
e-mail: jachym.cepicky at centrum.cz
URL: http://les-ejk.cz
GPG: http://les-ejk.cz/gnupg_public_key/
-----------------------------------------
OFFICE:
Department of Geoinformation Technologies
LDF MZLU v Brně
Zemědělská 3
613 00 Brno
e-mail: xcepicky at node.mendelu.cz
URL: http://mapserver.mendelu.cz
Tel.: +420 545 134 514
More information about the MapServer-users
mailing list