I remembered this interesting post from the MapServer list back in
July. I believe that they still ended up with lines and not polys, but
their methods may inspire a creative solution...<br>
<br>
_________________<br>
<br>
<pre id="nonprop"><span id="MSGHDR-PRE"><span id="MSGHDR-DATE-H-PRE"><b>Date:</b></span> <span id="MSGHDR-DATE-PRE">Fri, 22 Jul 2005 22:05:54 -0800</span><br><span id="MSGHDR-REPLY-TO-H-PRE"><b>Reply-To:</b></span>
<span id="MSGHDR-REPLY-TO-PRE">Rob Cermak <<a href="mailto:cermak@SFOS.UAF.EDU">cermak@SFOS.UAF.EDU</a>></span><br><span id="MSGHDR-SENDER-H-PRE"><b>Sender:</b></span> <span id="MSGHDR-SENDER-PRE">UMN MapServer Users List <
<a href="mailto:MAPSERVER-USERS@LISTS.UMN.EDU">MAPSERVER-USERS@LISTS.UMN.EDU</a>></span><br><span id="MSGHDR-FROM-H-PRE"><b>From:</b></span> <span id="MSGHDR-FROM-PRE">Rob Cermak <<a href="mailto:cermak@SFOS.UAF.EDU">
cermak@SFOS.UAF.EDU</a>></span><br><span id="MSGHDR-SUBJECT-H-PRE"><b>Subject:</b></span> <span id="MSGHDR-SUBJECT-PRE">MapHack for (Weather) data : Contouring (on-the-fly)</span><br></span><span id="MSGHDR-CONTENT-TYPE-H-PRE">
<b>Content-Type:</b></span> <span id="MSGHDR-CONTENT-TYPE-PRE">TEXT/PLAIN; charset=US-ASCII</span><p><br>Just happened to stumble on this quite by accident. We have been looking <br>for a way to dynamically contour data for quite some time. Googling
<br>surfaces a lot of people attempting various means. Mapping Hacks shows <br>how it can be done using GRASS.<br><br>Most OOS (Ocean Observing System) mapserver sites can plot data rather<br>easily...color filled circles for sea surface temperature... here is
<br>something to help you go the next step. I present a bare bones method to<br>preparing data for contouring any data you've got. Making it dynamic is<br>left as an exercise for the reader.<br><br>TOOLS:<br> Data (Longitude, Latitude, Value[Z])
<br> GMT (<a href="http://gmt.soest.hawaii.edu/">http://gmt.soest.hawaii.edu/</a>)<br> - compiled with NetCDF<br> - blockmean<br> - surface<br> - Unix util: awk<br> gdal_contour (from GDAL 1.2.5)<br> - must have: GMT (rw): GMT NetCDF Grid Format (as seen by --formats)
<br><br>Step 1: Get your data -- many many ways to do this. In this example, we <br>have long, lat, and sea level pressure(Z). The Z values can be anything.<br><br>170.000000,50.000000,1026.400000<br>170.500000,50.000000
,1026.320000<br>171.000000,50.000000,1026.210000<br>171.500000,50.000000,1026.080000<br>172.000000,50.000000,1025.910000<br>172.500000,50.000000,1025.720000<br>173.000000,50.000000,1025.510000<br>173.500000,50.000000,1025.260000
<br>174.000000,50.000000,1025.000000<br>174.500000,50.000000,1024.710000<br>[..snip..]<br><br>Step 2: Use GMT to convert the data to a GMT NetCDF gridded file<br><br># You will need to set the GRID REGION up properly for your use...
<br>GRID_REGION="-180/180/30/80";<br># blockmean option -I put all our nice data into 0.5 degree bins for <br># the surface program<br>awk -F, '{if ($3>0) printf "%f %f %f\n",$1,$2,$3}' ak.psl | blockmean -R$GRID_REGION -
I0.5 > ak.grd<br># Take output from blockmean and grid it into a 1/2 by 1/2 degree grid and <br># stuff info a GMT (netcdf) grid file<br>surface ak.grd -R$GRID_REGION -I0.5 -Gak2.grd<br><br>The normal next step is to plot it up using a GMT tool that produces
<br>postscript! Ick!<br><br>Step 3: Abuse of GDAL - gdal_contour<br><br>Gdal is really for operating on raster data, but with the GMT Netcdf Grid<br>driver, we can abuse it. Taking our above data, treat it like a raster or
<br>heights from a DEM. Now convert from GMT grid to shapefile.<br><br>gdal_contour -a elev ak2.grd ak2s.shp -i 1.0<br><br>From the shape file, you can import into PostGIS or go right to the <br>display in mapserver.<br>
<br> layer<br> GROUP slp<br> NAME something<br> DEBUG ON<br> DATA ak2s<br> TYPE LINE<br> STATUS ON<br> class<br> color 200 200 0<br> end<br> end<br><br> layer<br> GROUP slp<br> NAME label_slp
<br> DATA ak2s<br> LABELITEM elev<br> TYPE ANNOTATION<br> CLASS<br> LABEL<br> POSITION CC<br> SIZE TINY<br> COLOR 0 0 0 <br> END<br> END<br> end<br><br>So, that is the basic hack. To do:
<br><br>* allow the user to control the contour interval, color, ...<br> - allow variation of contours - dashed lines?<br>* control annotation the values 1004 vs 1004.000 (cleanup through PostGIS)<br>* properly identify the values -- its not just heights anymore...
<br><br>Now you can contour any parameter you want: sea level pressure, heights,<br>temperature, .... or about anything gridded in the oceanographic<br>/meteorologic world. For all the weather winnies of the world...<br>
<br>Sample image: Sea level pressure map as produced from NCEP's<br> Global Forecast System on a 1/2 deg x 1/2 deg grid.<br> <a href="http://ak.aoos.org/mapserver/aoos/slp.png">http://ak.aoos.org/mapserver/aoos/slp.png</a>
<br><br>If there is an easier way, please send it my way.<br><br>Happy Maphacking!<br>Rob<br>-- <br>Rob Cermak : 907-474-7948 : FAX 907-474-7204 : PGP = 0x75869A6E<br>Alaska Ocean Observing System : Data Management <br>School of Fisheries and Ocean Sciences : University of Alaska Fairbanks
<br><a href="mailto:cermak@sfos.uaf.edu">cermak@sfos.uaf.edu</a> : <a href="mailto:fnjrc1@uaf.edu">fnjrc1@uaf.edu</a><br><a href="mailto:ED9U1M7P01@alaska.edu">ED9U1M7P01@alaska.edu</a> : <a href="mailto:cermak@alaska.edu">
cermak@alaska.edu</a><br></p></pre>
<br><br><div><span class="gmail_quote">On 11/10/05, <b class="gmail_sendername">Frank Warmerdam</b> <<a href="mailto:warmerdam@pobox.com">warmerdam@pobox.com</a>> wrote:</span><blockquote class="gmail_quote" style="border-left: 1px solid rgb(204, 204, 204); margin: 0pt 0pt 0pt 0.8ex; padding-left: 1ex;">
On 11/10/05, Bryan Keith <<a href="mailto:bryan@geomega.com">bryan@geomega.com</a>> wrote:<br>> I use gdal_contour to create contours as lines. Does anyone know of a<br>> similar tool or have some code that will create polygons instead of
<br>> lines when contouring a surface? Each polygon would represent the area<br>> between n1 and n2. I have the ConRec code and thought about writing<br>> this myself, but it gets complicated with donut polygons. I bet
<br>> someone's done this before. Maybe there's a better mailing list for<br>> this question? Thanks.<br><br>Bryan,<br><br>An interesting problem, but I don't know the answer. It seems<br>like you should be able to use the contours, along with introduced
<br>edged lines and feed it into a topological analyser to build<br>polygons. You might want to look at GRASS to do this.<br><br>You might want to try the gislist, or freegis mailing list. Even<br>gdal-dev might be worth a try.
<br><br>Best regards,<br>--<br>---------------------------------------+--------------------------------------<br>I set the clouds in motion - turn up | Frank Warmerdam, <a href="mailto:warmerdam@pobox.com">warmerdam@pobox.com
</a><br>light and sound - activate the windows | <a href="http://pobox.com/~warmerdam">http://pobox.com/~warmerdam</a><br>and watch the world go round - Rush | Geospatial Programmer for Rent<br><br>_______________________________________________
<br>FWTools mailing list<br><a href="mailto:FWTools@lists.maptools.org">FWTools@lists.maptools.org</a><br><a href="http://lists.maptools.org/mailman/listinfo/fwtools">http://lists.maptools.org/mailman/listinfo/fwtools</a>
<br><a href="http://fwtools.maptools.org/">http://fwtools.maptools.org/</a><br></blockquote></div><br>