############################################################################
#
#	File:     tgrprep.icn
#
#	Subject:  Program to prepare TIGER line chains
#
#	Authors:  Gregg M. Townsend and William S. Evans
#
#	Date:     June 9, 2000
#
############################################################################
#
#  Tgrprep writes files of "line chain" data extracted from Census Bureau
#  TIGER/Line data files.  The main purpose of this is to prepare
#  input for the "tgrmap" program.
#
#  Usage:  tgrprep rec1file rec2file
#
#	rec1file:  tgr*.rt1 file containing Type 1 (chain) data
#	rec2file:  tgr*.rt2 file containing Type 2 (shape point) data
#
############################################################################
#
#  Output consists of:
#
#  Line 1: the smallest longitude and latitude found in the Type 1 file
#  Line 2: the largest  longitude and latitude found in the Type 1 file
#  Line 3-n: chain specifications, one per line, in the following format
#
#	(3 chars)  census feature class code (CFCC)
#	(1 char)   boundary code (see below)
#	(varies)   optional feature name, delimted by '|' (see below)
#	(4 chars)  max dimension (latitude or longitude units), max 9999
#	(7 chars)  starting longitude, fraction E of Greenwich meridian
#	(7 chars)  starting latitude, fraction S of North Pole
#	(4 chars)  delta longitude to first point, same units, plus 5000
#	(4 chars)  delta latitude to first point, same units, plus 5000
#  followed by any number (zero or more) of additional delta pairs. 
#  Output lines may be arbitrarily long.
#
#  Boundary codes are:
#	9	national boundary  (not used)
#	8	state boundary  (not used)
#	7	county boundary
#	5	city limits
#	0	other, unknown, not a boundary
#
#  Feature name is the concatenation of the following TIGER fields
#	"Feature Direction, Prefix"
#	"Feature Name"
#	"Feature Type"
#	"Feature Direction, Suffix"
#  The concatenation is surrounded by vertical bars "|" unless it is empty
#  (all spaces).  Any "|" within a TIGER field is replaced by "!".
#
#  For input formats and the definition of CFCC codes, see
#	TIGER/Line Files, 1998 Technical Documentation
#	Bureau of the Census, Washington, DC, 1998.
#	http://www.census.gov/geo/www/tiger/tiger98.pdf
############################################################################


global minlon, maxlon	# min/max longitude seen (in input terms)
global minlat, maxlat	# min/max latitude seen (in input terms)

global curlon, curlat	# current longitude/latitude for output
global deltas		# string of deltas for output


procedure main(args)
   local details, file1, file2, n

   *args = 2 | stop("usage: ", &progname, " rec1file rec2file")
   file1 := open(args[1]) | stop("can't open ", args[1])
   file2 := open(args[2]) | stop("can't open ", args[2])

   write(&errout, "prescanning ", args[1])
   n := llrange(file1)
   write(&errout, right(n, 10), " chain records")
   write("        ", rz(cvlon(minlon)), rz(cvlat(maxlat)))
   write("        ", rz(cvlon(maxlon)), rz(cvlat(minlat)))

   write(&errout, "prescanning ", args[2])
   details := dtindex(file2)
   write(&errout, right(*details, 10), " supplemental sets")

   write(&errout, "scanning ", args[1])
   n := scan(file1, file2, details)
   write(&errout, right(n, 10), " supplements used")
   write(&errout, "done")
end


#  scan(file1, file2, details) -- scan records and write output.
#
#  returns the number of supplements referenced.

procedure scan(file1, file2, details)
   local line, tlid, cfcc, lon, lat, n, l, r
   local startlon, startlat, endlon, endlat, dim, bound, fename

   n := 0
   seek(file1, 1)
   while line := read(file1) do line ? {
      ="1" | next
      tab(6)
      tlid := move(10)
      tab(18)

      fename := ""
      fename ||:= " " || ("" ~== trim(move(2) \ 1))	# direction prefix
      fename ||:= " " || ("" ~== trim(move(30) \ 1))	# name
      fename ||:= " " || ("" ~== trim(move(4) \ 1))	# type
      fename ||:= " " || ("" ~== trim(move(2) \ 1))	# direction prefix
      if fename ~== "" then 
         fename := "|" || map(fename[2:0], "|", "!") || "|"

      tab(56)
      cfcc := move(3)
      bound := "0"

      tab(135)
      l := move(3)		#left county code
      r := move(3)		#right county code
      if l ~== r then		#different --> county boundary
         bound := "7"
      else {
         tab(161)
         l := move(5)		#left city code
         r := move(5)		#right city code
         if l ~== r then	#different --> city boundary
            bound := "5"
         }

      tab(191)
      startlon := curlon := minlon := maxlon := cvlon(move(10))
      startlat := curlat := minlat := maxlat := cvlat(move(9))
      endlon := cvlon(move(10))
      endlat := cvlat(move(9))

      deltas := ""
      if seek(file2, \details[tlid]) then {
         n +:= 1
         while line := read(file2) do line ? {
            tab(6)
            =tlid | break
            tab(19)
            every 1 to 10 do
               drawto(cvlon(0 ~= move(10)), cvlat(0 ~= move(9)))
            }
         }
      drawto(endlon, endlat)

      dim := startlon - minlon
      dim <:= maxlon - startlon
      dim <:= startlat - minlat
      dim <:= maxlat - startlat
      dim >:= 9999

      write(cfcc, bound, fename,
         right(dim, 4), rz(startlon), rz(startlat), deltas)
      }
   return n
end


#  drawto(lon, lat) -- append deltas, updating curlon/curlat

procedure drawto(lon, lat)
   local dlon, dlat

   dlon := lon - curlon
   dlat := lat - curlat

   if abs(dlon | dlat) >= 5000 then {
      drawto(curlon + dlon / 2, curlat + dlat / 2)
      drawto(lon, lat)
      }
   else {
      deltas ||:= rz(dlon + 5000, 4)
      deltas ||:= rz(dlat + 5000, 4)
      curlon := lon
      curlat := lat
      minlon >:= lon
      maxlon <:= lon
      minlat >:= lat
      maxlat <:= lat
      }
   return
end



#  rz(v, n) -- right-justify value in n digits with zero fill

procedure rz(v, n)
   /n := 7
   return right(v, n, "0")
end



#  cvlon(n) -- convert longitude to output form
#
#  (Fraction of circle east of Greenwich, as 0000000 to 9999999).

procedure cvlon(n)
   static m
   initial m := 9999999 / 360.0 / 1000000

   n := integer(n)
   if n < 0 then
      n +:= 360000000
   return integer(m * n)
end


#  cvlat(n) -- convert latitude to output form
#
#  (Fraction of semicircle south of North Pole, as 0000000 to 9999999).

procedure cvlat(n)
   static m
   initial m := 9999999 / 180.0 / 1000000
   return integer(m * (90000000 - n))
end



#  dtindex(f) -- return table of record indices by TLID from file f

procedure dtindex(f)
   local details, line, w

   details := table()
   seek(f, 1)
   while (w := where(f)) & (line := read(f)) do line ? {
      ="2" | next
      tab(6)
      /details[move(10)] := w
   }
   return details
end



#  llrange(f) -- scan f to set min/max lon/lat, returning record count

procedure llrange(f)
   local line, n, lon, lat

   minlon := +180000000
   maxlon := -180000000
   minlat := +90000000
   maxlat := -90000000
   n := 0

   seek(f, 1)
   while line := read(f) do line ? {
      ="1" | next
      n +:= 1
      tab(191)
      every 1 | 2 do {
         lon := integer(move(10))
         lat := integer(move(9))
         minlon >:= lon 
         maxlon <:= lon 
         minlat >:= lat 
         maxlat <:= lat 
         }
      }
   return n
end