floating point grass rasters and mapserver - a resume
Neil Best
nbest at LANWORTH.COM
Wed Jul 25 10:25:40 EDT 2007
Also, is this work posted somewhere on the web as suggested? I did not
find it in the GRASS wiki, but maybe I have overlooked it. Thanks again.
Neil Best wrote:
> Can anyone comment on whether this information is still current? Are there
> any plans to accommodate GRASS rasters in Mapserver? Thanks for any info or
> tips.
>
> Neil
>
>
> On Sun, 16 Oct 2005 12:43:44 +0200, Jachym Cepicky
> <jachym.cepicky at CENTRUM.CZ> wrote:
>
>> 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