/* Routine to output a series of river long profiles past stations /* AND FIX THEM /* This macro runs under the ARC/INFO GRID module. /* It makes unix calls "&sys echo" and "&sys cat /dev/null", so people with /* Redmond-based operating systems should make sure that their unix emulators /* are in good order. /* /* Copyleft Harvey Greenberg (hgreen@u.washington.edu) University of Washington /* and Rolf Aalto (aalto@u.washington.edu) University of Washington /* First we trace upstream without writing to output. We trace /* until we reach the cell adjacent to the drainage divide, or until the /* discharge is less than a specified threshhold. /* Then we trace downhill, fixing elevations. /* things to do: corrections for latlong grids, /* accepting arguments instead of having a dialog, running recursively to /* fix an entire watershed (a kludgy solution would be to make a list of /* channel heads and let this program trace down to the basin outlet n times) /* The program begins with a dialog. Hit a carriage return to accept the /* author's default values. /* For this version: input is /* name of DEM as it existed before "sinks" were filled. /* name of flowdirection grid /* name of discharge grid (if none, repeat the name of flowaccumulation grid) /* name of flowaccumulation grid /* minimum discharge value for tracing upstream /* name of file of "basinname" x y /* max number of pixels to follow downstream from starting point /* minimum stream slope (multiply by factor if zunits ne xyunits) /* misc grid to sample (optional) /* another misc grid to sample (optional) /* each line of the input file contains a string and X and Y, e.g. /* b01,435720,2454460 /* b02,434500,2453660 /* ... /* The coordinate must fall on a calculated stream cell. If you choose a channel /* head, the program will execute faster. If you choose a downstream point, the /* program will first trace uphill in the direction of maximum discharge (or /* until discharge falls below a certain value), then trace back downhill. You /* can specify the maximum distance to trace downhill below the starting point. /* The program assumes the "sinks" in the DEM are spurious artifacts. If a cell /* downstream (as calculated from a filled DEM) has a higher elevation, that cell /* is lowered so that the river maintains a minimum slope. This fix is reflected /* both in the profile, and in a file of XYZ points that can be used in a separate /* step to fix your DEM. /* A brief summary is writen to the file "profile.report". Old files by this /* name will be overwritten. /* The program creates a file named "fixer.gen" contains XYZ values that can be /* used to improve your DEM. Old files by this name will be overwritten. /* For each input line, a file is written containing a line for /* each grid cell on the channel trace. The name of this file is the fist item /* in the input file, with the extension _prf.xy (by default) appended. Each line /* contains: /* distance downstream from the channelhead /* fixed DEM value (usually the same as the next field) /* elevation from the input DEM /* contributing area /* discharge /* fisrt misc value (if requested by user) /* second misc value (if requested by user) /* X /* Y /* check to make sure program is in GRID &if [show program] ne GRID &then grid &type proffix.aml, version 1.00 /*Input DEM, flowdirections grid, stations, and output. &do &until [query 'Continue?' .FALSE.] &s dem = [response 'Input DEM name (elv)' elv] /* must be well masked? &s flowd = [response 'Input flowdirection grid (optional)(FLOWDIR)' FLOWDIR] &ty ' The discharge grid is used to trace uphill.\~ If you don't have one, use flowaccumulation.' &s disch = [response 'Input discharge grid (discharge)' discharge] &s flowak = [response 'Input flowaccumulation grid (optional )(FLOWACC)' FLOWACC] &s thresh = [response 'Minumum discharge to trace up channel (0)' 0] &s stations_f = [response 'Input station file (outlets.txt)' outlets.txt] &s count = [response '# of pixels in downstream (99999)' 99999] &s outfile = [response 'Outfile name appendix? (_prf.xy)' _prf.xy] &s dz = [response 'Minimum stream slope (.0001)' .0001] &s misc1 = [response 'Misc grid #1 (none)' none] &s misc2 = [response 'Misc grid #2 (none)' none] &type Entered values are %dem%, %flowd%, %disch%, %flowak%, %stations_f%, %count%, %outfile%, %dz%, %misc1%, and %misc2%. &end &type This program runs slowly. Each '.' echoed to the screen represents one grid cell. /* Read in *.txt Data Files &sys echo opening: &sys ls -l %stations_f% &s readu = [open [unquote %stations_f%] status1 -read] &ty Opened %stations_f% with status %status1% setwindow %dem% &s sp ' ' &describe %dem% /* input information for dem &s dx = %GRD$DX% &s dd [calc %dx% * 1.414213462] &s dy = %GRD$DY% &s xmax = %GRD$XMAX% &s xmin = %GRD$XMIN% &s ymax = %GRD$YMAX% &s ymin = %GRD$YMIN% &s ak ' ' &s m1 ' ' &s m2 ' ' &type Beginning main program /* begin main loop &s status2 0 &label bigloop &s textline = [unquote [read %readu% status2 ]] &ty %textline% %status2% &if %status2% ne 0 &then; &ret &s station [extract 1 %textline%] &if [exists %station%%outfile% -file] &then &type [delete %station%%outfile%] &sys cat /dev/null >! %station%%outfile% /*prepare file and coords for inner loops &sys cat /dev/null >! fixer.gen &sys cat /dev/null >! profile.report &type Processing %station% : &describe %station% &s xsta = [extract 2 %textline%] &s ysta = [extract 3 %textline%] /* Going upstream just to find starting point. &s x = %xsta% /* reset variables for next loop &s y = %ysta% &s i = 0 &s dist 0.0 &s inbounds = .TRUE. &type /& ---> upstream &do &while %inbounds% /*upstream loop &s dc = [show cellvalue %disch% %x% %y% ] /* read in flow accumulation data &s c1 = [show cellvalue %disch% [calc %x% + %dx%] %y%] &s c2 = [show cellvalue %disch% [calc %x% + %dx%] [calc %y% - %dy%]] &s c3 = [show cellvalue %disch% %x% [calc %y% - %dy%]] &s c4 = [show cellvalue %disch% [calc %x% - %dx%] [calc %y% - %dy%]] &s c5 = [show cellvalue %disch% [calc %x% - %dx%] %y%] &s c6 = [show cellvalue %disch% [calc %x% - %dx%] [calc %y% + %dy%]] &s c7 = [show cellvalue %disch% %x% [calc %y% + %dy%]] &s c8 = [show cellvalue %disch% [calc %x% + %dx%] [calc %y% + %dy%]] /* reset cells that don't flow into current cell &if [type %c1%] gt 0 &then;&s c1 0;&else &if [show cellvalue %flowd% [calc %x% + %dx%] %y%] ne 16 &then &s c1 = 0 &if [type %c2%] gt 0 &then;&s c2 0;&else &if [show cellvalue %flowd% [calc %x% + %dx%] [calc %y% - %dy%]] ne 32 &then &s c2 = 0 &if [type %c3%] gt 0 &then;&s c3 0;&else &if [show cellvalue %flowd% %x% [calc %y% - %dy%]] ne 64 &then &s c3 = 0 &if [type %c4%] gt 0 &then;&s c4 0;&else &if [show cellvalue %flowd% [calc %x% - %dx%] [calc %y% - %dy%]] ne 128 &then &s c4 = 0 &if [type %c5%] gt 0 &then;&s c5 0;&else &if [show cellvalue %flowd% [calc %x% - %dx%] %y%] ne 1 &then &s c5 = 0 &if [type %c6%] gt 0 &then;&s c6 0;&else &if [show cellvalue %flowd% [calc %x% - %dx%] [calc %y% + %dy%]] ne 2 &then &s c6 = 0 &if [type %c7%] gt 0 &then;&s c7 0;&else &if [show cellvalue %flowd% %x% [calc %y% + %dy%]] ne 4 &then &s c7 = 0 &if [type %c8%] gt 0 &then;&s c8 0;&else &if [show cellvalue %flowd% [calc %x% + %dx%] [calc %y% + %dy%]] ne 8 &then &s c8 = 0 /* order discharge and select largest &s dc2 = [extract 8 [sort %c1% %c2% %c3% %c4% %c5% %c6% %c7% %c8% -NUMERIC]] /* &ty xy = %x%,%y%, cs are %c1% %c2% %c3% %c4% %c5% %c6% %c7% %c8%, dc2 = %dc2% &if %dc2% eq 0 &then ; &goto continue &select %dc2% /* determine which cell this is, flowdirectionwise &when %c1% &s us = 1 &when %c2% &s us = 2 &when %c3% &s us = 4 &when %c4% &s us = 8 &when %c5% &s us = 16 &when %c6% &s us = 32 &when %c7% &s us = 64 &when %c8% &s us = 128 &otherwise &goto continue &end &if %us% in { 1, 2, 128 } &then &s x = [calc %x% + %dx% ] &if %us% in {8, 16, 32} &then &s x = [calc %x% - %dx% ] &if %us% in {32, 64, 128} &then &s y = [calc %y% + %dy% ] &if %us% in {2, 4, 8} &then &s y = [calc %y% - %dy% ] &s ds = [show cellvalue %flowd% %x% %y% ] &if %ds% in {1,4,16,64} &then /* harvey added this (4 lines) &s dist [calc %dist% - %dx%] &else &s dist [calc %dist% - %dd%] &if %x% > %xmax% or %x% < %xmin% &then &s inbounds = .FALSE. &if %y% > %ymax% or %y% < %ymin% &then &s inbounds = .FALSE. &if %dc2% le %thresh% &then &s inbounds = .FALSE. /* end if last upstream pixel &s i [calc %i% - 1] &end /*end of upstream loop &label continue &type traced %dist% units upstream to top. &sys echo traced %dist% units upstream to top. >> profile.report &type Now tracing downhill in earnest from %x% %y%. &s dc = [show cellvalue %disch% %x% %y% ] &s inbounds = .TRUE. &s dist 0.0 &s zold = 99999 &s olddir 1 &type /& ---> downstream &do &while %inbounds% /*downstream loop &s z = [show cellvalue %dem% %x% %y% ] &if [type %z%] gt 0 &then ;&goto breakloop /* hit NODATA &s ds = [show cellvalue %flowd% %x% %y% ] &if [type %ds%] gt 0 &then ;&goto breakloop /* hit NODATA &s dc = [show cellvalue %disch% %x% %y% ] &if [exist %flowak% -grid] &then ;&s ak = [show cellvalue %flowak% %x% %y% ] &if [exist %misc1% -grid] &then ;&s m1 = [show cellvalue %misc1% %x% %y% ] &if [exist %misc2% -grid] &then ;&s m2 = [show cellvalue %misc2% %x% %y% ] &if %ds% in {1,4,16,64} &then &s dxy %dx% &else &s dxy %dd% &if %z% ge %zold% &then &s fixed [calc %zold% - ( %dz% * %dxy% ) ] &else &s fixed = %z% &sys echo %dist% %fixed% %z% %ak% %dc% %m1% %m2% %x% %y% >> %station%%outfile% &if %z% ne %fixed% &then ; &sys echo 1 %x% %y% %fixed% >> fixer.gen &s zold = %fixed% &s i = [calc %i% + 1] &if [type %ds%] gt 0 &then; &goto breakloop &if %ds% in {1,4,16,64} &then /* to the downward cell &s dxy %dx% &else &s dxy %dd% &if %ds% in { 1, 2, 128 } &then &s x = [calc %x% + %dx% ] &if %ds% in {8, 16, 32} &then &s x = [calc %x% - %dx% ] &if %ds% in {32, 64, 128} &then &s y = [calc %y% + %dy% ] &if %ds% in {2, 4, 8} &then &s y = [calc %y% - %dy% ] &s dist [calc %dist% + %dxy%] &if %x% > %xmax% or %x% < %xmin% &then &s inbounds = .FALSE. &if %y% > %ymax% or %y% < %ymin% &then &s inbounds = .FALSE. &if %i% > %count% &then &s inbounds = .FALSE. &sys echo -n . &end &label breakloop &goto bigloop &type That is all folks!