path: root/ipl/gpacks/tiger/tgrmap.icn

############################################################################
#
#	File:     tgrmap.icn
#
#	Subject:  Program to generate map from TIGER files
#
#	Authors:  Gregg M. Townsend and William S. Evans
#
#	Date:     July 29, 2000
#
############################################################################
#
#  Tgrmap draws maps based on TIGER data files from the Census Bureau.
#  Data files must be in "line chain" (.lch) format as written by the
#  associated "tgrprep" program.
#
#  Usage:  tgrmap [file.lch ...]
#	Input is zero or more files of chains created by "tgrprep.icn".
#
#  All manipulation is done by mouse actions, keyboard shortcuts, or
#  window resizing.  There are no menus (although they would be nice
#  to have.)
#
#  Mouse actions:
#	Sweeping an area with the left mouse button zooms the image
#	to display that area better.  To cancel a sweep, just reduce
#	the swept height or width to less than 10 pixels.
#
#	Clicking the center mouse button pops up the Layers dialog,
#	which selects the categories of data to be shown or hidden.
#	No other actions are accepted while the dialog box is up.
#	(This only works on 8-bit displays, a vanishing breed.)
#
#	Clicking the right mouse button when the cursor is on a line
#	brings up a subwindow that shows the name and type of the
#	line.
#
#  Keyboard shortcuts:
#	F	find a feature by name
#	L	bring up the layers dialog
#	M	create PPM image
#	O	open a new file
#	P	create PostScript file for printing
#	Q	quit
#	R	refresh the display
#	S	save the displayed data to file
#	+	zoom in by a factor of 2
#	-	zoom out by a factor of 2
#	2-9	zoom to factor given
#	1	reset original map (centered and unzoomed)
#	arrow	arrow keys shift the center of the displayed area
#
#  Window resizing is allowed.
#
############################################################################
#
#  Requires:  Version 9 graphics
#
############################################################################
#
#  Links:  clipping, graphics, numbers, pscript, strings
#
############################################################################

#  Ideas for future changes:
#	Add menu alternatives to keyboard shortcuts
#	Write *color* PostScript, at least as an option
#	   (the programming is easy; tuning the colors is the hard part)

link clipping
link graphics
link numbers
link pscript
link strings

$include "keysyms.icn"

$ifndef _X_WINDOW_SYSTEM
   $define Key_KP_Up Key_Up
   $define Key_KP_Down Key_Down
   $define Key_KP_Left Key_Left
   $define Key_KP_Right Key_Right
$endif

$define MARGIN 5			# margin around full-sized map
$define CLIP 2				# clipping margin, allowing for linewdth
$define SHIFTBY 32			# number of pixels to shift at once
$define PSSCALE 100			# scaling from pixels to PostScript
$define MAXDRAW 4000			# maximum (even) args to avoid error 301
$define EPSILON 2.5			# how close is enough when clicking

# file values
global ifileList, fnameList		# input files and their names
global lonmin, lonmax, latmin, latmax 	# input range

# windows
global wintbl				# table of GCs by type
global msgwin				# base window for notices
global title				# window title
global tmpwin				# temp window for PPM snapshots

# window parameters
global dx, dy				# current translation values
global fullx, fully			# scaling for zoom-1 display

# display parameters
global ctrlon, ctrlat			# longitude/latitude of center
global curzoom, xscale, yscale		# current zoom factor and scaling
global lonrange, latrange		# distance from center to edge

# inquiry parameters
global litName				# string to match against feature name


# the classification list finds things based on line type
record class(				# classification record
   prefix,				#    CFCC prefix
   psgray,				#    PostScript graylevel
   pswidth,				#    PostScript linewidth
   label,				#    label
   color,				#    color
   width,				#    line width
   index,				#    mutable color index
   vispfx)				#    prefix code, but only if visible
global clist				# list of classification records
global ctable				# table keyed by two-char prefix



procedure main(args)
   local e, xywh, lon1, lat1, lon2, lat2, hilightOnly

   Window("size=870,870", "bg=pale moderate brown", args)

   &error := 1
   WAttrib("resize=on")
   &error := 0

   Font(("Helvetica" | "Univers" | "LucidaSans") || ",bold,12")	# may fail
   WAttrib("pointer=cross" | "pointer=crosshair")		# may fail
   msgwin := Clone(&window)		# save shaded window for notices
   setclasses()

   if *args > 0 then
      setfiles(args) | exit()
   else
      setfiles() | exit()

   setwindow()
   setcenter()
   setzoom()
   drawmap()

   repeat case e := Event() of {

      !"01": {
         setregion(lonmin, lonmax, latmin, latmax)
         drawmap()
         }

      !"23456789":	{ setzoom(e); drawmap() }
      !"+=":		{ setzoom(2.0 * curzoom); drawmap() }
      "-":		{ setzoom(0.5 * curzoom); drawmap() }

      !"Oo": {
         if setfiles() then {
            setwindow()
            setcenter()
            setzoom()
            drawmap()
            }
         }

      !"Ll":	setlayers()
      !"Rr":	drawmap()

      !"Mm": {
         tmpwin := WOpen("canvas=hidden",
            "width=" || WAttrib("width"), "height=" || WAttrib("height"))
         if /tmpwin then {
            Notice("can't open temporary canvas", "for PPM snapshot")
            break
            }
         CopyArea(&window, tmpwin)
         writefile(writeppm, "Write PPM file:", "Writing PPM file...")
         WClose(tmpwin)
         tmpwin := &null
         }

      !"Pp":	
         writefile(writeps, "Write PostScript file:", "Writing PostScript...")

      !"Ss":	
         writefile(writelch, "Save displayed portion as:", "Saving...")

      !"Ff": {
         if /litName then {
            hilightOnly := 1
            litName := ""
            }
         else {
            hilightOnly := &null
         }
         if TextDialog("Find features named:", , litName) == "Okay" then {
            litName := map(dialog_value[1])
            if litName == "" then litName := &null
            drawmap(hilightOnly)
         }
      }

      QuitEvents():	break

      Key_Left | Key_KP_Left:	shift(e, +SHIFTBY, 0)
      Key_Right | Key_KP_Right:	shift(e, -SHIFTBY, 0)
      Key_Up | Key_KP_Up:	shift(e, 0, +SHIFTBY)
      Key_Down | Key_KP_Down:	shift(e, 0, -SHIFTBY)

      &lpress: {
         xywh := Sweep()
         if xywh[3|4] < 10 then
            next
         lon1 := ctrlon + (get(xywh) - 0.5) / xscale
         lat1 := ctrlat + (get(xywh) - 0.5) / yscale
         lon2 := lon1 + (get(xywh) + 0.5) / xscale
         lat2 := lat1 + (get(xywh) + 0.5) / yscale
         setregion(lon1, lon2, lat1, lat2)
         drawmap()
         }

      &mrelease: {
         setlayers()
         }

      &rrelease: {
         identify(&x, &y)
         }

      &resize: {
         resize()
         }
   }
end


procedure writefile(proc, caption, message)
   local oname, ofile

   repeat case OpenDialog(msgwin, caption) of {
      "Okay": {
         if *dialog_value = 0 then
            next
         if close(open(oname := dialog_value)) then
            case TextDialog(msgwin, "Overwrite existing file?", , , ,
               ["Yes", "No", "Cancel"]) of {
                  "Yes":	&null
                  "No":		next
                  "Cancel":	fail
               }
         if ofile := open(oname, "w") then
            break
         case TextDialog(msgwin, "Cannot open " || oname) of {
            "Okay":	next
            "Cancel":	fail
            }
         }
      "Cancel":
         fail
      }

   Popup(msgwin, , , 32 + TextWidth(msgwin, message), 32,
      popmsg, message, proc, ofile)
   close(ofile)
   return
end

procedure popmsg(message, proc, ofile)
   CenterString(WAttrib("clipw") / 2, WAttrib("cliph") / 2, message)
   return proc(ofile)
end


procedure setfiles(L)
   local f, fname

   /L := list()
   fnameList := list()
   every close(!(\ifileList))
   ifileList := list()
   prescan()			# reset lonmin,lonmax,latmin,latmax
   until *fnameList > 0 do {
      until *L > 0 do {
         case OpenDialog(msgwin, "Input file(s):") of {
            "Okay":    put(L, words(dialog_value))
            "Cancel":  fail
         }
      }
      while fname := get(L) do {
         if not (f := open(fname)) then {
            Notice(msgwin, "Cannot open " || fname)
            next
         }
         if not (prescan(f)) then {
            Notice(msgwin, "Invalid format: " || fname)
            close(f)
            next
         }
         put(fnameList, fname)
         put(ifileList, f)
      }
   }
   return
end



#  prescan(f) -- verify that f is a valid file, setting globals if so

procedure prescan(f)
   local line, alon, alat, blon, blat
   if /f then {
      lonmin := latmin := 9999999
      lonmax := latmax := 0
      return
   }
   line := read(f)		| fail
   line ? {
      ="        "		| fail
      alon := move(7)		| fail
      alat := move(7)		| fail
   }
   line := read(f)		| fail
   line ? {
      ="        "		| fail
      blon := move(7)		| fail
      blat := move(7)		| fail
   }
   if alon > blon then {
      alon :=: blon
      alat :=: blat
   }
   lonmin >:= alon
   latmin >:= alat
   lonmax <:= blon
   latmax <:= blat
   return
end



procedure setwindow()
   local ww, wh, xstr, ystr, latsin, raspr, waspr

   ww := WAttrib("width")
   wh := WAttrib("height")
   dx := ww / 2
   dy := wh / 2

   xstr := "dx=" || (dx := WAttrib("width") / 2)
   ystr := "dy=" || (dy := WAttrib("height") / 2)
   every WAttrib(&window | !wintbl, xstr, ystr)

   # calculate aspect ratio of file region
   latsin := sin(((latmax + latmin) / 2.0) * (&pi / 9999999))
   raspr := real(lonmax - lonmin) / real(latmax - latmin) * latsin * (360 / 180)

   # calculate aspect ratio of window
   waspr := real(ww - 2 * MARGIN) / real(wh - 2 * MARGIN)

   # calculate scaling for zoom factor of 1.0
   if waspr > raspr then {
      # window is too wide
      fully := real(wh - 2 * MARGIN) / (latmax - latmin)
      fullx := fully * latsin * (360 / 180)
      }
   else {
      # window is too tall
      fullx := real(ww - 2 * MARGIN) / (lonmax - lonmin)
      fully := fullx / latsin / (360 / 180)
      }
   return
end



procedure setcenter(lon, lat)
   ctrlon := round(\lon | (lonmin + lonmax) / 2.0)
   ctrlat := round(\lat | (latmin + latmax) / 2.0)
   return
end



procedure setzoom(n)
   local x1, y1, x2, y2

   curzoom := \n | 1.0
   xscale := curzoom * fullx
   yscale := curzoom * fully
   lonrange := integer(dx / xscale + 0.5)
   latrange := integer(dy / yscale + 0.5)

   # clip out-of-bounds data because it's probably incomplete
   x1 := integer((lonmin - ctrlon) * xscale - 0.5) - CLIP
   x2 := integer((lonmax - ctrlon) * xscale + 0.5) + CLIP
   y1 := integer((latmin - ctrlat) * yscale - 0.5) - CLIP
   y2 := integer((latmax - ctrlat) * yscale + 0.5) + CLIP

   # limit clipping bounds to sensible values, else X gets confused
   x1 <:= -dx
   x2 >:= dx
   y1 <:= -dy
   y2 >:= dy

   # clip only drawing windows; NOT &window, used for copying and erasing!
   every Clip(!wintbl, x1, y1, x2 - x1, y2 - y1)
   return
end



procedure resize()
   local dxold, dyold, xshift, yshift

   dxold := dx			# save old translation values
   dyold := dy

   setwindow()			# set window parameters for new size

   xshift := dx - dxold
   yshift := dy - dyold

   # move to realign existing map with new window center
   CopyArea(-dx - xshift, -dy - yshift, 2 * dx, 2 * dy, -dx, -dy)
   if xshift > 0 then EraseArea(dx - xshift, -dy)
   if yshift > 0 then EraseArea(-dx, dy - yshift)

   # restore scaling and clipping
   setzoom(xscale / fullx)	# don't change zoom, but reset other globals

   return
end



procedure shift(e, nx, ny)

   while Pending()[1] === e do
      Event()				# consume duplicate shift events

   setcenter(ctrlon - nx / xscale, ctrlat - ny / yscale)
   CopyArea(-dx, -dy, 2 * dx, 2 * dy, -dx + nx, -dy + ny)
   if (nx > 0) then EraseArea(-dx, -dy, nx, 2 * dy)
   if (ny > 0) then EraseArea(-dx, -dy, 2 * dx, ny)
   if (nx < 0) then EraseArea(dx + nx, -dy)
   if (ny < 0) then EraseArea(-dx, dy + ny)
   settitle()			# reset center coords in title
   setzoom(curzoom)		# reset clipping

   drawmap()

   return
end



procedure drawmap(hilightOnly)
   local line, w, worig, lon, lat, dlon, dlat, a, bdy, dim, class, fename, f
   local drawProc, litFeature

   WAttrib("pointer=wait" | "pointer=watch")

   if /hilightOnly then {
      EraseArea()
      settitle()
   }
   litFeature := list()

   every f := !ifileList do {
      seek(f, 1)
      read(f)		# skip minima line
      read(f)		# skip maxima line

      while line := read(f) do line ? {

         if *Pending() > 0 then {
	    WAttrib("pointer=cross" | "pointer=crosshair")
	    return
	 }

         w := \wintbl[class := move(2)] | next
         move(1)
         bdy := move(1)
         if ="|" then {
            fename := tab(upto('|'))
            move(1)
         }
         else {
            fename := &null
         }
         dim := integer(move(4))
         lon := move(7) - ctrlon
         lat := move(7) - ctrlat

# quick clip
         if dim < 9999 & 
            (lon - dim > lonrange | lon + dim < -lonrange |
             lat - dim > latrange | lat + dim < -latrange) then
            next

         a := [xscale * lon, yscale * lat]
         while (dlon := move(4) - 5000) & (dlat := move(4) - 5000) do
            put(a, xscale * (lon +:= dlon), yscale * (lat +:= dlat))

# if beyond valid X range (with dx/dy margin), use library clipper
         if (!a > 32000) | (!a < -32000) then
            drawProc := DrawClipped
         else
            drawProc := DrawLine

         push(a, w)	# add graphics context

         if find(\litName, map(\fename)) then {
            put(litFeature, drawProc, a)
         }

         if /hilightOnly then {
            if any('57', bdy) & (a[1] := \wintbl["Y" || bdy]) then {
               drawProc ! a		# draw boundary indicator
               if any('F', class) then
                  next			# chain is ONLY a boundary
               a[1] := w
            }
            drawProc ! a			# draw line itself
         }
      }
   }
   if w := \wintbl["LL"] then {
      repeat {
         drawProc := get(litFeature) | break
         a := get(litFeature) | break
         a[1] := w	# replace graphics context
         drawProc ! a
      }
   }

   WAttrib("pointer=cross" | "pointer=crosshair")
   return
end



procedure identify(x, y)
   local line, lon, lat, dlon, dlat, dim, s, f
   local fename, cfcc, bndry, x0, y0, w, h
   local features

   WAttrib("pointer=wait" | "pointer=watch")

   # calculate region of interest in lat/lon coordinates
   x := (x - EPSILON) / xscale
   y := (y - EPSILON) / yscale
   w := (1 + 2 * EPSILON) / xscale
   h := (1 + 2 * EPSILON) / yscale

   features := set()
   every f := !ifileList do {
      seek(f, 1)
      read(f)		# skip minima line
      read(f)		# skip maxima line

      while line := read(f) do line ? {

         if *Pending() > 0 then {
	    WAttrib("pointer=cross" | "pointer=crosshair")
	    return
	 }

         cfcc := move(3)
         bndry := move(1)
         if ="|" then {		# get feature name
            fename := tab(upto('|'))
            move(1)
         }
         else {
            fename := ""
         }
         dim := integer(move(4))
         lon := move(7) - ctrlon
         lat := move(7) - ctrlat
         if dim < 9999 &
            (lon - dim > lonrange | lon + dim < -lonrange |
             lat - dim > latrange | lat + dim < -latrange) then
            next

         x0 := lon
         y0 := lat
         while (dlon := move(4) - 5000) & (dlat := move(4) - 5000) do {
            lon +:= dlon
            lat +:= dlat
            if ClipLine([x0, y0, lon, lat], x, y, w, h) then {
               s := case bndry of {
                  "9":" (national boundary) "
                  "8":" (state boundary) "
                  "7":" (county boundary) "
                  "5":" (city limit) "
                  "0":" "
               }
               insert(features, cfcc || s || fename)
            }
            x0 := lon
            y0 := lat
         }
      }
   }
   WAttrib("pointer=cross" | "pointer=crosshair")
   Popup(, , , WAttrib("leading") * (0 ~= *features), popList, sort(features))
   return
end



procedure popList(l)
   WAttrib("row=1", "col=1")
   every WWrite(!l)
   until Active()
end



procedure settitle()
   local lon, lat

   lon := ctrlon * (360.0 / 9999999)
   if lon > 180.0 then
      lon -:= 360.0
   lat := 90.0 - ctrlat * (180.0 / 9999999)
   title := fnameList[1]
   if *fnameList > 1 then
      title ||:= "..."
   title ||:= ":  " || dms(lon, "W", "E") || "  " || dms(lat, "S", "N")
   WAttrib("label=" || title)
   return
end



procedure dms(n, s1, s2)
   local deg, min, sec

   if n < 0 then
      n := -n
   else
      s1 := s2

   deg := integer(n)
   n := (n - deg) * 60
   min := integer(n)
   n := (n - min) * 60
   sec := integer(n + 0.5)

   return deg || "\260" || right(min, 2, "0") || "'" ||
      right(sec, 2, "0") || "\"" || s1
end



procedure setregion(lomin, lomax, ltmin, ltmax)
   local xzoom, yzoom

   setcenter((lomin + lomax + 1) / 2, (ltmin + ltmax + 1) / 2)
   xzoom := ((dx - MARGIN) * 2 / fullx) / (lomax - lomin)
   yzoom := ((dy - MARGIN) * 2 / fully) / (ltmax - ltmin)
   if xzoom < yzoom then
      setzoom(xzoom)
   else
      setzoom(yzoom)
   return
end



#  setclasses() -- initialize table of classifications
#
#  The order used here is reflected in the Layers dialog box.

procedure setclasses()
   local c, w, mcolors, m

   clist := [				# classification list
      # prefix, psgray&w, label,   color,  width
      class("A1", .0, 4, "roads", "black", 3),		# freeway/tollway
      class("A2", .0, 2, "roads", "black", 2),		# primary road
      class("A3", .0, 1, "roads", "black"),		# secondary road
      class("A4", .0, 0, "roads", "white"),		# local road
      class("A",  .0, 0, "roads", "white"),		# other road
      class("B1", .4, 2, "railroads", "deep green", 2),	# railroad line
      class("B",  .4, 1, "railroads", "deep green"),	# r.r. spur, yard, etc.
      class("H",  .7, 1, "water", "dark cyanish blue"),	# water
      class("Y7", .9, 5, "major boundaries", "orange", 3),	# county
      class("Y5", .9, 3, "major boundaries", "orange", 2),	# city
      class("E",  .9, 1, "minor boundaries", "light orange"),	# visible
      class("F",  .9, 1, "minor boundaries", "light orange"),	# invisible
      class("D",  .0, 1, "landmarks", "dark red"),		# landmark
      class("C",  .5, 1, "piplines & power", "purple"),		# pipe, power
      class("LL", .2, 2, "highlighted feature", "yellow", 10),	# hilit feature
      class("T0", .8, 3, "GPS track", "dark greenish cyan", 2),	# Track data
      class("X",  .8, 1, "unclassified", "purple")]		# unclassified

   every c := !clist do
      c.vispfx := c.prefix		# initially, all layers visible

   ctable := table()
   every c := !clist do
      if *c.prefix = 1 then
         every /ctable[c.prefix || !"0123456789"] := c
      else
         ctable[c.prefix] := c

   wintbl := table()			# global window table
   mcolors := table()			# local table of mutable colors

   every c := !clist do {

      w := Clone() | stop("can't clone window for ", c.label)
      /mcolors[c.color] := NewColor(w, c.color)		# may fail
      c.index := mcolors[c.color]
      Fg(w, \c.index | c.color) | stop("can't set color for ", c.label)

      WAttrib(w, "linewidth=" || \c.width)
      wintbl[c.prefix] := w
      if *c.prefix = 1 then
         every /wintbl[c.prefix || (0 to 9)] := w
      }
   return
end



#  setlayers() -- bring up layers dialog

procedure setlayers()
   local c, i, defaults, buttons, choice, lset
   static labels, values

   initial {
      lset := set()
      labels := list()
      values := list()
      every c := !clist do
         if \c.index & not member(lset, c.label) then {
            insert(lset, c.label)
            put(labels, c.label)
            put(values, 1)
            }
      }

   if *labels = 0 then {
      Notice("No layer control available")
      fail
      }

   while choice ~=== "Okay" do {	# loop when "Apply" selected

      defaults := values
      buttons := ["Okay", "Apply", "Cancel"]
      choice := ToggleDialog(msgwin, "Layers:", labels, defaults, buttons)
      if choice == "Cancel" then
         fail
      values := dialog_value

      #  change mutable color for every item that changed in the dialog
      every i := 1 to *values do
         if values[i] ~=== defaults[i] then
            every c := !clist do
               if c.label == labels[i] then {
                  if \values[i] then {
                     Color(\c.index, c.color)
                     c.vispfx := c.prefix
                     }
                  else {
                     Color(\c.index, Bg())
                     c.vispfx := &null
                     }
                  }
      }
   return
end



procedure writelch(ofile)
   local line, dim, lon, lat, f, a, b, x, y, dlon, dlat, nlon, nlat, w, head
   local startlon, startlat, minlon, minlat, maxlon, maxlat, deltas, class

   write(ofile, "        ",
      right(ctrlon - lonrange, 7), right(ctrlat - latrange, 7))
   write(ofile, "        ",
      right(ctrlon + lonrange, 7), right(ctrlat + latrange, 7))

   every f := !ifileList do {
      seek(f, 1)
      read(f)			# skip minima line
      read(f)			# skip maxima line

      while line := read(f) do line ? {
         w := \wintbl[class := move(2)] | next
         head := class || move(2)
         if ="|" then {
            head ||:= "|" || tab(upto('|') + 1)
         }
         dim := integer(move(4))
         lon := move(7) - ctrlon
         lat := move(7) - ctrlat
# quick clip
         if dim < 9999 & 
            (lon - dim > lonrange | lon + dim < -lonrange |
             lat - dim > latrange | lat + dim < -latrange) then
            next

         a := [xscale * lon, yscale * lat]
         while (dlon := move(4) - 5000) & (dlat := move(4) - 5000) do
            put(a, xscale * (lon +:= dlon), yscale * (lat +:= dlat))
         a := Coalesce(ClipLine(w, a)) | next
         every b := !a do {
            deltas := ""
            startlon := minlon := maxlon := lon :=
               round(get(b) / xscale) + ctrlon
            startlat := minlat := maxlat := lat :=
               round(get(b) / yscale) + ctrlat
            while nlon := round(get(b) / xscale) + ctrlon do {
               nlat := round(get(b) / yscale) + ctrlat
               deltas ||:= right(nlon - lon + 5000, 4, "0")
               deltas ||:= right(nlat - lat + 5000, 4, "0")
               lon := nlon
               lat := nlat
               maxlon <:= lon
               minlon >:= lon
               maxlat <:= lat
               minlat >:= lat
            }
            dim := startlon - minlon
            dim <:= maxlon - startlon
            dim <:= startlat - minlat
            dim <:= maxlat - startlat
            dim >:= 9999

            write(ofile, head, right(dim, 4), right(startlon, 7, "0"),
                  right(startlat, 7, "0"), deltas)
         }
      }
   }
   return
end

#  writeppm(ofile) -- write PPM image to ofile
#
#  comments note latitude and longitude bounds in arc-seconds

procedure writeppm(ofile)
   local w, h, rw, rh, s, lon, lat, dlon, dlat

   w := WAttrib("width")
   h := WAttrib("height")
   rw := real(w)
   rh := real(h)

   lon := ctrlon * (360.0 / 9999999)
   if lon > 180.0 then
      lon -:= 360.0
   lat := 90.0 - ctrlat * (180.0 / 9999999)

   dlon := lonrange * (360.0 / 9999999)
   dlat := latrange * (180.0 / 9999999)

   write(ofile, "P6")
   write(ofile, "#RTIN")
   write(ofile, "#lon,lat:",
      arcs(lon - dlon), arcs(lat - dlat), arcs(lon - dlon), arcs(lat + dlat),
      arcs(lon + dlon), arcs(lat + dlat), arcs(lon + dlon), arcs(lat - dlat))
   write(ofile, "#x,y: 0.0 0.0 0.0 ", rh, " ", rw, " ", rh, " ", rw, " 0.0")
   write(ofile, w, " ", h)
   write(ofile, 255)

   every writes(ofile, rgb24(Pixel(tmpwin)))

   return
end



#  arcs(n) -- format latitude or longitude in arc-seconds with leading space

procedure arcs(n)
   return " " || (n * 3600.0)
end



#  rgb24(k) -- return 24-bit (3-byte) r-g-b value for k

procedure rgb24(k)
   local s, r, g, b
   static t
   initial t := table()

   if s := \t[k] then
      return s

   (ColorValue(k | Color(k)) | fail) ? {
      s := char(tab(upto(',')) / 256)
      move(1)
      s ||:= char(tab(upto(',')) / 256)
      move(1)
      s ||:= char(tab(0) / 256)
      }
   t[k] := s
   return s
end



#  writeps(ofile) -- write Encapsulated PostScript to ofile

procedure writeps(ofile)
   local x1, x2, y1, y2, line, prevcode, code, dim, lon, lat, a, n, c, f

   x1 := integer((lonmin - ctrlon) * xscale - 0.5) - CLIP
   y1 := integer((latmin - ctrlat) * yscale - 0.5) - CLIP
   x2 := integer((lonmax - ctrlon) * xscale + 0.5) + CLIP
   y2 := integer((latmax - ctrlat) * yscale + 0.5) + CLIP
   x1 <:= -dx
   y1 <:= -dy
   x2 >:= dx
   y2 >:= dy
   epsheader(ofile, (dx + x1) * PSSCALE, (dy - y2) * PSSCALE,
      (x2 - x1) * PSSCALE, (y2 - y1) * PSSCALE, "r")

   every write(ofile, ![
      "/m { moveto } bind def",
      "/r { rlineto } bind def",
      "/s { stroke } bind def",
      "/w { .00666667 mul inch setlinewidth setgray stroke } bind def"])

   every c := !clist do
      write(ofile, "/", left(\c.vispfx, 2),
         " { ", c.psgray, " ", c.pswidth, " w } bind def")

   every f := !ifileList do {
      seek(f, 1)
      read(f)			# skip minima line
      read(f)			# skip maxima line

      while line := read(f) do line ? {
         code := \ctable[move(2)].vispfx | next
         move(2)
# skip feature name
         if ="|" then
            tab(upto('|') + 1)
         dim := integer(move(4))
         lon := xscale * (move(7) - ctrlon)
         lat := yscale * (move(7) - ctrlat)
         if dim < 9999 &
            (lon - dim > lonrange | lon + dim < -lonrange |
             lat - dim > latrange | lat + dim < -latrange) then
            next
         else {
            writes(ofile, integer(PSSCALE * (dx + lon)), " ",
               integer(PSSCALE * (dy - lat)), " m")
            while lon := move(4) - 5000 & lat := move(4) - 5000 do
               writes(ofile, "\n", integer(PSSCALE * xscale * lon), " ",
                  integer(-PSSCALE * yscale * lat), " r")
            write(ofile, " ", (prevcode ~===:= code) | "s")
         }
      }
   }
   write(ofile, "showpage")
   return
end