[FWTools] contours as polygons

Thu Nov 10 12:04:21 EST 2005

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...

_________________

*Date:*         Fri, 22 Jul 2005 22:05:54 -0800
*Reply-To:*     Rob Cermak <cermak at SFOS.UAF.EDU>
*Sender:*       UMN MapServer Users List <MAPSERVER-USERS at LISTS.UMN.EDU>
*From:*         Rob Cermak <cermak at SFOS.UAF.EDU>
*Subject:*      MapHack for (Weather) data : Contouring (on-the-fly)
*Content-Type:* TEXT/PLAIN; charset=US-ASCII

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

>From the shape file, you can import into PostGIS or go right to the
display in mapserver.

  layer
    GROUP slp
    NAME something
    DEBUG ON
    DATA ak2s
    TYPE LINE
    STATUS ON
    class
      color 200 200 0
    end
  end

  layer
    GROUP slp
    NAME label_slp
    DATA ak2s
    LABELITEM elev
    TYPE ANNOTATION
    CLASS
      LABEL
        POSITION CC
        SIZE TINY
        COLOR 0 0 0
      END
    END
  end

So, that is the basic hack.  To do:

* allow the user to control the contour interval, color, ...
  - allow variation of contours - dashed lines?
* control annotation the values 1004 vs 1004.000 (cleanup through PostGIS)
* properly identify the values -- its not just heights anymore...

Now you can contour any parameter you want:  sea level pressure, heights,
temperature, ....  or about anything gridded in the oceanographic
/meteorologic world.  For all the weather winnies of the world...

Sample image: Sea level pressure map as produced from NCEP's
  Global Forecast System on a 1/2 deg x 1/2 deg grid.
  http://ak.aoos.org/mapserver/aoos/slp.png

If there is an easier way, please send it my way.

Happy Maphacking!
Rob
--
Rob Cermak : 907-474-7948 : FAX 907-474-7204 : PGP = 0x75869A6E
Alaska Ocean Observing System : Data Management
School of Fisheries and Ocean Sciences : University of Alaska Fairbanks
cermak at sfos.uaf.edu : fnjrc1 at uaf.edu
ED9U1M7P01 at alaska.edu : cermak at alaska.edu

On 11/10/05, Frank Warmerdam <warmerdam at pobox.com> wrote:
>
> On 11/10/05, Bryan Keith <bryan at geomega.com> wrote:
> > I use gdal_contour to create contours as lines. Does anyone know of a
> > similar tool or have some code that will create polygons instead of
> > lines when contouring a surface? Each polygon would represent the area
> > between n1 and n2. I have the ConRec code and thought about writing
> > this myself, but it gets complicated with donut polygons. I bet
> > someone's done this before. Maybe there's a better mailing list for
> > this question? Thanks.
>
> Bryan,
>
> An interesting problem, but I don't know the answer. It seems
> like you should be able to use the contours, along with introduced
> edged lines and feed it into a topological analyser to build
> polygons. You might want to look at GRASS to do this.
>
> You might want to try the gislist, or freegis mailing list. Even
> gdal-dev might be worth a try.
>
> Best regards,
> --
>
> ---------------------------------------+--------------------------------------
> I set the clouds in motion - turn up | Frank Warmerdam,
> warmerdam at pobox.com
> light and sound - activate the windows | http://pobox.com/~warmerdam
> and watch the world go round - Rush | Geospatial Programmer for Rent
>
> _______________________________________________
> FWTools mailing list
> FWTools at lists.maptools.org
> http://lists.maptools.org/mailman/listinfo/fwtools
> http://fwtools.maptools.org/
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://lists.maptools.org/pipermail/fwtools/attachments/20051110/9a6abe5f/attachment.html