Displaying GRIB data with MapServer

I recently had the opportunity to prototype WMS visualization of meteorological data.  MapServer, GDAL and Python to the rescue!  Here are the steps I took to make it happen.

The data, (GRIB), is a GDAL supported format, so MapServer can handle processing as a result.  The goal here was to create a LAYER object.  First thing was to figure out the projection, then figure out the band pixel values/ranges and correlate to MapServer classes (in this case I just used a simple greyscale approach).

Here’s the hack:

import sys
import osgeo.gdal as gdal
import osgeo.osr as osr

if len(sys.argv) < 3:
 print 'Usage: %s <file> <numclasses>' % sys.argv[0]
 sys.exit(1)

cvr = 256  # range of RGB values
numclasses = int(sys.argv[2])  # number of classifiers

ds = gdal.Open(sys.argv[1])

# get proj4 def and write out PROJECTION object
p = osr.SpatialReference()
s = p.ImportFromWkt(ds.GetProjection())
p2 = p.ExportToProj4().split()

print '  PROJECTION'
for i in p2:
 print '   "%s"' % i.replace('+','')
print '  END'

# get band pixel data ranges and classify
band = ds.GetRasterBand(1)
min = band.GetMinimum()
max = band.GetMaximum()

if min is None or max is None:  # compute automagically
 (min, max) = band.ComputeRasterMinMax(1)

# calculate range of pixel values
pixel_value_range = float(max - min)
# calculate the intervals of values based on classes specified
pixel_interval = pixel_value_range / numclasses
# calculate the intervals of color values
color_interval = (pixel_interval * cvr) / pixel_value_range

for i in range(numclasses):
 print '''  CLASS
  NAME "%.2f to %.2f"
  EXPRESSION ([pixel] >= %.2f AND [pixel] < %.2f)
  STYLE
   COLOR %s %s %s
  END
 END''' % (min, min+pixel_interval, min, min+pixel_interval, cvr, cvr, cvr)
 min += pixel_interval
 cvr -= int(color_interval)

Running this script outputs various bits for MapServer mapfile configuration.  Passing more classes to the script creates more CLASS objects, resulting in a smoother looking image.

Here’s an example GetMap request:

Meteorological data in GRIB format via MapServer WMS

Users can query to obtain pixel values (water temperature in this case) via GetFeatureInfo.  Given that these are produced frequently, we can use the WMS GetMap TIME parameter to create time series maps of the models.

4 Comments so far »

  1. David said,

    Wrote on June 18, 2010 @ 08:06:13

    Google Chrome 5.0.375.55 Windows XP

    Hi!

    Thanks for the guide.

    I have to do a similar test.
    Will be possible to have a .zip file containing the mapfile and file layers?

    Thank you very much!

    david

    Posted from United States
    Google Chrome 5.0.375.55 Windows XP
  2. tomkralidis said,

    Wrote on June 18, 2010 @ 13:00:20

    Mozilla Firefox 3.6.3 Mac OS X 10

    You can find some sample data at http://dd.weatheroffice.ec.gc.ca/model_gem_regional/high_resolution/grib1/00/.

    As for a mapfile def, simply follow the MapServer Raster Data Guide in the docs.

    To generate the classes, just run the script above.

    Posted from Canada Canada
    Mozilla Firefox 3.6.3 Mac OS X 10
  3. David said,

    Wrote on June 30, 2010 @ 03:05:01

    Google Chrome 5.0.375.55 Windows XP

    Hi!

    Now, I’m able to display correctly my grib file.
    Processing the band number 4 and using the getfeatureinfo (OpenLayer) I get the pixel value.

    It’s so difficult to find documentation about it.
    Now, I need to get for all timesteps within the grib the pixel values, it seems impossible.
    Do you have any idea?

    Thanks!

    Posted from United States
    Google Chrome 5.0.375.55 Windows XP
  4. tomkralidis said,

    Wrote on June 30, 2010 @ 22:01:48

    Mozilla Firefox 3.6.6 Mac OS X 10

    I’m not exactly sure if you can do that without scripting or post-processing. Your best bet is to ask the mapserver-users mailing list.

    Posted from Canada Canada
    Mozilla Firefox 3.6.6 Mac OS X 10

Comment RSS · TrackBack URI

Leave a Comment

Name: (Required)

E-mail: (Required)

Website:

Comment:

Modified: 4 June 2010 19:43:09 EST