/* Given coverages of contour lines, spot elevations, and rivers,
/* we want river slope, which incorporates an elevation for river
/* junctions dervied by interpolating elevation along the main stem.
/* vertical and xy dimensions may not match, but zfactor makes slope correct


/* Before running this aml:
  /*  Run flipper.aml to create the arc/node cover wholerivs.
/* This aml creates splitrivs and rivspots.
/* It creates and deletes splitconts, tmpcov, and tmpcov2.
/* If the grid %firstgrid% (defined below) does not exist, it will
/* be created at a cell size of %cellsize% units.
/* This assumes that x,y, and z coords are in the same units.
/* For dangling river segments, endpoint elevations are assigned from
/* the grid.  If that value is implausible (dz < 0 or dz > interval, we
/* assume that the dangling segment has the same slope as its neighbor.
/* Note the constants which are set in the next few lines.
/* version 2.0 is based on vresion 1.08, rewritten to match flipper
submit off
&set concov contours
&set netcov wholerivs
&set firstgrid FIRSTGRID
&set elv_item DXF-ELEVATION
&set cellsize 5  /* for grid used to interpolate dangling node elevation
&set interval 20 /* intervals on contour files 
&set maxdangle 2
&set zfactor 1

&TYPE topgrida.aml, Version 2.01 Wed Sep 25 17:00:01 PDT 1996
&type Using the following constants:\concov %concov%\netcov %netcov%\elv_item %elv_item%\cellsize %cellsize%\interval %interval%\maxdangle %maxdangle%\zfactor %zfactor%

&if ^ [exist %concov% -cover] &then
  &return The coverage %concov% does not exist here.
&if ^ [exist %netcov% -cover] &then
  &return The coverage %netcov% does not exist here.
&if ^ [iteminfo %concov% -ARC JORDER -exists] &then
  additem %concov%.aat %concov%.aat JORDER 4 4 I
&if ^ [iteminfo %concov% -ARC SHREVE -exists] &then
  additem %concov%.aat %concov%.aat SHREVE 4 4 I
&if ^ [iteminfo %concov% -ARC STRAHLER -exists] &then
  additem %concov%.aat %concov%.aat STRAHLER 2 2 I
&if ^ [iteminfo %netcov% -ARC %elv_item% -exists] &then
  additem %netcov%.aat %netcov%.aat %elv_item% 8 8 N 2 

arcedit
ec %netcov%;ef arc;sel all;put tmpcov
ec %concov%;ef arc;sel all;put tmpcov;yes
quit

clean tmpcov tmpcov2 %maxdangle% .001 line
kill tmpcov
 
Arcedit
ec tmpcov2
ef arc
sel all
calc tmpcov2-id = tmpcov2#
sel strahler eq 0
put splitconts  
nselect
put splitrivs   
QUIT
y
y
&type The coverages splitrivs and splitconts have been created.

/* going to the use of NEAR instead of tracking topology
renode  splitrivs
nodepoint splitrivs rivspots
near rivspots splitconts line .001
/* now we have to go back to topology
INDEXITEM splitconts.aat splitconts#
INDEXITEM rivspots.pat   splitconts#
INDEXITEM splitrivs.aat FNODE#
INDEXITEM splitrivs.aat TNODE#
INDEXITEM rivspots.pat  RIVSPOTS#
/* These are all the relates for this aml
relate add
scont;   splitconts.aat;INFO;splitconts#;splitconts#;linear;ro
topt;    rivspots.pat;  INFO;TNODE#;     RIVSPOTS#;  linear;ro
frompt;  rivspots.pat;  INFO;FNODE#;     RIVSPOTS#;  linear;ro
[unquote '']
tables
additem rivspots.pat elv  4    12     F      3
sel rivspots.pat
calc elv = -99
resel scont//%elv_item% ge 0  /* otherwise we lose the -99
calc elv = scont//%elv_item%
&if ^ [iteminfo splitrivs -ARC uid -exists] &then /* unique ID
  additem splitrivs.aat uid 4 5 B
&if ^ [iteminfo splitrivs -ARC elvf -exists] &then
  additem splitrivs.aat elvf 4 10 f 3
&if ^ [iteminfo splitrivs -ARC elvt -exists] &then
  additem splitrivs.aat elvt 4 10 f 3 
&if ^ [iteminfo splitrivs -ARC dz -exists] &then
  additem splitrivs.aat dz 4 10 f 3 
&if ^ [iteminfo splitrivs -ARC slope -exists] &then
  additem splitrivs.aat slope 4 10 f 3
sel splitrivs.aat
calc elvf = -99
calc elvt = -99
resel frompt//elv ge 0
calc elvf = frompt//elv
sel splitrivs.aat
resel topt//elv ge 0
calc elvt = topt//elv
sel splitrivs.aat
calc uid = splitrivs-id
quit  /* tables


kill tmpcov2
&if ^ [exist %firstgrid% -grid] &then; &do
  topogrid %firstgrid% %cellsize%   /* this require setenv GRIDALLOCSIZE 177
  contour %concov% %elv_item%
  /* (if you have point coverage)  point spot %elv_item%
  end
  &if ^ [exist %firstgrid% -grid] &then
    &return The grid %firstgrid% was not built.  Aborting.
 &end
&else
  &type Using exisiting version of %firstgrid%
latticespot %firstgrid% rivspots

/* here comes the tricky part, now in tables
tables
sel splitrivs.aat
&set oldnum 9999
&label hugeloop
resel elvf < 0 /* sub-top junctions arcs (plus dangling sources.)
resel SHREVE > 1 /* exclude dangling sources and uncoded streams
resel elvt > 0
&set num [show number select]
&type looking at %num% arcs
&if %num% eq %oldnum% &then
  &return ABORTING: infinite hugeloop
&set oldnum %num%

/* Except for complex river topology between contours ,
/* the aml will drop out of the loop on the 2nd iteration.
&if  %num% = 0 &then; &goto interpdone
UNLOAD tricktmp $recno # init
&sys cat -n tricktmp | tr -d ' ' | tr '\011' ' ' | sed 's/^/\&set .id/' >! tricktmp2
&RUN tricktmp2
/*T &do index = 1 &to %num% /* using global vars for debugging
/*T   &s .id%index% = [show arc [show select %index%] item UID]
/*T &end

&s .index 0
&label bigloop /* I use gotos because aml nested doloops are unreliable
  &s .index [calc %.index% + 1 ]
  &set .id = [value .id%.index%]
  &TYPE index %.index%  %.id%
  sel splitrivs.aat
  resel $recno = %.id%
  /* we are now at a contour line or a previously coded node
  &s .totlength [show record %.id% item LENGTH]
  &s .stringnum = 1
  &s .stringid%.stringnum% %.id%
  &s .lowelv [show record %.id% item elvt]
  &label searchup
  sel splitrivs.aat
  resel TNODE# = [show record %.id% item FNODE#] /* all upstream arcs
  &s .unum [show number select]     /* ususally 2
  &if %.unum% = 0 &then ; &goto endloop  /* dangler
  UNLOAD tricktmp $recno # init
  &sys cat -n tricktmp | tr -d ' ' | tr '\011' ' ' | sed 's/^/\&set .uid/' >! tricktmp2
  &RUN tricktmp2
/*T   &do .uind = 1 &to %.unum%
/*T      &s .uid%.uind% = [show arc [show select %.uind%] item UID]
/*T   &end
  &s .shrevemax [show record %.uid1% item shreve]
  &s .lengmax [show record %.uid1% item LENGTH]
  &if [show record %.uid1% item elvf] le 0 &then
    &s .lengmax 0 /* v. 1.04
  &s .maxind 1
  &do .uind = 2 &to %.unum% /* ususally 2
    &s .shreve = [show record [value .uid%.uind%] item shreve] 
    &if %.shreve% eq %.shrevemax% &then ; &do   /* if equal shreve,
     &if [show record [value .uid%.uind%] item LENGTH] > %.lengmax% and  ~
      [show record [value .uid%.uind%] item elvf] gt 0 &then ; &do   /* v. 1.04
       &s .shrevemax %.shreve%
       &s .maxind %.uind%
     &end
    &end
    &if %.shreve% > %.shrevemax% &then ; &do
       &s .shrevemax %.shreve%
       &s .maxind %.uind%
    &end
  &end 
/*T  sel UID = [value .uid%.maxind%] /* the upstream arc with highest shreve
  sel splitrivs.aat
  resel $recno = [value .uid%.maxind%]
  &s .totlength = [calc %.totlength% + [show record [value .uid%.maxind%] item LENGTH]]
  &s .stringnum = %.stringnum% + 1
  &s .stringid%.stringnum% [value .uid%.maxind%]
  &s .highelv = [show record [value .uid%.maxind%] item elvf]
  &if  %.highelv% < 0 &then;&DO
    &TYPE searching up %.highelv% %.stringnum%
    &set .id [value .uid%.maxind%]
    &goto searchup
  &end
  /* We have traced up to a contour
  &type Filling nodes on %.stringnum% arcs from %.lowelv% to %.highelv%
  &s subleng 0
  &s .delv = %.highelv% - %.lowelv%
  &do .uind = 1 &to [ calc %.stringnum% - 1 ]
    sel splitrivs.aat
    resel $recno = [value .stringid%.uind%]
    &s subleng = [ calc %subleng% + [show record [value .stringid%.uind%] item LENGTH] ]
    &s .here = [ calc %subleng% / %.totlength% ]
    &s .elv = [ calc %.here% * %.delv% + %.lowelv% ]
    &if [show record [value .stringid%.uind%] item elvf] < 0 &then
      CALC elvf = %.elv%
    &ty %.index% calculated elv of %.elv% on arc [value .stringid%.uind%]
/*     &if %.uind% < %.stringnum% &then ; &do
      &set upnode = [show record [value .stringid%.uind%] item FNODE#] /* all upstream arcs
      sel splitrivs.aat
      resel TNODE# = %upnode%
      CALC elvt = %.elv%
/*    &end
  &end
&label endloop
/* &end /* .index
&if %.index% < %num% &then; &goto bigloop
 &goto hugeloop
&label interpdone
quit /* out of tables
/* end of the tricky part

ARCEDIT
ec splitrivs
ef arc
&type fixing  to nodes
sel  elvf - topt//SPOT gt %interval% /* the grid would give us a large dz
asel  elvf le topt//SPOT /* the grid would give us a negative slope
resel elvt lt 0 /* dangling mouths
&set num [show number select]
&if %num% > 0 &then;&do
  &do .index = 1 &to %num%
    &s .id%.index% = [show arc [show select %.index%] item UID]
  &end
  &do .index = 1 &to %num%
    &set .id = [value .id%.index%]
    sel UID = %.id%
    sel connect
    &if [show number select] ne 1 &then
      &type WARNING: skipped assigning elevation to arc %.id%, [show number select] connected arcs.
    &else;&do
      &set dz = [show arc [show select 1] item elvf] - [show arc [show select 1] item elvt]
      &set slope = [calc %dz%  / [show arc [show select 1] item length]]
/* no correction, everything is meters
      sel UID = %.id%
      calc elvt = elvf - %slope% * length
      calc dz = elvf - elvt
      &set dz = [show arc [show select 1] item dz]
      &if %dz% > %interval% &then;&do
        calc dz = %interval% * .8
        calc elvt = elvf - dz
      &end
      &if %dz% le 0 &then;&do
        calc dz = %interval% * .2
        calc elvt = elvf - dz
      &end
    &end
  &end
&end
sel elvt le 0
&if [show number select] > 0 &then
  calc elvt = topt//SPOT

&type fixing  from nodes
sel elvt ge frompt//SPOT
asel frompt//SPOT - elvt gt %interval%
resel elvf lt 0 /* dangling sources
&if [show number select] > 0 &then;&do
  &set num [show number select] 
  &do .index = 1 &to %num%
    &s .id%.index% = [show arc [show select %.index%] item UID]
  &end
  &do .index = 1 &to %num%  
    &set .id = [value .id%.index%]
    sel UID = %.id% 
    sel connect 
    &if [show number select] ne 1 &then 
      &type WARNING: skipped assigning elevation to arc %.id%, [show number select] connected arcs.
    &else;&do 
      &set dz = [show arc [show select 1] item elvf] - [show arc [show select 1] item elvt]
      &set slope = [calc %dz%  / [show arc [show select 1] item length]] /* no correction, everything is meters 
      sel UID = %.id%
      calc elvf = elvt + %slope% * length 
      calc dz = elvf - elvt
      &set dz = [show arc [show select 1] item dz]
      &if %dz% > %interval% &then;&do
        calc dz = %interval% * .8
        calc elvf = elvt + dz
      &end
      &if %dz% le 0 &then;&do
        calc dz = %interval% * .2
        calc elvf = elvt + dz
      &end
    &end
  &end
&end
sel elvf lt 0
&if [show number select] > 0 &then 
  calc elvf = frompt//SPOT

select all
calc dz = elvf - elvt
/* vertical and xy dimensions may not match, but zfactor makes slope correct
calc slope = dz  / length   * %zfactor%
quit
y
y

/* delete unneeded coverages
kill splitconts
/* kill %firstgrid%