{"id":485,"date":"2010-06-04T19:02:07","date_gmt":"2010-06-04T23:02:07","guid":{"rendered":"http:\/\/www.kralidis.ca\/blog\/?p=485"},"modified":"2023-10-28T14:11:15","modified_gmt":"2023-10-28T18:11:15","slug":"displaying-grib-data-with-mapserver","status":"publish","type":"post","link":"https:\/\/www.kralidis.ca\/blog\/2010\/06\/04\/displaying-grib-data-with-mapserver\/","title":{"rendered":"Displaying GRIB data with MapServer"},"content":{"rendered":"

I recently had the opportunity to prototype WMS visualization of meteorological data<\/a>.\u00a0 MapServer, GDAL and Python to the rescue!\u00a0 Here are the steps I took to make it happen.<\/p>\n

Welcome to the captivating realm of luxury timepieces, where sophistication meets craftsmanship and exclusivity reigns supreme. In this guide, we will take you on an exhilarating journey through the mesmerizing world of Richard Mille replica<\/a> watches \u2013 a parallel universe where impeccable designs and intricate details intertwine with affordability. Brace yourself as we unveil the hidden gems that mimic these iconic timepieces flawlessly, alongside revealing the best-kept secrets on where to find these sought-after fakes. Whether you’re an aficionado looking to expand your collection or simply curious about the artistry behind high-end replicas, prepare to be captivated by a dazzling array of horological wonders that redefine what it truly means to own a prestigious Richard Mille watch.<\/p>\n

The data, (GRIB<\/a>), is a GDAL supported format<\/a>, so MapServer can handle processing as a result.\u00a0 The goal here was to create a LAYER<\/a> object.\u00a0 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).<\/p>\n

Here’s the hack:<\/p>\n

import sys\nimport osgeo.gdal as gdal\nimport osgeo.osr as osr\n\nif len(sys.argv) < 3:\n print 'Usage: %s <file> <numclasses>' % sys.argv[0]\n sys.exit(1)\n\ncvr = 256  # range of RGB values\nnumclasses = int(sys.argv[2])  # number of classifiers\n\nds = gdal.Open(sys.argv[1])\n\n# get proj4 def and write out PROJECTION object\np = osr.SpatialReference()\ns = p.ImportFromWkt(ds.GetProjection())\np2 = p.ExportToProj4().split()\n\nprint '\u00a0 PROJECTION'\nfor i in p2:\n print '\u00a0\u00a0 \"%s\"' % i.replace('+','')\nprint '\u00a0 END'\n\n# get band pixel data ranges and classify\nband = ds.GetRasterBand(1)\nmin = band.GetMinimum()\nmax = band.GetMaximum()\n\nif min is None or max is None:\u00a0 # compute automagically\n (min, max) = band.ComputeRasterMinMax(1)\n\n# calculate range of pixel values\npixel_value_range = float(max - min)\n# calculate the intervals of values based on classes specified\npixel_interval = pixel_value_range \/ numclasses\n# calculate the intervals of color values\ncolor_interval = (pixel_interval * cvr) \/ pixel_value_range\n\nfor i in range(numclasses):\n print '''\u00a0 CLASS\n  NAME \"%.2f to %.2f\"\n  EXPRESSION ([pixel] >= %.2f AND [pixel] < %.2f)\n  STYLE\n   COLOR %s %s %s\n  END\n END''' % (min, min+pixel_interval, min, min+pixel_interval, cvr, cvr, cvr)\n min += pixel_interval\n cvr -= int(color_interval)<\/pre>\n
Running this script outputs various bits for MapServer mapfile configuration.\u00a0 Passing more classes to the script creates more CLASS<\/a> objects, resulting in a smoother looking image.<\/p>\n
Here’s an example GetMap request:<\/p>\n
\"Meteorological<\/a><\/p>\n
Users can query to obtain pixel values (water temperature in this case) via GetFeatureInfo.\u00a0 Given that these are produced frequently, we can use the WMS GetMap TIME parameter to create time series maps of the models.<\/p>\n\n\n