diff options
Diffstat (limited to 'ipl/gprogs/dlgvu.icn')
| -rw-r--r-- | ipl/gprogs/dlgvu.icn | 1900 |
1 files changed, 1900 insertions, 0 deletions
diff --git a/ipl/gprogs/dlgvu.icn b/ipl/gprogs/dlgvu.icn new file mode 100644 index 0000000..63a933d --- /dev/null +++ b/ipl/gprogs/dlgvu.icn @@ -0,0 +1,1900 @@ +############################################################################ +# +# File: dlgvu.icn +# +# Subject: Program to display USGS DLG map files +# +# Authors: Gregg M. Townsend and William S. Evans +# +# Date: October 2, 2005 +# +############################################################################ +# +# This file is in the public domain. +# +############################################################################ +# +# Contributor: Frank Glandorf +# +############################################################################ +# +# Dlgvu displays and prints USGS digital map data. +# +# usage: dlgvu [options] file... +# +# Each file argument is one of: +# a directory containing DLG files in SDTS format +# a ZIP format archive of such files (requires "unzip" utility) +# a text file containing coordinates of paths and features +# +# All interaction is via mouse actions and keyboard shortcuts. +# The display window may be resized as desired. +# +############################################################################ +# +# Command options: +# +# -c display coverage map only, without loading data +# -d print some debugging data +# -n no display; just report statistics, and then quit +# (this still requires X-Windows, unfortunately) +# -p use slower, more precise coordinate conversion +# -q quiet mode: no commentary to stdout +# -t load only maps traversed by paths, ignoring others +# -o logfile specify output log to use instead of standard output +# +# -l abcd display only layers a, b, c, d +# -x abcd exclude layers a, b, c, d +# +# For -l and -x, the following layer codes are used. +# (USGS equivalents are given in parentheses.) +# +# b boundaries (BD: boundaries) +# c contour lines (HP: hypsography) +# d sand, gravel, lava (NV: nonvegetative features) +# f feature labels read from text files +# g GPS paths read from text files +# l land sections (PL: public lands) +# m markers (SM: survey markers) +# n file names +# o other unknown layers from non-DLG data +# r roads (RD: roads) +# s structures (MS: manmade structures) +# t train tracks (RR: railroads) +# u utilities (MT: miscellaneous transportation) +# v vegetation (SC: surface cover) +# w water (HY: hydrology) +# +# Additionally, the standard Window() options are accepted; +# in particular, "-F fgcolor" sets the color used for drawing +# unrecognized ("other") layers, as from USGS "National Atlas" +# digital files, and "-G 800x500" sets the initial window size. +# +# Typical usage is simply +# dlgvu dir1 [dir2 ...] +# to display one or more adjacent maps. The -x option can speed +# things up by excluding unwanted layers; the contour layer is +# especially slow. +# +# A ZIP archive can replace a directory name if Icon can open +# the unzip program via a pipe. For example: +# dlgvu woodside.zip palo_alto.zip +# +############################################################################ +# +# Mouse actions: +# +# To zoom to a particular region, sweep out the region using the +# left mouse button. To cancel a sweep, reduce its width or height +# to fewer than ten pixels. +# +# If nothing appears to be happening after zooming in, the program +# is probably drawing offscreen. It's not smart about that. Be +# patient, and it will soon display the visible region. +# +# To display the latitude / longitude of a location, and a scale bar, +# hold down the right mouse button. +# +# To record a labeled feature, shift-click the left mouse button. +# Enter a name in the pop-up dialog box. The location and name are +# written to the log file and added to the feature layer of the map. +# +# To record an anonymous location to the log file, shift-click with +# the right mouse button instead. No dialog box appears. A sequence +# of anonymous locations can be read as a path by a subsequent +# program run. +# +############################################################################ +# +# Keyboard actions: +# +# + or = zoom in +# - or _ zoom out +# 0 or Home zoom to initial view +# arrow keys pan the display (hold Shift key for smaller pan) +# +# b, c, d, etc. toggle display of specified layer +# a display all loaded layers including n (file names) +# x display no layers (just an empty window) +# +# Esc stop drawing (any unrecognized key does this) +# space or Enter redraw screen (e.g. after inadvertent interrupt) +# q quit +# +# p or PrntScrn print visible portion to PostScript file +# +# The file produced by PrntScrn is an Encapsulated PostScript file +# suitable either for direct printing ("lpr file.ps") or for import +# into another document. +# +############################################################################ +# +# Input files: +# +# In directories and archives, only files with names ending in .ddf +# or .DDF are read; others are ignored. These files must be in SDTS +# (Spatial Data Transfer Standard) format, which is used by the USGS +# for all new DLG files. +# +# Text files supply coordinates for features or paths. GPS receivers +# are one possible source for such data. A text file can supply +# paths, features, or both. +# +# Paths are specified by sequences of lines that end with two decimal +# values. The values are interpreted as latitude and longitude, in +# that order. An interruption in the sequence (such as a blank line) +# indicates a break between paths. +# +# Features, or waypoints, are given by lines that *begin* with two +# decimal values. The rest of the line is taken as a label, which +# must not be empty and must not end with two decimal values. +# +# Any other line in a text file breaks a path sequence but is +# otherwise ignored. +# +############################################################################ +# +# About DLG files: +# +# Dlgvu was written to display digitized topographic maps produced +# by the United States Geological Survey (USGS). The current file +# format is based on the Spatial Data Transfer Standard (SDTS). +# Some older files are available in other formats (including +# "standard" and "optional") not supported by this program. +# +# DLG files are available free from the USGS at this web page: +# http://edc.usgs.gov/doc/edchome/ndcdb/ndcdb.html +# Coverage is incomplete. 24K maps, the most detailed, are available +# for only some areas, and many maps lack some of the data layers. +# +# Each map is represented by a collection of gzipped tar files +# (one for each map layer) that are unpacked for display. Multiple +# files versions may be available, and not all layers are available +# for all maps. +# +# IMPORTANT: Do not blindly unpack all the tar files of a map into +# the same directory; due to the use of duplicate file names in the +# transportation layers, some files will be overwritten. Instead, +# unpack the roads, railroads, and miscellaneous transportation +# layers separately, each time renaming the TR*.DDF files to RD*.DDF, +# RR*.DDF, and MT*.DDF respectively. +# +# Dlgvu has mainly been tested and tuned using "large scale" DLG +# files (1:24000 scale, covering 7.5 minute quadrangles). Other +# scales produce less attractive displays, partly due to different +# encodings: For example, the same residential streets may be encoded +# as "Class 3 Roads" in 24K files but "Class 4 Roads" in 100K files. +# +# Dlgvu does not presume to understand ISO 8211, DDF, STDS, and TVP +# in their full complexity and generality. Undoubtedly it is making +# some unwarranted assumptions based on observed practice. The file +# renaming recommended above is contrary to specification but allows +# a complete map to be stored in a single flat directory. +# +# For more information, and some sample data files, visit: +# http://www.cs.arizona.edu/icon/oddsends/dlgvu/ +# +############################################################################ +# +# Displayed features: +# +# DLG files are rich in detail. Dlgvu displays only some of this +# encoded data. +# +# Put simply, dlgvu understands point and line features but not +# area features. It draws a small square for a structure location, +# or draws the outline of a large building, but it does not color in +# an "urban area" in which individual structures are not plotted. +# It displays the shoreline of a river, and any islands, but does +# not understand enough to color the river area itself blue. +# +# Dlgvu recognizes some line features for special display. For +# example, major roads are drawn with wide lines, and trails are +# drawn with dashed red lines. Lines with unrecognized attributes, +# or no attributes, are drawn in a default style. Point features +# ("school", "windmill", etc.) are not distinguished. +# +# Area features are drawn only in outline. The most obvious of +# these are vegetated areas and urban areas. Land section and +# civil boundaries also delimit area features. +# +# Colors are assigned as follows (with layer codes on the left): +# +# b boundaries gold +# c contour lines tan +# f feature labels black +# g GPS path bold pink over "highlighter" +# l land sections pale red +# m survey markers blue +# n file names green +# o other data brown (can override with -F option) +# r roads, class 1-3 black or dark gray +# r roads, class 4-5 dashed dark gray +# r trails dashed red +# s structures brownish gray +# t railroads rust +# t rapid transit rails dark blue +# u pipelines dashed purple +# u power lines purple +# u airport runways gray +# v vegetation light green +# w water light blue +# x sand, gravel, lava greenish gray +# +# Dlgvu uses a simple rectangular projection that is satisfactory +# for small areas like 24K quadrangles but less suitable for large +# areas such as whole states. +# +############################################################################ +# +# The loading process: +# +# Data is loaded in two phases. A quick preloading phase determines +# the available layers and the geographic areas covered. A status +# line is output for each file. For example: +# +# bcl-r-tu-w N27 15 C66 42a 93a ia/ames-w +# +# The first field shows which layers were found. N27 declares that +# coordinates use the NAD 1927 geodetic datum; N83 for NAD 1983 is +# the other likely value. 15 is the UTM zone number; state maps with +# latitude/longitude data show "LL" here. C66 means that the data +# appears to have been projected using the Clarke 1866 ellipsoid; the +# other likely value is "G80" for the GRS 1980 ellipsoid. Dlgvu uses +# this to infer the datum, because the declared datum value is a less +# reliable indicator. +# +# "42a 93a" gives the coordinates of the southeast corner of the map, +# in degrees North and West, with letters "a" through "h" indicating +# fractions from 0/8 through 7/8. The final field is the file name. +# +# If the layers in a file are inconsistent (for example, in the +# inferred ellipsoid), multiple lines appear with a "*" prefix. +# If display of a file is suppressed by the "-t" option, an X +# prefixes the line. +# +# For text files, a notation such as "3:489+0" replaces layer +# indicators, counting continuous segments, total points, and +# feature labels. Coordinate values are assumed to use the +# WGS 1984 ("W84") datum and ellipsoid. +# +# The display appears during the longer second loading phase. For +# each layer of each input file, bounds are drawn and a progress bar +# changes as data is read. The color of the label indicates the +# layer being loaded. +# +############################################################################ +# +# Tiling multiple maps: +# +# Multiple maps are displayed in proper relation. To quickly see +# how the maps of a set will join, use "dlgvu -c". +# +# Small gaps or overlaps occasionally appear along boundaries when +# maps are tiled; these are symptomatic of inconsistent datums, and +# they reflect the true relationships of the maps to the earth. +# +# Dlgvu loads all necessary data into memory, so there is a very +# real limit to the amount of data that can be displayed. Contour +# lines, especially, take lots of memory, but they can be excluded +# by calling "dlgvu -xc". A 128MB Linux system can typically +# display three to five complex 24K quadrangles simultaneously +# without thrashing. +# +############################################################################ +# +# Known problems: +# +# On Linux, we have seen occasional crashes of the XFree86 server, +# especially under conditions of tight memory and/or extreme zooming. +# +# Colors on printed maps vary somewhat from those seen onscreen, +# depending on the printer. Printed maps do not include the "n" +# (file name) layer. +# +# While data is being loaded from a ZIP file, interrupting with ^Z +# can disrupt the "unzip" pipe and cause the program to crash or to +# display artifacts after resumption. +# +# Some 100K USGS maps come with multiple sets of boundary files, +# leading to file name collisions for which no workaround has been +# found. +# +############################################################################ +# +# Requires: Version 9 graphics +# +############################################################################ +# +# Links: cartog, clipping, ddfread, geodat, graphics, io, mapnav, +# numbers, options, pscript, random, strings, wildcard, zipread +# +############################################################################ + + + +$include "keysyms.icn" + +link cartog +link clipping +link ddfread +link geodat +link graphics +link io +link mapnav +link numbers +link options +link pscript +link random +link strings +link wildcard +link zipread + + + +$define DLG_LAYERS "boclvdwsrtum" # all "real" layers, in loading order + +$define WSIZE "size=1000,1000" # default window size +$define ZOOMF 1.5 # zoom factor + +$define MAXFONT 18 # maximum font size +$define MINFONT 8 # minimum font size +$define MINBOLD 10 # minimum bold font size + +$define MEGABYTE (1024 * 1024) # how many bytes in a megabyte? +$define STRSIZE (1 * MEGABYTE) # default string region size +$define BLKSIZE (16 * MEGABYTE) # default block region size +$define MAXDRAW 4000 # maximum (even) args to avoid error 301 + +$define DEGSCALE 1.0e+6 # divisions per degree for integer scale + +# parameters for displaying progress during loading +$define PINTERVAL 100 # progress interval +$define PSQUARES 8 # number of progress squares +$define PSQSIZE 10 # size of progress squares +$define PSQGAP 2 # size of gap between squares + +# PostScript output parameters +$define PSSCALE 10 # scaling from pixels to PS units +$define PSPT 24 # size of point feature in PS units +$define PSLWI 120 # linewidth scaling factor + + + +record arg ( # command-line argument, excluding options + name, # file or directory name + type, # "dir", "zip", or "txt" + wanted, # if a wanted layer (null if suppressed by -t option) + ltable, # table of layer records, indexed by layer code + pcount # progress bar counter + ) + +record layer ( # one layer in one directory (or zip file) + lcode, # layer code character + arg, # corresponding arg record + files, # list of file names + zone, # UTM zone, or -1 if data is in lat/lon , from XREF file + xscale, yscale, # scaling factors for file values + datum, # stated coordinate datum, from XREF file + ellipsoid, # inferred geodetic ellipsoid + icorners, # map corners in input terms + ocorners, # map corners as projected to lat,lon + px, py, # progress reporting coordinates + wd # width of layer in screen units + ) + +record attrec ( # line drawing attributes: + seq, # drawing sequence + lcode, # layer code + key, # table key (layer or attribute code) + width, # line width + color, # line color + style, # line style + segs # list of segments (list of paths) + ) + +record feature ( # feature or waypoint + lat, # latitude + lon, # longitude + label # label + ) + + + +global arglist # list of arg records +global opts # command options +global chosen # cset of chosen layers + +global xmin, xmax, ymin, ymax # data range +global aspect # input coordinate aspect ratio + +global attrib # attribute style table +global slist # list of style records w/ seg lists +global pcolors # list of path background colors + +global features # list of feature records + +global logfile # feature log file, if any + + + +# main program + +procedure main(args) + local a, c, e, g, i, r, s, t, v + + # use large region sizes for better efficiency + collect(2, STRSIZE) # string region + collect(3, BLKSIZE) # block (heap) region + + # open window first, to validate and remove any window options + Window("label=--", "gamma=1.5", "bg=white", "fg=brown", + "resize=on", "canvas=hidden", WSIZE, args) + + randomize() + initattrib() + + # process command options + opts := options(args, "o:l:x:cdnpqt") + if \opts["o"] then { + if opts["o"] == "-" then + logfile := &output + else + logfile := open(opts["o"], "w") | stop("cannot write ", opts["o"]) + } + else + logfile := &output + + chosen := cset(\opts["l"]) | (&lcase -- 'n') # start with explicit layers + chosen ++:= 'go' # add paths & other data, if loaded + chosen --:= cset(\opts["x"]) # now apply exclusions + + # any remaining arguments are directory names + if *args = 0 then + stop("usage: ", &progname, " [options] dir...") + + # build list of arg records, classifying each filename or directory + arglist := [] + every s := !args do { + if directory(s) then + t := "dir" + else if iszip(s) then + t := "zip" + else + t := "txt" + put(arglist, arg(s, t, 1)) + } + + # scan text files first, because we haven't really done any validation + # (any unrecognized file is classified as a text file) + features := [] + every (a := !arglist) & (a.type == "txt") do + rdtext(a) + + # take inventory of DLG directories and files, and load XREF/NPnn info + every (a := !arglist) & (a.type ~== "txt") do { + inventory(a) + every r := !a.ltable do { + loadref(r) + if r.zone >= 0 then + loadcorners(r) + else + loadbounds(r) + } + if \opts["t"] & not traversed(!a.ltable) then + a.wanted := &null + lstats(a) + } + + if \opts["n"] then + return + (*(!arglist).ltable > 0) | stop("no data") + + # show initial screen + winit() + mapinit(draw, , xmin, xmax, ymax, ymin, aspect) + if WAttrib("label") == "--" then # set window label, if not specified + WAttrib("label=" || args[1]) + WAttrib("canvas=normal") # make window visible + Font("sans,bold,72") + Fg("pale yellowish gray") + DrawString(60, 120, "LOADING...") + + if \opts["c"] then # if just coverage wanted + chosen := 'n' # turn on names, turn off loaded paths + else { + + # finally: load in the data + alllabels() # show coverage while loading + every c := !DLG_LAYERS do # load by layers + every a := !arglist do + if \a.wanted then + loadlayer(\a.ltable[c]) + + # report memory usage + every put(c := [], &collections) + collect() + every put(a := [], &storage) + if /opts["q"] then { + write(" ", (a[2] + a[3] + MEGABYTE / 2) / MEGABYTE, + " MB loaded (", c[3], "+", c[4], " GC)") + } + } + + # put segment lists in order for drawing + # shuffle segments of each list to minimize "dead time" drawing offscreen + every put(slist := [], !attrib) + slist := sortf(slist, field(attrec, "seq")) + every g := (!slist).segs do + every !g :=: ?g # imperfect but good enough shuffle + + # report attribute counts, if -d given + if \opts["d"] then { + write() + every e := !slist do + if *e.segs > 0 then + write(right(e.seq, 3), ". ", e.lcode, " ", + left(e.key, 8), right(*e.segs, 7)) + write() + } + + # consume any events that may have occurred during loading + while *Pending() > 0 do + Event() + + # draw initial screen + EraseArea() + mapgen() + + # process interactive commands + repeat case e := Event() of { + &shift & &lpress: { logfeat(e) } + &shift & &rpress: { logfeat(e) } + &rpress: { locate() } + !"\n\r ": { mapgen() } + !"pP" | Key_PrSc: { print(); Bg("white") } + !"aA": { chosen := &lcase; mapgen() } + !"xX": { chosen := ''; EraseArea(); mapgen() } + !"qQ": { exit() } + any(&letters, e) & e: { + e := map(e) + if any(chosen, e) then { + chosen --:= e + EraseArea() + mapgen() + } + else { + chosen ++:= e + mapgen() + } + } + default: { mapevent(e) } + } +end + + + +# rdtext(arg) -- read a text file of paths and features + +procedure rdtext(arg) + local f, i, n, r, s, t, w, line + local lat, lon, alt, segs, points, nsegs, npts, nfeat + local xmn, xmx, ymn, ymx + static npaths + initial npaths := 0 + + f := open(arg.name) | stop("cannot open: ", arg.name) + s := "g" || (npaths % *pcolors + 1) + npaths +:= 1 + + segs := attrib[s].segs + nsegs := *segs + npts := 0 + nfeat := 0 + xmn := ymn := +180 * DEGSCALE + xmx := ymx := -180 * DEGSCALE + + points := [] + while line := read(f) do { # read line + every put(w := [], words(line)) # break into fields + # check first for path entry + if (lat:=real(w[-3])) & (lon:=real(w[-2])) & (alt:=real(w[-1])) & + (-90.<=lat<=90.) & (-180.<=lon<=180.) & (-1400<alt<30000) then { + npts +:= 1 + lon *:= DEGSCALE + lat *:= DEGSCALE + put(points, integer(lon), integer(lat)) + xmn >:= lon + ymn >:= lat + xmx <:= lon + ymx <:= lat + } + else if (lat := real(w[-2])) & (lon := real(w[-1])) & + (-90. <= lat <= 90.) & (-180. <= lon <= 180.) then { + npts +:= 1 + lon *:= DEGSCALE + lat *:= DEGSCALE + put(points, integer(lon), integer(lat)) + xmn >:= lon + ymn >:= lat + xmx <:= lon + ymx <:= lat + } + else { + # interrupt path sequence + if *points > 0 then { + put(segs, points) + points := [] + } + # check for feature (waypoint) label + if (lat := real(get(w))) & (lon := real(get(w))) & + (-90. <= lat <= 90.) & (-180. <= lon <= 180.) then { + nfeat +:= 1 + lon *:= DEGSCALE + lat *:= DEGSCALE + xmn >:= lon + ymn >:= lat + xmx <:= lon + ymx <:= lat + s := "" + while s ||:= " " || get(w) + put(features, feature(lat, lon, s[2:0])) + } + } + } + if *points > 0 then + put(segs, points) + + nsegs := *segs - nsegs + if nsegs = 0 & nfeat = 0 then + stop("no data: ", arg.name) + + r := layer("g", arg) + r.zone := -1 + r.datum := "WGS84" + r.ellipsoid := "WGS84" + r.icorners := r.ocorners := [xmn, ymn, xmn, ymx, xmx, ymx, xmx, ymn] + t := table() + t["g"] := r + arg.ltable := t + + n := 0 + every n +:= *segs[-nsegs to -1] + if /opts["q"] then + write(right(nsegs || ":" || npts || "+" || nfeat, 14), " ", lsumm(r)) + + close(f) + return +end + + + +# ddpopen(r, p) -- generate open DDF files from layer r matching pattern p + +procedure ddpopen(r, p) + local a, f, d, s, fname + + a := r.arg + every fname := !r.files do { + if not (map(fname) ? wild_match(p)) then + next + s := a.name || "/" || fname + f := &null + if a.type == "zip" then + f := zipfile(a.name, fname) + else + f := open(s, "ru") + d := ddfopen(\f) | stop("cannot open as DDF: ", s) + suspend d + } + fail +end + + + +# inventory(a) -- inventory arg entry a + +procedure inventory(a) + local b, c, f, fname, m, flist, trcount + + # load filenames into list, because we need to scan it twice + flist := [] + if a.type == "zip" then + f := zipdir(a.name) + else + f := open(a.name) + while put(flist, read(f)) + close(f) + + # count TR01LE??.DDF files + trcount := 0 + every fname := !flist do + if map(fname) ? (tab(-12) & ="tr01le") then + trcount +:= 1 + + # classify files and save the ones we want + a.ltable := table() + every fname := !flist do { + map(fname) ? { + while tab(upto('/') + 1) + pos(-12) | next + move(8) | next + =".ddf" | next + } + b := fname[-12:-4] | next + + every c := !lcodes(b, trcount) do { + if any(chosen, c) then { + # this is a wanted file in a wanted layer; remember it + /a.ltable[c] := layer(c, a, []) + put(a.ltable[c].files, fname) + } + } + } + + return +end + + + +# lcodes(basename, trcount) -- deduce layer code(s) from file basename + +procedure lcodes(basename, trcount) + local n, s, tr + + map(basename) ? { + if move(4) & ="a" & move(2) & any('f') then { + # xxxxAllF.DDF is layer ll attribute file + s := move(-2) + } + else if ="tr01" & =("le" | "np") & (n := integer(move(2))) then { + # TR01LEnn.DDF (or NPnn) is a transportation layer in a 100K map + if trcount > 12 then + s := ["mt", "rd", "rd", "rd", "rd", "rr"] [(n + 3) / 4] + else + s := ["mt", "rd", "rr"] [(n + 3) / 4] + } + else if move(2) & ="tr" & =("le" | "ne") & (n := integer(move(2))) then { + # xxTRLEnn.DDF (or NExx) is a transportation layer in state xx 250K map + s := ["mt", "rd", "rr"] [n % 3 + 1] + } + else { + move(2) + if any(&letters) then + s := move(2) # xxllyyyy is layer ll for state xx + else + s := move(-2) # ll01xxxx is layer ll otherwise + } + } + + return case s of { + "bd": "b" # boundaries (BD: boundaries) + "hp": "c" # contours (HP: hypsography) + "nv": "d" # sand etc. (NV: nonvegetative features) + "pl": "l" # land sections (PL: public lands) + "sm": "m" # markers (SM: survey markers) + "rd": "r" # roads (RD: roads) + "ms": "s" # structures (MS: manmade structures) + "rr": "t" # train tracks (RR: railroads) + "mt": "u" # utilities (MT: miscellaneous transportation) + "tr": "rtu" # transportatn (TR: transportation, shared by r/t/u) + "sc": "v" # vegetation (SC: surface cover) + "hy": "w" # water (HY: hydrology) + default: "o" # other + } +end + + + +# getdata(r, p, l) -- get data vector l of layer r using file pattern p + +procedure getdata(r, p, l) + local ddfile, d, e, zone + + ddfile := ddpopen(r, p) | + stop("no file ", p, " for layer ", r.lcode, ": ", r.arg.name) + while d := ddfread(ddfile) do + every e := !d do + if e[1] == l then + break break + ddfclose(ddfile) + return e +end + + + +# loadref(r) -- load XREF and IREF files for layer r of arg a + +procedure loadref(r) + local e + + e := getdata(r, "*iref.ddf", "IREF") + until get(e) == "BI32" + r.xscale := real(get(e)) + r.yscale := real(get(e)) + + e := getdata(r, "*xref.ddf", "XREF") + case e[5] of { + "NAS": r.datum := "NAD27" # North American 1927 + "NAX": r.datum := "NAD83" # North American 1983 + "WGA": r.datum := "WGS60" # World Geodetic System 1960 + "WGB": r.datum := "WGS66" # World Geodetic System 1966 + "WGC": r.datum := "WGS72" # World Geodetic System 1972 + "WGE": r.datum := "WGS84" # World Geodetic System 1984 + default: r.datum := "?????" # unrecognized + } + if e[4] == "UTM" then + r.zone := integer(e[6]) + else + r.zone := -1 + + return +end + + + +# loadbounds(r) -- load SPDM file to determine range of locations +# +# (SPDM files are used with 250K DLG layers) + +procedure loadbounds(r) + local e, xmn, xmx, ymn, ymx + + e := getdata(r, "*spdm.ddf", "DMSA") + get(e) + xmn := get(e) * r.xscale * DEGSCALE + ymn := get(e) * r.yscale * DEGSCALE + xmx := get(e) * r.xscale * DEGSCALE + ymx := get(e) * r.yscale * DEGSCALE + r.ellipsoid := "Clarke66" + r.icorners := r.ocorners := [xmn, ymn, xmn, ymx, xmx, ymx, xmx, ymn] + return +end + + + +# loadcorners(r) -- load NPnn file to determine corner points +# +# (NPnn files are used with 24K and 100K DLG layers) + +procedure loadcorners(r) + local ddfile, d, e, i, x, y, L, C66, G80 + + every ddfile := ddpopen(r, "*np??.ddf") do { + L := [] + while d := ddfread(ddfile) do + every e := !d do + if get(e) == "SADR" then + while put(L, get(e)) + ddfclose(ddfile) + r.icorners := cmerge(r.icorners, L) + } + + if /r.icorners then + stop("no NPnn file for layer ", r.lcode, ": ", r.arg.name) + + # infer ellipsoid of UTM projection + L := [] + every i := 1 to *r.icorners by 2 do { + x := (r.icorners[i] * r.xscale - 500000.0) + y := (r.icorners[i+1] * r.yscale) + put(L, r.zone, x, y) + } + C66 := project(invp(utm("Clarke66")), L) + G80 := project(invp(utm("GRS80")), L) + + if quadfit(C66) < quadfit(G80) then { + r.ellipsoid := "Clarke66" + r.ocorners := project(molodensky("NAD27", "NAD83"), C66) + } + else { + r.ellipsoid := "GRS80" + r.ocorners := G80 + } + + every !r.ocorners *:= DEGSCALE + return +end + + + +# cmerge(A, B) -- merge two corners lists +# +# Assumes that the corner order is [SW, NW, NE, SE] +# and takes the more extreme value for each coordinate. + +procedure cmerge(A, B) + local C + + if /A | /B then return \A | \B + C := [] + + if A[1] + A[2] < B[1] + B[2] then + put(C, A[1], A[2]) + else + put(C, B[1], B[2]) + + if A[3] - A[4] < B[3] - B[4] then + put(C, A[3], A[4]) + else + put(C, B[3], B[4]) + + if A[5] + A[6] > B[5] + B[6] then + put(C, A[5], A[6]) + else + put(C, B[5], B[6]) + + if A[7] - A[8] > B[7] - B[8] then + put(C, A[7], A[8]) + else + put(C, B[7], B[8]) + + return C +end + + + +# quadfit(L) -- proximity of coordinate in L to multiple of 1/8 + +procedure quadfit(L) + local i, mn, mx, a, b + + mn := 1.0 + every i := 1 to *L by 2 do { + a := L[i] * 8 + b := L[i+1] * 8 + mx := max(abs(a - round(a)), abs(b - round(b))) + mn := min(mn, mx) + } + return mn +end + + + +# lstats(a) -- report statistics for the layers of arg a + +procedure lstats(a) + local c, d, g, k, l, n, r, v, z + + if \opts["q"] then + return + + # group by identical projection attributes + g := table('') + every r := !a.ltable do { + k := lsumm(r) + g[k] ++:= r.lcode + } + + # report consistent layers together on one line + l := sort(g, 3) + while k := get(l) do { + v := get(l) + writes(if /a.wanted then "X" else " ") + writes(if *g = 1 then " " else "*") + every c := !cset(DLG_LAYERS) do # list alphabetically + writes(if upto(v, c) then c else "-") + write(" ", k) + } + return +end + + + +# lsumm(r) -- return one-line layer info summary + +procedure lsumm(r) + return r.datum[1] || r.datum[-2:0] || " " || + (if r.zone < 0 then "LL" else right(r.zone, 2)) || " " || + r.ellipsoid[1] || r.ellipsoid[-2:0] || " " || + right(degc(r.ocorners[-1]), 3) || " " || + left(degc(r.ocorners[-2]), 4) || " " || + r.arg.name +end + + + +# degc(d) -- code degree measurement as nnnx where x is a-h for 0/8 to 7/8 + +procedure degc(d) + local n, x + + d := abs(d / DEGSCALE) + 0.0625 # 1/16 for rounding + n := integer(d) + x := "abcdefgh" [1 + integer(8 * (d - n))] + return n || x +end + + + +# field(constr, key) -- given record constructor, find index of named field + +procedure field(constr, key) + local i, r + image(constr) ? ="record constructor" | fail + r := constr() + every i := 1 to *r do + r[i] := i + return r[key] +end + + + +# traversed(r) -- check whether layer r is traversed by a path + +procedure traversed(r) + local k, i, segs, pts, xmin, xmax, ymin, ymax + + xmin := xmax := r.ocorners[1] + every xmin >:= r.ocorners[3 | 5 | 7] + every xmax <:= r.ocorners[3 | 5 | 7] + ymin := ymax := r.ocorners[2] + every ymin >:= r.ocorners[4 | 6 | 8] + every ymax <:= r.ocorners[4 | 6 | 8] + + every k := key(attrib) do + if k ? (="g" & tab(many(&digits)) & pos(0)) then + every pts := !attrib[k].segs do + every i := 1 to *pts by 2 do + if xmin < pts[i] < xmax & ymin < pts[i+1] < ymax then + return + + fail +end + + + +# loadlayer(r) -- load one layer of files + +procedure loadlayer(r) + local p, attid, ddfile + + setdraw(attrib[r.lcode]) + drawlabel(r) + + attid := table() + every ddfile := ddpopen(r, "*a??f.ddf") do { + loadatts(ddfile, r, attid) + ddfclose(ddfile) + } + + every ddfile := ddpopen(r, "*ne??.ddf" | "*le??.ddf") do { + loadpts(ddfile, r, attid) + ddfclose(ddfile) + } + + return +end + + + +# loadatts(ddfile, r, attid) -- load attribute ID table + +procedure loadatts(ddfile, r, attid) + local d, e, i, k, n, s, v + + n := -1 + if r.lcode == "t" then + i := [1, 7] # for RR, append tunnel and rapid transit flags + else + i := [] + + while d := ddfread(ddfile) do { + k := &null + every e := !d do { + s := get(e) + if s == "ATPR" then + k := get(e) || get(e) + else if s == "ATTP" then { + v := get(e) + every \v ||:= (" " ~== e[!i]) + attid[\k] := v + if (n +:= 1) % PINTERVAL = 0 then + progress(r) + } + } + } + return +end + + + +# loadpts(ddfile, r, attid) -- load coordinate file into memory + +procedure loadpts(ddfile, r, attid) + local a, d, e, i, k, m, n, p, s, v, vv, x, y + local lcode, zone, coords, arec + + lcode := r.lcode + zone := r.zone + + if zone >= 0 then { # if not already in lat/lon form + if /opts["p"] then { # if no -p option + p := pptrans(r.icorners, r.ocorners) # use approx, faster projection + zone := &null # indicate such for code below + } + else { + p := invp(utm(r.ellipsoid)) # use full inverse-UTM projection + if r.ellipsoid == "Clarke66" then # and if needed, + p := compose(molodensky("NAD27", "NAD83"), p) # datum conversion + } + } + + n := 0 + while d := ddfread(ddfile) do { + a := lcode || "-" + v := [] + coords := [] + every e := !d do { + if *e < 3 then + next + s := get(e) + if s == "ATID" then { + k := get(e) || get(e) + while k[4] ~== "F" do + k := get(e) || get(e) | break + a := \attid[k] | lcode + } + else if s == "SADR" then { + if /p then { + # latitude/longitude input + while x := get(e) & y := get(e) do + put(v, x * r.xscale * DEGSCALE, y * r.yscale * DEGSCALE) + } + else if /zone then { + # using approximate projection, which includes scaling + while x := get(e) & y := get(e) do + put(coords, x, y) + } + else { + # full inverse UTM projection + while x := get(e) & y := get(e) do + put(coords, zone, x * r.xscale - 500000.0, y * r.yscale) + } + } + } + + if \p then { # if projection needed + coords := project(p, coords) # project UTM to lat/lon + m := if /zone then 1 else DEGSCALE # select multiplier + while put(v, integer(m * get(coords))) # convert to scaled integer + } + + if *v = 0 then + next + + if not (arec := \attrib[a]) then { + # add unrecognized attribute code to table + arec := copy(attrib[lcode]) + arec.key := a || "*" # "*" indicates unregistered attribute + arec.segs := [] + attrib[a] := arec + } + + while *v > MAXDRAW do { # break too-large path into pieces + vv := [] + every 3 to MAXDRAW by 2 do + put(vv, get(v), get(v)) # move out of v + put(vv, v[1], v[2]) # leave one point for overlap + put(arec.segs, vv) # store extracted piece + } + + # loops are rare in the data, but can crash XFree86 server if dashed + if v[1] = v[-2] & v[2] = v[-1] then { # if loop + put(v, v[3], v[4]) # overshoot to 2nd point again + } + put(arec.segs, v) # store what's left of original + + if (n +:= 1) % PINTERVAL = 0 then + progress(r) + + } + return +end + + + +# logfeat() -- record current location to log file + +procedure logfeat(e) + local ll, lat, lon, locn, label + + until Event() === (&lrelease | &rrelease) # wait for button up + ll := project(invp(mapproj()), [&x + 0.5, &y + 0.5]) # cvt coords to lat/lon + lon := get(ll) / DEGSCALE + lat := get(ll) / DEGSCALE + locn := frn(lat, 0, 6) || " " || frn(lon, 0, 6) + label := "" + + if e === &lpress then { # if named (not anonymous), ask for label + setdraw(attrib["DIALOG"]) + VSetFont() + case TextDialog( + ["Enter label for", locn || ":"], , , 30, ["Okay", "Cancel"]) of { + "Okay": label := " " || get(dialog_value) + "Cancel": fail + } + put(features, feature(DEGSCALE * lat, DEGSCALE * lon, label[2:0])) + if any(chosen, "f") then + allfeats(mapproj(), Pending) # redraw feats to display label + } + + write(logfile, locn, label) + flush(logfile) + return +end + + + +# locate() -- display location while right button is held down + +$define BOXW 265 # popup box width +$define BOXH 90 # popup box height +$define SMAX (BOXW - 40) # maximum scalebar length + +procedure locate() + + setdraw(attrib["DIALOG"]) # set colors and font for drawing + Font("mono,bold,14") + + if &x < BOXW + 40 & &y < BOXH + 40 then + Popup(20, WAttrib("height") - BOXH - 20, BOXW, BOXH, locproc, mapproj()) + else + Popup(20, 20, BOXW, BOXH, locproc, mapproj()) + + return +end + + + +# locate(wproj) -- calculate scale and location using caller's projection + +procedure locproc(wproj) + local d, e, m, s, u, cx, dx, dy, ll, lat, lon, winv + + winv := invp(wproj) # get projection from screen to lat/lon + dx := WAttrib("dx") # get popup box coordinate system + dy := WAttrib("dy") + + # compute a reasonably round length that works for a scale bar + u := 90 * DEGSCALE / 1e7 # one meter, in latitude units + m := sbsize(wproj, xmin, ymin, u, SMAX) + + # draw the scale bar + ll := project(wproj, [xmin, ymin, xmin + m * u, ymin]) + d := ll[3] - ll[1] + cx := BOXW / 2 + FillRectangle(cx - d / 2, 55, d, 8) + + if m >= 1000 then + s := (m / 1000) || " km" + else + s := m || " m" + CenterString(cx, 70, s) + + # give coordinates of mouse location until button released + until e === &rrelease do { + + ll := project(winv, [&x + 0.5, &y + 0.5]) # cvt screen coords to lat/lon + lon := get(ll) / DEGSCALE # and scale from integer to real + lat := get(ll) / DEGSCALE + + GotoRC(1, 1) + WWrites("\n ", dms(lat, "S", "N"), frn(lat, 13, 6)) + WWrites("\n ", dms(lon, "W", "E"), frn(lon, 13, 6)) + + e := Event() # get next event + &x +:= dx # remove effect of popup box coordinate system + &y +:= dy + } + + return +end + +procedure dms(n, s1, s2) + local deg, min, sec + + if n < 0 then + n := -n + else + s1 := s2 + + n +:= 1 / 7200. # rounding + deg := integer(n); n := (n - deg) * 60 + min := integer(n); n := (n - min) * 60 + sec := integer(n) + + return s1 || right(deg, 4) || "\260" || right(min, 2, "0") || "'" || + right(sec, 2, "0") || "\"" +end + + + +# draw(win, pjn) -- draw all selected map layers, without erasing first + +procedure draw(win, pjn) + local a, d, v, arec + + every (arec := !slist) & any(chosen, arec.lcode) do { + setdraw(arec) | next + every d := !arec.segs do { + v := project(pjn, d) # project to window x/y coords + every !v >:= 30000.0 # clamp to legal X values allowing dx/dy + every !v <:= -30000.0 # clamp as floating to avoid lgint bug + if *v = 2 then + FillRectangle(v[1] - 1, v[2] - 1, 3, 3) + else + DrawLine ! v + if *Pending() > 0 then + return + } + } + + # draw feature (waypoint) labels + if any(chosen, "f") then + allfeats(pjn, Pending) + + # draw pseudo-layer "n" + if any(chosen, "n") then + alllabels(Pending) + + collect() # do this now, while awaiting input + return +end + + + +# winit() -- initialize window configuration + +procedure winit() + local a + + xmin := ymin := +180 * DEGSCALE + xmax := ymax := -180 * DEGSCALE + every a := !arglist do + if \a.wanted then { + every xmin >:= (!a.ltable).ocorners[1 | 3] + every xmax <:= (!a.ltable).ocorners[5 | 7] + every ymin >:= (!a.ltable).ocorners[2 | 8] + every ymax <:= (!a.ltable).ocorners[4 | 6] + } + aspect := cos(dtor((ymax + ymin) / (2 * DEGSCALE))) + return +end + + + +# allfeats(pjn, p) -- draw feature labels +# +# p is Pending procedure, if to check and quit early + +procedure allfeats(pjn, p) + local f, x, y, xy, xy2 + + xy := [] + every f := !features do + put(xy, f.lon, f.lat) + xy := project(pjn, xy) + xy2 := copy(xy) + + Font("sans, bold, 10") + setdraw(attrib["f"]) + Fg("white") # draw offset backgrounds in white + every f := !features do { + DrawString(get(xy2) + 4, get(xy2) + 5, f.label) + if *(\p)() > 0 then + break + } + + setdraw(attrib["f"]) # draw labels in black + every f := !features do { + x := get(xy) + y := get(xy) + FillRectangle(x - 1, y - 1, 3, 3) + DrawString(x + 5, y + 4, f.label) + if *(\p)() > 0 then + break + } + + return +end + + + +# alllabels(p) -- draw labels for all layers in standard color +# +# p is Pending procedure, if to check and quit early + +procedure alllabels(p) + local a, r + + setdraw(attrib["n"]) + every a := !arglist do { + if \a.wanted then { + drawlabel(!a.ltable) # pick any layer + if \opts["c"] then { + drawcoverage(a) + setdraw(attrib["n"]) + } + } + if *(\p)() > 0 then + break + } + return +end + + + +# drawlabel(r) -- draw label for layer r in current color +# +# sets r.px, r.py to progress bar position and r.wd to layer width + +procedure drawlabel(r) + local x, y, w, h, n, d, u, s, tw, tmax, v, wproj + static lc, uc + initial { + lc := string(&lcase) + uc := string(&ucase) + } + + # draw the bounding box + wproj := mapproj() + v := copy(r.ocorners) + put(v, r.ocorners[1], r.ocorners[2]) + v := project(wproj, v) # project to window x/y coords + every !v >:= 30000.0 # clamp to legal X values allowing dx/dy + every !v <:= -30000.0 # clamp as floating to avoid lgint bug + DrawLine ! v + + # find the center and range + x := (v[1] + v[3] + v[5] + v[7]) / 4 + y := (v[2] + v[4] + v[6] + v[8]) / 4 + w := (v[5] + v[7] - v[1] - v[3]) / 2 + h := (v[4] + v[6] - v[2] - v[8]) / 2 + + # trim the name + s := r.arg.name + while s[-1] == "/" do + s := s[1:-1] + s ? { + while tab(upto('/') + 1) + s := map(tab(0), lc, uc) + } + if s[-4:0] == (".ZIP" | ".GPS" | ".RTE" | ".TRK") then + s := s[1:-4] + + # draw the label + Font("sans,bold," || MAXFONT) + tw := TextWidth(s) + tmax := .90 * w + if tw > tmax then { + n := integer(MAXFONT * tmax / tw) + if n <:= MINFONT then { + # it doesn't fit, and will overlap neighbors with minimum font size; + # add pseudorandom vertical offset to mitigate overlap + d := abs(r.ocorners[7] / DEGSCALE) # SE corner longitude + u := integer(8 * d + 0.5) # 1/8-degree units + u +:= integer(2 * d + 0.5) # half-degree units + y -:= 0.20 * h * (1.5 - u % 4) + } + if n < MINBOLD then + Font("sans," || n) + else + Font("sans,bold," || n) + } + CenterString(x, y, s) + + r.px := integer(x) + r.py := integer(y + 0.75 * WAttrib("fheight")) + r.wd := w + return +end + + + +# progress(r) -- draw progress square for layer r + +procedure progress(r) + local a, x + + a := r.arg + a.pcount := (\a.pcount + 1) | 0 + x := r.px + PSQSIZE * (a.pcount % PSQUARES - PSQUARES / 2) + FillRectangle(x, r.py, PSQSIZE - PSQGAP, PSQSIZE - PSQGAP) + if (a.pcount / PSQUARES) % 2 = 1 then + EraseArea(x + 1, r.py + 1, PSQSIZE - PSQGAP - 2, PSQSIZE - PSQGAP - 2) + return +end + + + +# drawcoverage(a) -- draw coverage indicators for arg entry a + +procedure drawcoverage(a) + local c, r, x, y, w + + r := \!a.ltable | return + w := r.wd / *DLG_LAYERS + w >:= PSQSIZE + w <:= 2 + x := r.px - (w * *DLG_LAYERS) / 2 + y := r.py + + every c := !cset(DLG_LAYERS) do { + if r := \a.ltable[c] then { + setdraw(attrib[r.lcode]) + FillRectangle(x, y, w, w) + } + x +:= w + } + return +end + + + +# print() -- print visible portion to file + +procedure print() + local psname, psfile + + Bg("pale weak brown") + VSetFont() + setdraw(attrib["DIALOG"]) # set reasonable colors for dialog + repeat case OpenDialog("Print to file:") of { + "Okay": { + if *dialog_value = 0 then + next + if close(open(psname := dialog_value)) then + case TextDialog("Overwrite existing file?", , , , + ["Yes", "No", "Cancel"]) of { + "Yes": &null + "No": next + "Cancel": fail + } + if psfile := open(psname, "w") then + break + case TextDialog("Cannot write " || psname) of { + "Okay": next + "Cancel": fail + } + } + "Cancel": + fail + } + + Popup(, , 300, 50, + popwrite, [psfile, mapproj(), WAttrib("width"), WAttrib("height")]) + close(psfile) + return +end + +procedure popwrite(psargs) + CenterString(150, 25, "Writing PostScript...") + return writeps ! psargs +end + +procedure writeps(psfile, projn, wwidth, wheight) + local arec, color, style, width, ptoff, xmax, ymax, xmul, ymul + local a, b, f, m, w, h, pj, d, s, u, v, x, y, dx, dy, fx, fy, ll + + b := project(invp(projn), [0, 0, wwidth, wheight]) + xmax := PSSCALE * wwidth + ymax := PSSCALE * wheight + xmul := xmax / (b[3] - b[1]) + ymul := ymax / (b[2] - b[4]) + pj := rectp(b[1], b[4], 0, 0, xmul, ymul) # set projection + + ptoff := PSPT / 2 + s := " 0 " || PSPT || " rlineto" + s ||:= " " || PSPT || " 0 rlineto" + s ||:= " 0 -" || PSPT || " rlineto" + + epsheader(psfile, 0, 0, PSSCALE * wwidth, PSSCALE * wheight, "r") + every write(psfile, ![ + "1 setlinecap", + "/cdivr { 65535 div 3 1 roll } bind def", + "/color { cdivr cdivr cdivr setrgbcolor } bind def", + "/solid { [] 0 setdash } bind def", + "/dashed { [ .04 inch dup ] 0 setdash } bind def", + "/m { moveto } bind def", + "/r { rlineto } bind def", + "/s { rlineto stroke } bind def", + "/p { moveto" || s || " fill } bind def", + "/f { 2 copy p moveto 48 -36 rmoveto show } bind def", + ]) + + every (arec := !slist) & any(chosen, arec.lcode) do { + if *arec.segs = 0 | arec.width < 0 then + next + if color ~===:= arec.color then + write(psfile, map(ColorValue(arec.color), ",", " "), " color") + if width ~===:= arec.width then + write(psfile, arec.width / real(PSLWI), " inch setlinewidth") + if style ~===:= arec.style then + write(psfile, style) + every d := !arec.segs do { + v := project(pj, d) + if *v = 2 then { + x := integer(get(v)) + y := integer(get(v)) + if (0 <= x < xmax) & (0 <= y < ymax) then + write(psfile, x - ptoff, " ", y - ptoff, " p") + next + } + v := Coalesce(ClipLine(v, 0, 0, xmax, ymax)) | next + every a := !v do { + x := integer(get(a)) + y := integer(get(a)) + fy := integer(pull(a)) + fx := integer(pull(a)) + write(psfile, x, " ", y, " m") + while dx := integer(get(a) - x) do { + dy := integer(get(a) - y) + write(psfile, dx, " ", dy, " r") + x +:= dx + y +:= dy + } + write(psfile, fx - x, " ", fy - y, " s") + } + } + } + + # write features + if *features > 0 & any(chosen, "f") then { + write(psfile) + write(psfile, "/Times-Roman findfont 120 scalefont setfont") + write(psfile, "0 0 0 color") + every f := !features do { + a := project(pj, [f.lon, f.lat]) + x := integer(get(a)) + y := integer(get(a)) + if (0 <= x <= xmax) & (0 <= y <= ymax) then + write(psfile, "(", psprotect(f.label), ") ", + x - ptoff, " ", y - ptoff, " f") + } + } + + # write scale bar + u := 90 * DEGSCALE / 1e7 # one meter, in latitude units + m := sbsize(pj, xmin, ymin, u, 2000) + ll := project(pj, [xmin, ymin, xmin + m * u, ymin]) + d := ll[3] - ll[1] + if m >= 1000 then + s := (m / 1000) || " km" + else + s := m || " m" + + every write(psfile, ![ + "", + "0 0 0 color", + "0 0 m 0 120 r " || d || " 0 r 0 -120 r fill", + "/Helvetica findfont 100 scalefont setfont", + "65535 65535 65535 color", + integer(d / 2 - 120) || " 25 m (" || s || ") show", + ]) + + write(psfile, "showpage") + return +end + + + + +# initattrib() -- initialize drawing attributes +# +# IMPORTANT: Map entities are drawn in the order of the def() calls below. + +procedure initattrib() + local i, s + +$define ROUTE "magenta-red" # path foreground color + pcolors := [ # path background colors + "yellow", # yellow + "light green", # green + "light bluish cyan", # blue + "reddish yellow", # orange + "pale purple", # purple + "pale red-yellow", # peach + "pale moderate green", # greenish gray + "pale moderate cyan", # bluish gray + ] + pull(pcolors) # remove trailing null + + attrib := table() + deflayer(" ", "black") + def("SWEEP", 3, "reddish orange") # interactive sweeping with mouse + def("DIALOG", 1, "black") # dialog boxes + + every i := 1 to *pcolors do { + s := "g" || i + deflayer(s, ROUTE) # paths (first drawing) + def(s || "b", 10, pcolors[i]) # faint, wide highlighter background + def(s || "f", 2, ROUTE) # bold foreground + } + + deflayer("b", "light reddish yellow") # boundaries (wide, so draw first) + def("b", 3) + + deflayer("o", Fg()) # unknown other data; use specified Fg + def("o") + + deflayer("c", "light red-yellow") # contour lines (hypsography) + def("c-", , "pale moderate red-yellow") # deemphasize unattributed segments + def("c") # contour line + def("0200205", , "light moderate bluish-cyan") # bathymetric contour + def("0200206", , "light moderate bluish-cyan") # depth curve + def("0200210", , "light moderate bluish-cyan") # suppl bathymetric contour + def("0200207", , "deep red-yellow") # watershed (e.g. continental) divide + + deflayer("l", "pale whitish red") # land sections + def("l") + + deflayer("v", "light green") # vegetation (surface cover) + def("v") # surface cover + + deflayer("d", "light weak green") # gravel etc. (nonvegetative features) + def("d") + + deflayer("w", "bluish cyan") # water (hydrology) + def("w-", , "pale bluish cyan") # deemphasize unattributed segments + def("0500415", , , "dashed") # aqueduct or water pipeline + def("0500200", 2) # shoreline + def("0500201", 2) # manmade shoreline + def("w") # unspecified hydrology + def("0500412") # stream + + deflayer("s", "weak reddish yellow") # manmade structures + def("2000299", , "pale reddish yellow") # processing line + def("s") + def("s-") # uattributed, incl building outlines + def("2000400") # buildings as point nodes + def("2000202", , "light moderate reddish yellow") # wall + def("2000206", , "light moderate reddish yellow") # fence + + deflayer("r", "deep gray") # roads and trails + def("r-", , "pale gray") # deemphasize unattributed segments + def("1700201", 3, "black") # road, primary, undivided + def("1700202", 3, "black") # road, primary, divided + def("1700203", 2, "black") # road, primary, one of divided paths + def("1700204", 2, "black") # road, primary, one-way + def("1700205", 2, "black") # road, secondary + def("1700206", 2, "black") # road, secondary + def("1700207", 2, "black") # road, secondary + def("1700208", 2, "black") # road, secondary + def("1700214", 1, "black", "dashed") # ferry route + def("1700218") # road, class 3, divided + def("1700209") # road, class 3, undivided + def("1700402") # entrance ramp + def("r") # unspecified road or trail + def("1700210", , , "dashed") # road, class 4 + def("1700219", , , "dashed") # road, class 4, one-way + def("1700212", , , "dashed") # road, class 5, 4WD + def("1700211", , "dark red", "dashed") # trail + def("1700213", , "dark red", "dashed") # footbridge + + deflayer("t", "dark orange") # railroads + def("t-", , "pale weak orange") # deemphasize unattrib segments + def("t") # unspecified railroad + def("1800201", 2) # railroad main + def("1800201E", 2) # railroad main elevated + def("1800201R", 2) # railroad main on drawbridge + def("1800201T", 2, , "dashed") # railroad main in tunnel + def("1800207", 1, , "dashed") # railroad ferry route + def("1800208", 1) # railroad siding + def("1800209", 2) # railroad yard + def("1800400", 1) # railroad station +$define TRANSIT "dark blue" + def("1800201Y", 2, TRANSIT) # rapid transit rail main + def("1800201EY", 2, TRANSIT) # rapid transit main elevated + def("1800201RY", 2, TRANSIT) # rapid transit main on drawbrg + def("1800201TY", 2, TRANSIT, "dashed") # rapid transit main in tunnel + def("1800202Y", 2, TRANSIT) # rapid transit main in road + def("1800202RY", 2, TRANSIT) # rapid transit, in road on drawbridge + def("1800208Y", 1, TRANSIT) # rapid transit siding + def("1800400Y", 1, TRANSIT) # rapid transit station + + deflayer("u", "light gray") # misc transpt: power, pipe, airport + def("u") + def("u-", , "white-gray") # unattrib segments incl airport runways +$define UTILITY "strong purple-magenta" + def("1900201", 1, UTILITY, "dashed") # petroleum pipeline + def("1900202", 1, UTILITY) # power line + def("1900203", 1, UTILITY) # phone line + def("1900400", 1, UTILITY) # power plant + def("1900401", 1, UTILITY) # substation + def("1900402", 1, UTILITY) # hydro plant + def("1900403", 1, "light gray") # landing strip or airport + def("1900404", 1, "orange") # helipad + def("1900405", 1, "light gray") # launch complex + + deflayer("m", "blue") # survey markers + def("m-", , "pale weak blue") # deemphasize unattributed lines + def("m") + + deflayer("f", "black") # feature labels + def("f") + + deflayer("n", "deep green") # file labels + def("n") + + deflayer("g", ROUTE) # paths (retraced) + + every i := 1 to *pcolors do { # link ea GPS bg/fg/bg set to one list + s := "g" || i + def(s, 2) + attrib[s || "b"].segs := attrib[s].segs + attrib[s || "f"].segs := attrib[s].segs + } + + return +end + + + +# deflayer -- define layer code and default color for subsequent defs + +global layercode, layercolor + +procedure deflayer(lcode, color) + layercode:= lcode + layercolor := color + return +end + + + +# def(key, width, color, style) -- define style info for code or attribute +# +# default width is 1 +# default color is as last set by deflayer() +# default style is "solid" +# +# a key of "x" matches undefined attributes of layer x +# a key of "x-" matches segments without attributes +# +# a width of -1 means "don't draw" + +procedure def(key, width, color, style) + static seq + initial seq := 0 + + /width := 1 + /color := layercolor + /style := "solid" + attrib[key] := attrec(seq +:= 1, layercode, key, width, color, style, []) + return +end + + + +# setdraw(arec) -- set color, linewidth, linestyle based on attribute record +# +# fails if width is negative, meaning that drawing is to be suppressed + +procedure setdraw(arec) + if arec.width < 0 then + fail + WAttrib("fg=" || arec.color, + "linewidth=" || arec.width, "linestyle=" || arec.style) + return +end |