MapHack for (Weather) data : Contouring (on-the-fly)

[Rob Cermak]

Just happened to stumble on this quite by accident.  We have been looking 
for a way to dynamically contour data for quite some time.   Googling 
surfaces a lot of people attempting various means.  Mapping Hacks shows 
how it can be done using GRASS.

Most OOS (Ocean Observing System) mapserver sites can plot data rather
easily...color filled circles for sea surface temperature... here is
something to help you go the next step.  I present a bare bones method to
preparing data for contouring any data you've got.  Making it dynamic is
left as an exercise for the reader.

TOOLS:
 Data (Longitude, Latitude, Value[Z])
 GMT (http://gmt.soest.hawaii.edu/)
  - compiled with NetCDF
  - blockmean
  - surface
  - Unix util: awk
 gdal_contour (from GDAL 1.2.5)
  - must have: GMT (rw): GMT NetCDF Grid Format (as seen by --formats)

Step 1: Get your data -- many many ways to do this.  In this example, we 
have long, lat, and sea level pressure(Z).  The Z values can be anything.

170.000000,50.000000,1026.400000
170.500000,50.000000,1026.320000
171.000000,50.000000,1026.210000
171.500000,50.000000,1026.080000
172.000000,50.000000,1025.910000
172.500000,50.000000,1025.720000
173.000000,50.000000,1025.510000
173.500000,50.000000,1025.260000
174.000000,50.000000,1025.000000
174.500000,50.000000,1024.710000
[..snip..]

Step 2: Use GMT to convert the data to a GMT NetCDF gridded file

# You will need to set the GRID REGION up properly for your use...
GRID_REGION="-180/180/30/80";
# blockmean option -I put all our nice data into 0.5 degree bins for 
# the surface program
awk -F, '{if ($3>0) printf "%f %f %f\n",$1,$2,$3}' ak.psl | blockmean -R$GRID_REGION -I0.5 > ak.grd
# Take output from blockmean and grid it into a 1/2 by 1/2 degree grid and 
# stuff info a GMT (netcdf) grid file
surface ak.grd -R$GRID_REGION -I0.5 -Gak2.grd

The normal next step is to plot it up using a GMT tool that produces 
postscript!  Ick!

Step 3: Abuse of GDAL - gdal_contour

Gdal is really for operating on raster data, but with the GMT Netcdf Grid
driver, we can abuse it.  Taking our above data, treat it like a raster or
heights from a DEM.  Now convert from GMT grid to shapefile.

gdal_contour -a elev ak2.grd ak2s.shp -i 1.0