/* aml to test loading algorithms Michael Silberbauer Jan 1994 /* Edited by michael Feb 1994 08:00 Thursday /* revisited Michael Silberbauer Feb 2005 /* removed large AML "arrays" Michael Silberbauer Mar 2005 /* /* /prjws8/users/michael/aml/load.aml /* ------------------------------------------------------------------------------------ /* This AML was developed as part of the WaterMarque water quality assessment /* project of the Institute for Water Quality Studies /* Department of Water Affairs and Forestry /* Private Bag X313, PRETORIA, 0001 South Africa. /* tel: (012) 808-0374 fax: (012) 808-0338 /* The AML is in the public domain, but copyright is held by the Department /* of Water Affairs and the CSIR. The AML should not be re-distributed /* without permission. Where the AML is used for publication or reporting, /* the source should be acknowledged. /* Purpose: to experiment with load calculation algorithms. /* December 1993: Michael Silberbauer wrote first version /* March 2005 : Michael Silberbauer dusted off and updated paths /* 2005-07-14 : Michael Silberbauer variable arguments /* 2007-12-03 : Michael Silberbauer added _%ISOdate1%_%ISOdate2% to filename /* 2007-12-14 : Michael Silberbauer inserted small main routine for batch /* 2009-03-30 : Michael Silberbauer reselect %tlt% logq and logc > 0 for plotting /* 2009-04-02 : Michael Silberbauer colour-code load to regression points /* 2009-04-03 : Michael Silberbauer fix text angle on regression line /* 2009-04-06 : Michael Silberbauer cosmetic changes to load listing /* 2009-04-07 : Michael Silberbauer change >= 0 to > 0 for flow (was so for concentration) /* ------------------------------------------------------------------------------------ &args Station DateStart DateEnd ChemItem LeBug &if [null %Station%] &then &sv Station A2H027 &else &sv Station [translate %Station%] &if [null %DateStart% ] &then &sv DateStart 2002-01-01 &if [null %DateEnd% ] &then &sv DateEnd 2002-12-31 &if [null %ChemItem%] &then &sv ChemItem = TDS &else &sv ChemItem = [translate %ChemItem%] &if [null %LeBug%] &then &sv Debug = .FALSE. &else &sv Debug = .TRUE. &if [length %DateStart%] = 4 &then &sv Year1 = %DateStart% &else &sv Year1 = [substr %DateStart% 1 4] &if [length %DateEnd% ] = 4 &then &sv Year2 = %DateEnd% &else &sv Year2 = [substr %DateEnd% 1 4] &call main &if [show program] = ARCPLOT &then quit &return /* - - - - - - - - - - - - - - - - - &routine main &sv newColumn = .FALSE. &sv fill = .true. &sv closeall = [close -all] &sv FlowStn = %Station%A01 &sv Flow-File = /hri/db/iwqs/db/flow/[locase %FlowStn%] &if [exists %Flow-File% -file] &then &type Found %Flow-File% &else &stop Sorry, no flow data. /*&sv Chem-File = /hri/db/wmrq/wmdata/wq/inorganic.dat &sv Chem-File = /hri/db/wmrq/wmdata/wq/MACRO080815.dat &sv StnFile = /hri/db/cover/s-africa/nms_wms &sv tf = FLOWLL.TMP &sv tc = CHEM.TMP &sv tl = FLOWCHEM.TMP &sv tlt = FLOWCHEMT.TMP &sv Catchment = /hri/db/cover/s-africa/hca_4 &sv Rivers = /hri/db/cover/s-africa/wri_500 &sv Lake = /hri/db/cover/s-africa/wla_500 &sv ChemStn = %Station%Q01 &sv Primary = [substr %ChemStn% 1 1] &sv Secondary = [substr %ChemStn% 1 2] &type Chem = %ChemStn%, flow = %FlowStn% &call CheckDate ArcPlot /* Create a plotfile name, removing any funny characters: &sv PlotFile = [locase %Station%L_%ChemItem%_%ISOdate1%_%ISOdate2%] &sv PlotFile = [subst %PlotFile% ' ' ''] &sv PlotFile = [subst %PlotFile% '/' ''] &sv PlotFile = [subst %PlotFile% '(' ''] &sv PlotFile = [subst %PlotFile% ')' ''] &sv PlotFile = [subst %PlotFile% '.' ''] &sv PlotFile = [subst %PlotFile% '-' ''] &sv TextFile = %PlotFile%.txt &sv tlx = %PlotFile%q.txt &if [exists %TextFile% -file] &then &sv dtf = [delete %TextFile% -file] &sys echo %Station% > %TextFile% &sys echo %AML$FULLfile% >> %TextFile% &call ISOdate &sys echo %ISOdate% >> %TextFile% &if [exists load.log -file] &then &sys echo %Station% >> load.log &else &sys echo %Station% > load.log &sys echo %AML$FULLfile% >> load.log &sys echo %ISOdate% >> load.log /*Display 9999 4 Display 1040 %PlotFile% &type creating plotfile %PlotFile% PAGEUNITS cm &call PageSetup TEXTSET font TEXTSYMBOL 14 TEXTSCALE 2 TEXTCOLOR 4 CLEARSELECT %Chem-File% INFO RESELECT %Chem-File% INFO STATION = [quote %ChemStn%] RESELECT %Chem-File% INFO DATE >= %ISOdate1% and DATE <= %ISOdate2% &if [before [show select %Chem-File% INFO] ,] = 0 &then &do &type Sorry, no chem data available %ISOdate1% - %ISOdate2%. &return &end STATISTICS %Chem-File% INFO min NDATE max NDATE end &sv startchem = [calc [show statistic 1 1] ] &sv endchem = [calc [show statistic 2 1] ] &if %Debug% &then &do &type Using existing temporary flow INFO file MOVE %xtxtp% %ytxtp% textsymbol 14;textcolor 1;textscale 1 text [quote Using existing temporary flow INFO file] ul &call nextrow &end &else &call FlowFile RESELECT FLOWL.TMP INFO DATE >= %ISOdate1% and DATE <= %ISOdate2% &if [before [show select FLOWL.TMP INFO] ,] = 0 &then &do &type Sorry, no flow data available %ISOdate1% - %ISOdate2%. &return &end STATISTICS FLOWL.TMP INFO min NDATE max NDATE end &sv startflow = [calc [show statistic 1 1] ] &sv endflow = [calc [show statistic 2 1] ] &if ( ( %endchem% < %startflow% ) OR ( %startchem% > %endflow% ) ) &then &do &type Chemical and flow files do not overlap: &type Giving up... &return &end &sv startday = [max %startflow% %startchem%] &sv endday = [min %endflow% %endchem% ] RESELECT FLOWL.TMP INFO NDATE >= %startday% and NDATE <= %endday% &if [before [show select FLOWL.TMP INFO] ,] = 0 &then &do &type Sorry, no flow data available %ISOdate1% - %ISOdate2%. &return &end &if [exists %tf% -INFO] &then &sv delflow = [delete %tf% -INFO] INFOFILE FLOWL.TMP INFO %tf% DATE NDATE FLOW init /*list %tf% info /*&stop debug stop RESELECT %Chem-File% info NDATE >= %startday% and NDATE <= %endday% RESELECT %Chem-File% info %ChemItem% > 0 RESELECT %Chem-File% info STATION = [quote %ChemStn%] &sv nChems = [before [show select %Chem-File% INFO] ,] &if %nChems% = 0 &then &do &type Sorry, no chem data available %ISOdate1% - %ISOdate2%. &return &end &if [exists %tc% -INFO] &then &sv delchem = [delete %tc% -INFO] INFOFILE %Chem-File% INFO C.TMP STATION DATE TIME NDATE %ChemItem% init /* Ensure only one value per day: &data arc frequency C.TMP CHEM.TMP NDATE END END &end arc joinitem CHEM.TMP C.TMP CHEM.TMP NDATE /*INFOFILE %Chem-File% INFO %tc% STATION DATE TIME NDATE %ChemItem% init listoutput load.log &call LoadCalc MOVE %xtxtp% %ytxtp% textsymbol 14;textcolor 1;textscale 1 text [quote Using chem file %Chem-File%] ul &call nextrow TEXTSCALE 2 textoffset %rowsize% -%rowsize% BOX %xpgq1% %ypgq1% %xpgq2% %ypgq2% MOVE %xpgq1% %ypgq2%; text Flow\m3/sec ul LINECOLOR 2 BOX %xpgc1% %ypgc1% %xpgc2% %ypgc2% MOVE %xpgc1% %ypgc2%; text Conc\%ChemItem%\mg/L ul LINECOLOR 3 BOX %xpgl1% %ypgl1% %xpgl2% %ypgl2% MOVE %xpgl1% %ypgl2%; text Load\%ChemItem%\kg/day ul LINECOLOR 4 BOX %xpgr1% %ypgr1% %xpgr2% %ypgr2% LINECOLOR 5 BOX %xpgm1% %ypgm1% %xpgm2% %ypgm2% LINECOLOR 1 BOX %xtx1% %ytx1% %xtx2% %ytx2% textoffset 0 0 &call ISOdate MOVE %xpg2% %ypg1% TEXTSYMBOL 14 TEXTSCALE 2 TEXT [quote Plotted on %ISOdate%] lr TEXTSCALE 1 &type [show textsize] &call PlotMap /* * * * * * * * * * /* Plot flow /* Set the flow graph extent (force the date axis limits): LINESYMBOL 1 CLEARSELECT %tf% INFO reselect %tf% INFO flow > 0 &sv ndata = [extract 1 [show select %tf% INFO]] STATISTICS %tf% INFO mean FLOW end &sv MeanQ = [calc [show statistic 1 1] * 86400] &sv datafile = %tf% &sv datavar = FLOW &sv Percentile = 50 &call Percentile &data ARC INFO ARC SELECT [translate %tf%] SORT ON NDATE Q STOP q &end &sv MedQ = [calc %PcntVar% * 86400] &type %ndata% flows from %ISOdate1% to %ISOdate2%: Mean = [round %MeanQ%] Median = [round %MedQ%] m3/day GRAPHEXTENT %tf% INFO ndate flow &sv GExMin = [extract 1 [show graphextent]] &sv GEyMin = [extract 2 [show graphextent]] &sv GExMax = [extract 3 [show graphextent]] &sv GEyMax = [extract 4 [show graphextent]] GRAPHLIMITS %xgrq1% %ygrq1% %xgrq2% %ygrq2% &sv GLxcm = [calc %xgrq2% - %xgrq1%] &sv GExDays = [calc %GExMax% - %GExMin%] &sv cm_Day = [calc %GLxcm% / %GExDays%] &sv GLycm = [calc %ygrq2% - %ygrq1%] &sv GEyUnits = [calc %GEyMax% - %GEyMin%] &sv cm_yUnit = [calc %GLycm% / %GEyUnits%] &type Flow (i.e. q) Y max = %GEyMax% GRAPHEXTENT %.startday% 0 %.endday% %GEyMax% UNITS graph LINESYMBOL 1 LINECOLOR 1 AXIS vertical TEXTSCALE 2 AXISRULER TEXTOFFSET [calc -1 * ( 4 * %xsp% ) ] 0 TEXTANGLE 90 /*EXT m3/sec TEXTOFFSET 0 0 TEXTANGLE 0 TEXTSCALE 2 AXIS horizontal &r /rpd2/spek/prjws8/users/michael/aml/AXISDATE M Year M Day &sv deltaplotfile = %tl% &sv deltaplotvar = flow &sv deltaplotsym = fdelta &sv hue = 300 &call deltaPlot GRAPHPOINT %tf% INFO ndate flow 1 /* * * * * * * * * * /* Plot chem /* Set the chemical graph extent (force the date axis): LINESYMBOL 1 clearselect %tl% info reselect %tl% info %ChemItem% > 0 STATISTICS %tl% INFO mean %ChemItem% end &sv MeanC = [show statistic 1 1] &ty Mean %ChemItem% from %ISOdate1% to %ISOdate2% = %MeanC% GRAPHLIMITS %xgrc1% %ygrc1% %xgrc2% %ygrc2% &sv GLycm = [calc %ygrc2% - %ygrc1%] &sv GEyUnits = [calc %GEyMax% - %GEyMin%] &sv cm_yUnit = [calc %GLycm% / %GEyUnits%] GRAPHEXTENT %tl% INFO ndate %ChemItem% &sv ChemMax = [extract 4 [show graphextent]] &type Chem (i.e. c) Y max = %ChemMax% GRAPHEXTENT %.startday% 0 %.endday% %ChemMax% UNITS graph TEXTSCALE 2 AXIS vertical AXISRULER TEXTOFFSET [calc -1 * ( 4 * %xsp% ) ] 0 TEXTANGLE 90 /*AXISTEXT %ChemItem% TEXTOFFSET 0 0 TEXTANGLE 0 &sv deltaplotfile = %tl% &sv deltaplotvar = %ChemItem% &sv deltaplotsym = cdelta &sv hue = 120 &call deltaPlot clearselect %tc% info GRAPHPOINT %tc% INFO ndate %ChemItem% 1 AXIS horizontal &r /rpd2/spek/prjws8/users/michael/aml/AXISDATE M Year NoMonth NoDay /* * * * * * * * * * /* Plot load textsymbol 14;textcolor 1;textscale 1 MOVE %xtxtp% %ytxtp% text [quote %nl% loads calculated: Log flow range = [substr %logqMin% 1 4] to [substr %logqMax% 1 4] .] ul &call nextrow &if %nl% <= 0 &then &do &type No loads calculated! &return &end GRAPHLIMITS %xgrl1% %ygrl1% %xgrl2% %ygrl2% GRAPHEXTENT %.startday% 0 %.endday% %loadmax% UNITS GRAPH linesymbol 1; linecolor black AXIS vertical AXISRULER /*LINE %.startday% 0 %.startday% %loadmax% TEXTOFFSET [calc -1 * ( 4 * %xsp% ) ] 0 TEXTANGLE 90 TEXTSCALE 2 /*AXISTEXT kg/day TEXTOFFSET 0 0 TEXTANGLE 0 AXIS horizontal /*LINE %.startday% 0 %.endday% 0 &r /rpd2/spek/prjws8/users/michael/aml/AXISDATE M Year NoMonth NoDay LINESYMBOL 13 LINECOLOR 6 UNITS graph clearselect %tl% info reselect %tl% info LOAD >= 0 GRAPHBAR %tl% INFO NDATE LOAD &sv mlsym = 1 &do nyr = %Year1% &to %Year2% isodate2day %nyr%-01-01 .Day1 isodate2day %nyr%-12-31 .Day2 clearselect %tl% info reselect %tl% info ndate >= %.Day1% AND ndate <= %.Day2% &if [extract 1 [show select %tl% info]] = 0 &then &do move [calc 0.5 * ( %.Day1% + %.Day2% )] [calc 0.5 * %loadmax%] text 'No\data' &end &else &do nld = %.Day1% &to %.Day2% clearselect %tl% info reselect %tl% info ndate = %nld% &sv dc = -9 &if [extract 1 [show select %tl% info]] > 0 &then &sv dc = [show select %tl% info 1 item cdelta] &if %dc% = 0 &then &do &sv yl = [show select %tl% info 1 item load ] MOVE %nld% %yl% &sv nmrk = [calc %nld% - %.Day1% + 1] &sv txsc = [show textscale] textscale 0.5 textangle 90 textcolor %mlsym% TEXT %nmrk% cl &if %Debug% &then &type TEXT %nmrk% cc textscale %txsc% textangle 0 &end &end &sv mlsym = %mlsym% + 1 &if %mlsym% > 6 &then &sv mlsym = 1 &end LINESYMBOL 1 LINECOLOR BLACK TEXTCOLOR BLACK UNITS page GRAPHLIMITS %xgrr1% %ygrr1% %xgrr2% %ygrr2% MOVE %xgrrt% %ypgr1% textoffset 0 [calc %rowsize% / 2] TEXT 'log q' lc textoffset 0 0 MOVE %xpgr1% %ygrrt% TEXT 'log\c' cl &sv .nCases = %nl% clearselect %tl% info reselect %tl% info n > 0 reselect %tl% info logq >= -8.9 and logc >= -8.9 /* missing is -9 but real numbers don't equate exactly reselect %tl% info cdelta = 0 /* only non-filled data &r /rpd2/spek/prjws8/users/michael/aml/linregfil .rsq .Yint .slope LOGQ LOGC %tl% /* reads data from %tl% GRAPHEXTENT %logqMin% %logcMin% %logqMax% %logcMax% &sv graphxscale = [calc ( [extract 3 [show graphextent]] - [extract 1 [show graphextent]] ) / ~ ( [extract 3 [show graphlimits]] - [extract 1 [show graphlimits]] ) ] &sv graphyscale = [calc ( [extract 4 [show graphextent]] - [extract 2 [show graphextent]] ) / ~ ( [extract 4 [show graphlimits]] - [extract 2 [show graphlimits]] ) ] &sv graphaspect = [calc %graphxscale% / %graphyscale% ] UNITS graph LINESYMBOL 1 LINECOLOR BLACK AXIS horizontal AXISRULER AXIS vertical AXISRULER MARKERCOLOR 6 TEXTCOLOR 6 TEXTSYMBOL 14 TEXTSCALE 0.5 &sv mlsym = 1 &do nyr = %Year1% &to %Year2% &sv Date1 = %nyr%-01-01 &sv Date2 = %nyr%-12-31 CLEARSELECT %tl% INFO RESELECT %tl% info DATE >= %Date1% and DATE <= %Date2% and logq >= -8.9 and logc >= -8.9 reselect %tl% info cdelta = 0 /* only non-filled data &sv nflo = [extract 1 [show select %tl% info]] &type Found %nflo% records for %Date1% to %Date2% &if %nflo% >= 1 &then &do infofile %tl% info %tlt% init isodate2day %nyr%-01-01 .Day1 isodate2day %nyr%-12-31 .Day2 &sv nmrk = 0 markersymbol %mlsym% linesymbol %mlsym% &do nld = %.Day1% &to %.Day2% clearselect %tlt% info reselect %tlt% info ndate = %nld% /*reselect %tlt% info logq > -8 /*reselect %tlt% info logc > -8 &if [extract 1 [show select %tlt% info]] > 0 &then &do &sv xq = [show select %tlt% info 1 item logq] &sv yc = [show select %tlt% info 1 item logc] &type Day %nld% - logq=%xq% logc=%yc% MARKER %xq% %yc% MOVE %xq% %yc% textcolor %mlsym% /*&sv nmrk = %nmrk% + 1 &sv nmrk = [calc %nld% - %.Day1% + 1] /*TEXT [calc [unquote %nmrk%] * 1] TEXT %nmrk% textcolor 1 &end &end /*&echo &on clearselect %tlt% info reselect %tlt% info ndate > 0 &r /rpd2/spek/prjws8/users/michael/aml/linregfil .rsq .Yint .slope LOGQ LOGC %tlt% /* reads data from %tlt% &sv Yint1 = [calc ( %.slope% * %logqMin% ) + %.Yint% ] &sv Yint2 = [calc ( %.slope% * %logqMax% ) + %.Yint% ] &type logQmin %logqMin% Yint1 %Yint1% logQmax %logqMax% Yint2 %Yint2% LINE %logqMin% %Yint1% %logqMax% %Yint2% &sv txtpt1 = [show convert graph %logqMin% %Yint1% page] &sv txtpt2 = [show convert graph %logqMax% %Yint2% page] coordinate keyboard textangle * %txtpt1% %txtpt2% &ty graphextent [show graphextent] &ty graphlimits [show graphlimits] &ty graphxscale %graphxscale% graphyscale %graphyscale% graphaspect %graphaspect% &type Slope %.slope% arctan slope [atan %.slope%] Text angle [show textangle] move %logqMax% %Yint2% textcolor %mlsym% text %nyr% cr textangle 0 textcolor 1 /*&echo &off /*&stop debug angle &end &sv mlsym = %mlsym% + 1 &if %mlsym% > 6 &then &sv mlsym = 1 &end /*GRAPHPOINT LOAD.TMP INFO LOGQ LOGC TEXTCOLOR 1 TEXTSCALE 1 LINESYMBOL 1 UNITS page TEXTJUSTIFICATION ur TEXTSCALE 2.5 MOVE %xpgl2% %ypgl2% /*TEXT 'kg/yr' &sv slope+1 = [calc %.slope% + 1 ] &sv FlowFactor = [calc ( %nl% * %MeanQ% / %qSum% ) ** %slope+1% ] &sv LoadEst1 = [calc ( %qcSum% / %nl% ) ] &sv LoadEst2 = [calc %MeanQ% * %cSum% / %nl% ] &sv LoadEst3 = [calc %meanQ% * %qcSum% / %qSum% ] &sv LoadEst4 = [calc %LoadEst1% * %FlowFactor% ] &sv LoadEst1r = [round [calc %LoadEst1% / 1000 ]] &sv LoadEst2r = [round [calc %LoadEst2% / 1000 ]] &sv LoadEst3r = [round [calc %LoadEst3% / 1000 ]] &sv LoadEst4r = [round [calc %LoadEst4% / 1000 ]] textsymbol 14;textcolor 1;textscale 1 MOVE %xtxtp% %ytxtp% &sv txtlin [quote Method 1: Sum(q*c)/n = %LoadEst1r% kg %ChemItem%/ day] text %txtlin% ul &sys echo %txtlin% >> %TextFile% &call nextrow MOVE %xtxtp% %ytxtp% &sv txtlin [quote Method 2: (q_avg*c_sum)/n = %LoadEst2r% kg %ChemItem%/ day] text %txtlin% ul &sys echo %txtlin% >> %TextFile% &call nextrow MOVE %xtxtp% %ytxtp% &sv txtlin [quote Method 3: (q_avg*Sum(q*c))/n = %LoadEst3r% kg %ChemItem%/ day] text %txtlin% ul &sys echo %txtlin% >> %TextFile% &call nextrow &sv Yintr = [calc [round [calc 1000 * %.Yint% ]] / 1000] &sv Slopr = [calc [round [calc 1000 * %.slope%]] / 1000] MOVE %xtxtp% %ytxtp% &sv txtlin [quote Regression: logc=(%Slopr%*logq)+(%Yintr%) r^2=%.rsq%] text %txtlin% ul &sys echo %txtlin% >> %TextFile% &call nextrow MOVE %xtxtp% %ytxtp% &sv txtlin [quote Method 4: (Sum(q*c)/n)*(n*q_avg/Sum(q))^(slope+1) = %LoadEst4r% kg %ChemItem%/ day] text %txtlin% ul &sys echo %txtlin% >> %TextFile% &sv txtlin [quote year, n, mean_m3_day, %ChemItem%_kgObs, Lkg_method_4] &sys echo %txtlin% >> %TextFile% &sv newColumn = .TRUE. &call nextrow /*&call TableHead LINESYMBOL 1 LINECOLOR 5 GRAPHLIMITS %xgrl1% %ygrl1% %xgrl2% %ygrl2% GRAPHEXTENT %.startday% 0 %.endday% %loadmax% /*UNITS graph &do nyr = %Year1% &to %Year2% &sv Date1 = %nyr%-01-01 &sv Date2 = %nyr%-12-31 CLEARSELECT %tf% INFO RESELECT %tf% info DATE >= %Date1% and DATE <= %Date2% &sv nflo = [before [show select %tf% info] ,] CLEARSELECT %tl% INFO RESELECT %tl% info DATE >= %Date1% and DATE <= %Date2% and logq >= -8.9 and logc >= -8.9 reselect %tl% info cdelta = 0 /* only non-filled data &sv nload = [before [show select %tl% info] ,] &r /rpd2/spek/prjws8/users/michael/aml/isodate2day %nyr%-01-01 .Day1 &sv text_lm1 = NoData &sv text_lm4 = NoData &if %nflo% >= 1 &then &do STATISTICS %tf% info mean FLOW end &sv MeanQyr = [calc [show statistic 1 1] * 86400] MOVE %xtxtp% %ytxtp% &sv txtlin [quote %nyr%: Mean flow = [round %MeanQyr%] m3 / day] /*text %txtlin% ul &sv txtlin0 %txtlin% &sv FlowFactor = [calc ( %nl% * %MeanQyr% / %qSum% ) ** %slope+1% ] &sv LoadEstyr = [calc %LoadEst1% * %FlowFactor% ] &type Estimated load_4 for %nyr% = [calc %LoadEstyr% * 365 / 1000 ] kg / year TEXTJUSTIFICATION ul /*MOVE %.Day1% %loadmax% /*TEXT [round [calc %LoadEstyr% * 365 / 1000 ]] &if %nload% > 0 &then &do STATISTICS %tl% info sum LOAD sum QXC end &sv SumLyr = [show statistic 1 1] &sv SumQxCyr = [show statistic 2 1] &sv text_lcal = [round %SumLyr%] &end &else &sv text_lcal = NoData &sv text_mqd = [round %MeanQyr%] &sv text_lm1 = [round [calc %SumQxCyr% * 365.24 * 86400 / %nload% / 1000 ]] &sv text_lm4 = [round [calc %LoadEstyr% * 365.24 / 1000 ]] MOVE %xtxtp% %ytxtp% text %nyr% ul MOVE [calc %xtxtp% + ( 0.5 * %wTab% )] %ytxtp% text %nload% ul MOVE [calc %xtxtp% + %wTab%] %ytxtp% text %text_mqd% ul MOVE [calc %xtxtp% + ( 2 * %wTab% )] %ytxtp% text %text_lm1% ul /*text %text_lcal% ul MOVE [calc %xtxtp% + ( 3 * %wTab% )] %ytxtp% text %text_lm4% ul &sv txtlin [quote %nyr%, %nload%, %text_mqd%, %text_lcal%, %text_lm4%] /*text %txtlin% ul /*&sv txtlin [quote [unquote %txtlin0% %txtlin%]] &sys echo %txtlin% >> %TextFile% &call nextrow &end &else &do MOVE %xtxtp% %ytxtp% text %nyr% ul MOVE [calc %xtxtp% + %wTab%] %ytxtp% TEXT 'No flow' ul &call nextrow &end &end /*&sv delflow = [delete %tf% -INFO] /*&sv delchem = [delete %tc% -INFO] /*&sv delfc = [delete %tl% -INFO] /*&sv delfc = [delete %tlt% -INFO] /*&sv deltfc = [delete %tlx% -file] listoutput screen quit postscript %PlotFile% %PlotFile%.ps &sv PostPlot = /iwqs3/tmp/barcode/michael/in/%PlotFile%.ps &if [exists %PostPlot% -file] &then &sv df [delete %PostPlot% -file] mv %PlotFile%.ps %PostPlot% &type Acrobat Distiller can now convert %PostPlot% /*draw %PlotFile% 9999 3 &sv closeqc = [CLOSE -all] &return /* - - - - - - - - - - - - - - - - - &routine FlowFile /* As used in BARCODE.AML &if [exists FLOWL.TMP -info] &then &sv deli [delete FLOWL.TMP -info] &if [exists %Flow-File% -file] &then &do MOVE %xtxtp% %ytxtp% textsymbol 14;textcolor 1;textscale 1 text [quote Using flow file %Flow-File%] ul &call nextrow &DATA arc info ARC DFMT YMD-/ DEFINE FLOWL.TMP DATE , 8,11,D FLOW , 4, 9,F,3 CODE , 1, 1,C NDATE , 4,10,F,0 ADD FROM %Flow-File% CALC NDATE = DATE Q STOP &end &end &else &type Sorry - no file %Flow-File% &return /* - - - - - - - - - - - - - - - - - &routine CheckDate /* from Barcode.aml /* check and set date limits: &sv MonList = Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec &type Checking %Year1% and %Year2% &if [type %Year1%] > 0 &then /* assume 1984-01-01 format &do &sv ISOdate1 = %Year1% &sv Year1 = [substr %ISOdate1% 1 4] &sv Mon1 = [substr %ISOdate1% 6 2] /*&sv nMon = 0 /*&sv Month1 = Xxx /*&do Mon &list %MonList% /*&sv nMon = %nMon% + 1 /*&if %nMon% = %Mon1% &then &sv Month1 = %Mon% /*&end &select %Mon1% &when 1 &sv Month1 = Jan &when 2 &sv Month1 = Feb &when 3 &sv Month1 = Mar &when 4 &sv Month1 = Apr &when 5 &sv Month1 = May &when 6 &sv Month1 = Jun &when 7 &sv Month1 = Jul &when 8 &sv Month1 = Aug &when 9 &sv Month1 = Sep &when 10 &sv Month1 = Oct &when 11 &sv Month1 = Nov &when 12 &sv Month1 = Dec &end &sv Day1 = [substr %ISOdate1% 9 2] &sv Date1 = %Day1%-%Month1%-%Year1% &end &else /* assume 1984 format &do &sv Date1 = 1-Jan-%Year1% &sv ISOdate1 = %Year1%-01-01 &end &if [type %Year2%] > 0 &then /* assume 1984-01-01 format &do &sv ISOdate2 = %Year2% &sv Year2 = [substr %ISOdate2% 1 4] &sv Mon2 = [substr %ISOdate2% 6 2] &sv nMon = 0 &sv Month2 = Xxx &do Mon &list %MonList% &sv nMon = %nMon% + 1 &if %nMon% = %Mon2% &then &sv Month2 = %Mon% &end &sv Day2 = [substr %ISOdate2% 9 2] &sv Date2 = %Day2%-%Month2%-%Year2% &end &else /* assume 1984 format &do &sv Date2 = 31-Dec-%Year2% &sv ISOdate2 = %Year2%-12-31 &end &if %Year1% < 1000 or %Year2% < 1000 &then &stop Please use 4-digit years! /* Arc8 requires delimiters (/ or -) in dates: &sv DateRange = Date >= %ISOdate1% & date <= %ISOdate2% isodate2day %ISOdate1% .startday isodate2day %ISOdate2% .endday &return /* - - - - - - - - - - - - - - - - - &routine ISOdate /* date in ISO format &sv ISOdate [date -year]-[substr[date -usa] 1 2]-[substr[date -vmsdate] 1 2] - [before[date -vmstime] .] &return /* - - - - - - - - - - - - - - - - - &routine Percentile /* Calculate nth percentile. Algorithm supplied by John Carter 1993. /* Calculate the record containing the percentile... /* Sort the data for each variable and select the nth percentile. /* First calculate the (possibly) theoretical ExactRecord for the percentile, /* then find the two integers surrounding this value: /* Inputs: /* Percentile = the percentile sought, e.g. 95 /* ndata = the number of records /* datafile = the temporary INFO file with the data /* datavar = the variable to calculate the percentile on /* Output: /* PcntVar = the percentile value &sv PcntRec = [round [max 1 [calc %Percentile% / 100 * %ndata%]]] &if %PcntRec% > %ndata% &then &sv PcntRec = %ndata% &sv ExactRec = [calc ( %Percentile% * %ndata% ) / 100 ] &sv Recd1 [truncate %ExactRec%] &if %Recd1% < 1 &then &sv Recd1 = 1 &sv Recd2 = [calc %Recd1% + 1] &if %Recd2% > %ndata% &then &sv Recd2 = %ndata% &data ARC INFO ARC SELECT [translate [entryname %datafile%]] SORT ON [translate %datavar%] Q STOP q &end &sv FirstValue = [show select %datafile% INFO %Recd1% item %datavar%] &sv SecondValue = [show select %datafile% INFO %Recd2% item %datavar%] /* Interpolate between these two records, by the cunning method of finding out /* the exact percentile of the two records surrounding the percentile we want: &sv Val2-1 = [calc %SecondValue% - %FirstValue%] &sv Pcntile1 = [calc ( %Recd1% / %ndata% ) * 100] &sv Pcntile2 = [calc ( %Recd2% / %ndata% ) * 100] &if %Pcntile1% <> %Pcntile2% &then &do &sv PcntRatio = [calc ( %Percentile% - %Pcntile1% ) / ( %Pcntile2% - %Pcntile1% ) ] &sv PcntVar = [calc ( %PcntRatio% * %Val2-1% ) + %FirstValue%] &end &else &sv PcntVar = %FirstValue% &return /* - - - - - - - - - - - - - - - - - &routine nextrow /*%xtxtp% &sv ytxtp = [calc %ytxtp% - %rowsize%] &if ( [calc %ytxtp% - ( 2 * %rowsize% )] < %ytx1% ) or %newColumn% &then &do &sv xtxtp = [calc %xtxtp% + %wCol%] &sv ytxtp = [calc %ytx2% - %rowsize%] &sv newColumn = .FALSE. &call TableHead &end &return /* - - - - - - - - - - - - - - - - - &routine LoadCalc &if [exists %tlx% -file] &then &sv dtlx [delete %tlx% -file] &sv loadMax = 0 &sv qSum = 0 &sv cSum = 0 &sv qcSum = 0 &sv logqMin = 999999999 &sv logqMax = -999999999 &sv logcMin = 999999999 &sv logcMax = -999999999 &sv LoadTmp = [locase [show workspace]/%Station%_load.tmp] &sv Load-File = [OPEN %LoadTmp% fstat -w] MARKERSYMBOL 1 &sv nl = 0 &do nDay = %startday% &to %endday% clearselect %tc% info reselect %tc% info ndate = %nDay% &if [extract 1 [show select %tc% info]] = 1 &then &do &sv cvalue = [show select %tc% info 1 item %ChemItem%] &sv cdelta = 0 &sv cdate = [show select %tc% info 1 item date] &end &else /* do linear gap filling: &do &sv tx = %tc% &sv xItem = %ChemItem% &call Filler &sv cvalue = %xvalue% &sv cdelta = %xdelta% &end clearselect %tf% info reselect %tf% info ndate = %nDay% &if [extract 1 [show select %tf% info]] = 1 &then &do &sv fvalue = [show select %tf% info 1 item flow] &sv fdelta = 0 &sv fdate = [show select %tf% info 1 item date] &end &else /* do linear gap filling: &do &sv tx = %tf% &sv xItem = flow &call Filler &sv fvalue = %xvalue% &sv fdelta = %xdelta% &end &if %Debug% &then &type Fill data: C=%cvalue% (d=%cdelta%) and Q=%fvalue% (d=%fdelta%) /*&if ( ( %fvalue% >= 0 ) AND ( %cvalue% > 0 ) ) &then &if ( ( %fvalue% > 0 ) AND ( %cvalue% > 0 ) ) &then &do /* Calculate loading &sv nl = %nl% + 1 day2isodate %nDay% .ldate &sv fdate = %.ldate% &sv dflow = [calc %fvalue% * 86400] &sv load_n = [calc %dflow% * %cvalue% / 1000 ] &if %load_n% > %loadMax% &then &sv loadMax = %load_n% &sv x%nl% = %nDay% &sv qSum = [calc %dflow% + %qSum% ] &sv cSum = [calc %cvalue% + %cSum% ] &sv q_n = %dflow% &sv c_n = %cvalue% &sv qcSum = [calc %qcSum% + [calc ( %dflow% * %cvalue% ) ]] &if %dflow% > 0 &then &do &sv logq_n = [log10 %dflow% ] &if %logq_n% > %logqMax% &then &sv logqMax = %logq_n% &if %logq_n% < %logqMin% &then &sv logqMin = %logq_n% &end &else &sv logq_n = -9 &if %cvalue% > 0 &then &do &sv logc_n = [log10 %cvalue%] &if %logc_n% > %logcMax% &then &sv logcMax = %logc_n% &if %logc_n% < %logcMin% &then &sv logcMin = %logc_n% &end &else &sv logc_n = -9 &type [format~ '(%1,-4%) %2,10% %3,-12% %4,-8% m3/s %5,-12% m3/day %6,-10% log'~ %nl% %ChemStn% %fdate% %fvalue% [round %dflow%] %logq_n%] &type [format~ '[%1,-10%] = %2,-8% mg/l, Load = %3,-12% kg/day %4,-10% log'~ %ChemItem% %cvalue% %load_n% %logq_n% ] &if %nl% = 1 &then ~ &sys echo n,date,nday,dflow,fvalue,fdelta,cvalue,cdelta,load_n,logq_n,logc_n > %tlx% &sys echo %nl%,%fdate%,%nDay%,%dflow%,%fvalue%,%fdelta%,%cvalue%,%cdelta%,%load_n%,%logq_n%,%logc_n%,-9 >> %tlx% &end &else &type (%nDay%) Invalid Q or C: Q=%fvalue%, C=%cvalue% &end &if [exists %tl% -info] &then &sv dtl [delete %tl% -info] /*&if [exists LOAD.TMP -info] &then &sv deli [delete LOAD.TMP -info] &if [exists %tlx% -file] &then &do &sv datdir = [locase [show workspace]] &data ARC INFO ARC DFMT YMD-/ DEFINE %tl% N ,4, 5,B DATE ,8,11,D NDATE,4,10,F,0 DFLOW,4,12,F,0 FLOW ,4,12,F,3 FDELTA,4,5,B [translate %ChemItem%] ,4,12,F,3 CDELTA,4,5,B LOAD ,4,12,F,3 LOGQ ,4,12,F,3 LOGC ,4,12,F,3 QXC ,4,12,F,3 ADD FROM %datdir%/%tlx% err.tmp RESELECT FLOW GT 0 AND [TRANSLATE %ChemItem%] GT 0 CALC QXC = FLOW * [TRANSLATE %ChemItem%] Q STOP &end &end &else &type stop no load data in %tlx% &return /* - - - - - - - - - - - - - - - - - &routine Filler &if %fill% &then &do &call FindPrev &call FindNext &if %svalue0% = -9 OR %svalue1% = -9 &then &do &sv xvalue = -9 &sv xdelta = -9 &return &end &else &do &sv slope = [calc ( %svalue1% - %svalue0% ) / ( %sday1% - %sday0% )] &sv inter = [calc %svalue0% - ( %slope% * %sday0% )] &sv xvalue = [calc ( %slope% * %nDay% ) + %inter%] day2isodate %nDay% .xvaldate &sv xdate = [subst %.xvaldate% - ''] &sv xdelta = [min %sdelta0% %sdelta1%] &end &end &else &do &sv xvalue = -9 &sv xdelta = -9 &return &end &return /* - - - - - - - - - - - - - - - - - &routine FindPrev /* search for the previous complete record &sv svalue0 = -9 &sv sdate0 = -9 &sv sday0 = -9 &sv sdelta0 = -9 &do sDay = %nDay% &to %startday% &by -1 &sv sdelta0 = [calc %nDay% - %sDay%] &if %Debug% &then &type Finding previous %xItem%: %sdelta0% &message &off clearselect %tx% info reselect %tx% info ndate = %sDay% &message &on &if [extract 1 [show select %tx% info]] > 1 &then &do &type ***** Data excess! Multiple selection for %tx% (using first item) ***** list %tx% info &end &if [extract 1 [show select %tx% info]] >= 1 &then &do &sv svalue0 = [show select %tx% info 1 item %xItem%] &sv sdate0 = [show select %tx% info 1 item date] &sv sday0 = %sDay% &sv sdelta0 = [calc %nDay% - %sDay%] &return &end &end &return /* - - - - - - - - - - - - - - - - - &routine FindNext /* search for the next complete record &sv svalue1 = -9 &sv sdate1 = -9 &sv sday1 = -9 &sv sdelta1 = -9 &do sDay = %nDay% &to %endday% &by +1 &sv sdelta1 = [calc %sDay% - %nDay%] &if %Debug% &then &type Finding next %xItem%: %sdelta1% &message &off clearselect %tx% info reselect %tx% info ndate = %sDay% &message &on &if [extract 1 [show select %tx% info]] > 1 &then &type ***** Data excess! Multiple selection for %tx% (using first item) ***** &if [extract 1 [show select %tx% info]] >= 1 &then &do &sv svalue1 = [show select %tx% info 1 item %xItem%] &sv sdate1 = [show select %tx% info 1 item date] &sv sday1 = %sDay% &sv sdelta1 = [calc %sDay% - %nDay%] &return &end &end &return /* - - - - - - - - - - - - - - - - - &routine TableHead MOVE %xtxtp% %ytxtp% text Year ul MOVE [calc %xtxtp% + ( 0.5 * %wTab% )] %ytxtp% text [quote n Obs] ul MOVE [calc %xtxtp% + %wTab%] %ytxtp% text [quote Mean m3/d] ul MOVE [calc %xtxtp% + ( 2 * %wTab% )] %ytxtp% text [quote %ChemItem% kg Obs] ul MOVE [calc %xtxtp% + ( 3 * %wTab% )] %ytxtp% text [quote Method 4 ] ul &call nextrow &return /* - - - - - - - - - - - - - - - - - &routine deltaPlot /* Shade according to interpolation distance /* set up 100 shaded line symbols in HSV colour space: &sv nSat = 100 &do nSym = 1 &to 100 linesymbol 1 &sv light = [calc 101 - %nSym%] linecolor hsv %hue% %light% %nSat% lineput %nSym% /*&sv nSat = %nSat% - 1 &end STATISTICS %deltaplotfile% INFO min %deltaplotsym% max %deltaplotsym% end &sv deltasym0 = [calc [show statistic 1 1] ] &sv deltasym1 = [calc [show statistic 2 1] ] &do deltasym = %deltasym0% &to %deltasym1% &sv deltaplotcode = [round [calc 1 + ( %deltasym% * ( 100 / ( %deltasym1% + 0.001 ) ) )]] &if %deltaplotcode% > 95 &then &sv deltaplotcode = 95 /* avoid pure white &if %deltaplotcode% < 01 &then &sv deltaplotcode = 01 /* avoid zero symbol clearselect %deltaplotfile% INFO reselect %deltaplotfile% INFO %deltaplotsym% = %deltasym% GRAPHBAR %deltaplotfile% INFO ndate %deltaplotvar% %deltaplotcode% &end lineset plotter;lineset color linesymbol 1; linecolor black; lineput 1 clearselect %deltaplotfile% INFO &return /* - - - - - - - - - - - - - - - - - &routine PageSetup &sv xpg = 42 &sv ypg = 29.7 PAGESIZE %xpg% %ypg% textsize [calc %ypg% / 100] SHADESET colornames markerset mineral markerset color lineset plotter;lineset color &sv rowsize = [calc 1.0 * [extract 1 [show textsize]]] BOX 0 0 %xpg% %ypg% &sv xsp = [calc %xpg% / 84.0] &sv ysp = [calc %ypg% / 59.4] &sv xpg1 = %xsp% &sv ypg1 = %ysp% &sv xpg2 = [calc %xpg% - %xsp%] &sv ypg2 = [calc %ypg% - %ysp%] /* <--------- wtx ----------> /*+----------------------------------------------------------+ /*| Text (t) ^ | /*| htx | /*| v | /*+----------------------------------------------------------| /*| Load (l) | Locality | /*|<-- wgr -->| Map | /*| | (map) | /*| |<-- wmap -->| /*+--------------------------------------| | /*| Conc (c) | | /*| |-------------------| /*| | | /*| | | /*+--------------------------------------+ | /*| Flow (q) ^ | Regression | /*| hgr | Plot (reg) | /*| | | /*| v | | /*+--------------------------------------+-------------------+ &sv wall = %xpg2% - %xpg1% &sv hall = %ypg2% - %ypg1% &sv wgr = [calc ( 0.75 * %wall% ) - ( 3 * %xsp% ) ] &sv hgr = [calc ( 0.25 * %hall% ) - ( 2 * %ysp% ) ] &sv wmap = [calc ( 0.25 * %wall% ) - ( 0 * %xsp% ) ] &sv hmap = [calc ( 0.35 * %hall% ) - ( 3 * %ysp% ) ] &sv hreg = [calc ( 0.40 * %hall% ) - ( 2 * %ysp% ) ] &sv wtx = [calc ( 1.00 * %wall% ) - ( 2 * %xsp% ) ] &sv htx = [calc ( 0.25 * %hall% ) ] &sv nTxCols = 5 /* number of text columns in wtx &sv nTbCols = 4 /* number of table columns in wCol &sv wCol = [calc %wtx% / %nTxCols%] &sv wTab = [calc %wCol% / %nTbCols%] &sv xtx1 = [calc %xpg1% + %xsp% ] &sv xtx2 = [calc %xtx1% + %wtx% ] &sv ytx1 = [calc %ypg1% + ( 3 * %hgr% ) + ( 4 * %ysp% ) ] &sv ytx2 = [calc %ytx1% + %htx% ] &sv xtxtp = [calc %xtx1% + %xsp%] &sv ytxtp = [calc %ytx2% - %rowsize%] &sv xpgq1 = [calc %xpg1% + %xsp% ] &sv ypgq1 = [calc %ypg1% + %ysp% ] &sv xpgq2 = [calc %xpgq1% + %wgr% ] &sv ypgq2 = [calc %ypgq1% + %hgr% ] &sv xgrq1 = [calc %xpgq1% + ( %wgr% / 10 ) ] &sv ygrq1 = [calc %ypgq1% + ( %hgr% / 6 ) ] &sv xgrq2 = %xpgq2% &sv ygrq2 = %ypgq2% &sv xpgc1 = %xpgq1% &sv ypgc1 = [calc %ypgq2% + %ysp% ] &sv xpgc2 = %xpgq2% &sv ypgc2 = [calc %ypgc1% + %hgr% ] &sv xgrc1 = [calc %xpgc1% + ( %wgr% / 10 ) ] &sv ygrc1 = [calc %ypgc1% + ( %hgr% / 8 ) ] &sv xgrc2 = %xpgc2% &sv ygrc2 = %ypgc2% &sv xpgl1 = %xpgq1% &sv ypgl1 = [calc %ypgc2% + %ysp% ] &sv xpgl2 = %xpgq2% &sv ypgl2 = [calc %ypgl1% + %hgr% ] &sv xgrl1 = [calc %xpgl1% + ( %wgr% / 10 ) ] &sv ygrl1 = [calc %ypgl1% + ( %hgr% / 8 ) ] &sv xgrl2 = %xpgl2% &sv ygrl2 = %ypgl2% &sv xpgr1 = [calc %xpg1% + %wgr% + ( 2 * %xsp% ) ] &sv ypgr1 = [calc %ypg1% + %ysp% ] &sv xpgr2 = [calc %xpgr1% + %wmap% ] &sv ypgr2 = [calc %ypgr1% + %hreg% ] &sv xgrr1 = [calc %xpgr1% + ( %wmap% / 8 ) ] &sv ygrr1 = [calc %ypgr1% + ( %hreg% / 8 ) ] &sv xgrr2 = %xpgr2% &sv xgrrt = [calc ( %xgrr1% + %xgrr2% ) / 2] &sv ygrr2 = %ypgr2% &sv ygrrt = [calc ( %ygrr1% + %ygrr2% ) / 2] &sv xpgm1 = %xpgr1% &sv ypgm1 = [calc %ypgr2% + %ysp% ] &sv xpgm2 = %xpgr2% &sv ypgm2 = [calc %ypgm1% + %hmap% ] &return /* - - - - - - - - - - - - - - - - - &routine PlotMap CLEARSELECT %Catchment% POLY MAPLIMITS %xpgm1% %ypgm1% %xpgm2% %ypgm2% MAPPOSITION cen cen RESELECT %Catchment% POLY SECONDARY = [quote %Secondary%] MAPEXTENT POLY %Catchment% CLIPMAPEXTENT OFF RESELECT %Catchment% POLY SECONDARY = [quote %Secondary%] LINECOLOR 5 &do nOrder = 1 &to 7 &sv width = %nOrder% / 140 clearselect %Rivers% arc reselect %Rivers% arc order = %nOrder% linesize %width% arcs %Rivers% &end polygonshade %Lake% 45 LINECOLOR 2 ARCS %Catchment% CLEARSELECT %StnFile% POINT RESELECT %StnFile% POINT SECONDARY = [quote %Secondary%] RESELECT %StnFile% POINT STNTYPE = [quote H] MARKERSYMBOL 431 MARKERSCALE 0.5 POINTS %StnFile% CLEARSELECT %StnFile% POINT RESELECT %StnFile% POINT SECONDARY = [quote %Secondary%] RESELECT %StnFile% POINT STNTYPE = [quote R] MARKERSYMBOL 436 MARKERSCALE 0.5 POINTS %StnFile% CLEARSELECT %Catchment% POLY RESELECT %Catchment% POLY SECONDARY = [quote %Secondary%] NSELECT %Catchment% POLY POLYGONSHADE %Catchment% 17 MARKERSYMBOL 46 MARKERCOLOR 2 CLEARSELECT %StnFile% POINT RESELECT %StnFile% POINT STATION = [quote %ChemStn%] LIST %StnFile% POINT STATION, DESCRIPTION &sv Description = [SHOW SELECT %StnFile% POINT 1 ITEM DESCRIPTION] POINTSPOT CIRCLE %StnFile% .35 105 outline POINTTEXT %StnFile% STATION # lc TEXTJUSTIFICATION ul MOVE [calc %xpgm1% + %xsp%] %ypgm2% TEXTCOLOR 1 TEXTSCALE 3 textmask rectangle 0.01 TEXT [quote %Secondary%] textmask none BOX %xpgr1% %ypgr1% %xpgr2% %ypgr2% MOVE %xpg1% %ypg2% TEXTJUSTIFICATION ul &sv txtlin [quote %ChemStn%: %ChemItem% %Description%... %ISOdate1% to %ISOdate2%\] text %txtlin% TEXTSYMBOL 14 TEXTSCALE 2 &return