' Name: Maucha
'
' Title: Maucha test script to draw symbols showing major dissolved ions
'
' Description: Originally a Pascal program (1989), converted to C (1992) then AML (1993)
' and now ArcView, to plot a Maucha ionic diagram. This is a standalone version
' unlike MauchaGraphicsShape, which is called from WMAstep.
'
' Author: Michael Silberbauer, November 1996 [still not completed in November 1999!]
' ... got it going October 2000 ... more edits January 2001
' Institute for Water Quality Studies, Department of Water and Sanitation
' Private Bag X313 PRETORIA South Africa 0001
' SilberbauerM@dwaf.gov.za
' Please see end of script for detailed comments.
'-----------------------------------------------------------------------------
theView = av.GetActiveDoc
theDisplay = theView.GetDisplay
thePrj = theView.GetProjection
' Get the GraphicList for the View...
theGraphics = theView.GetGraphics
' Return the extent (the enclosing rectangle) of the View:
theMapExtent = theView.GetDisplay.ReturnVisExtent
' Extract the dimensions of the View from the extent:
'+----------------------------+x2m,y2m
'| |
'| |
'| |
'| |ywm
'| |
'| |
'| |
'+----------------------------+
'x1m,y1m xwm
theLowerLeft = theMapExtent.ReturnOrigin
x1m = theLowerLeft.Getx
y1m = theLowerLeft.Gety
xwm = theMapExtent.GetWidth
ywm = theMapExtent.GetHeight
x2m = x1m+xwm
y2m = y1m+ywm
theUpperRight = x2m@y2m
'MsgBox.Info( "x1m"++x1m.SetFormat( "d" ).AsString++" y1m"++y1m.SetFormat( "d" ).AsString++
' " x2m"++x2m.SetFormat( "d" ).AsString++" y2m"++y2m.SetFormat( "d" ).AsString,
'"Display extents" )
y = 1
x = y*(xwm/ywm)
theDiagonal = ((x^2)+(y^2)).sqrt
'Scale = theDiagonal / 100
Scale = (x2m-x1m)/100
' Hard-coded test data
'Scale = 0.008 'geographic
'Scale = 7500 'projected
xSymbol = 0.5
ySymbol = 0.5
Potassium = 0.51
Sodium = 2.6
Calcium = 31.8
Magnesium = 29.9
Chloride = 5.9
Sulphate = 9.4
TotalAlk = 189.6
logH = 7.75
TotalDisSalt = 317
theProject = av.getProject
DocumentList = theProject.GetDocs
TableList = List.Make
for each anItem in DocumentList
if (anItem.is(table)) then
TableList.Add(anItem)
end
end
_InorganicData = MsgBox.ListAsString(TableList,"Select a file with water quality data","DATA TABLES")
if (_InorganicData=nil) then
MsgBox.Info("No tables available in this project","ERROR: NO TABLE")
exit
end
_IDSymbol = MsgBox.YesNo("Label symbols with station code ?","Maucha: ID Symbols",FALSE)
IDtype = "station"
if (_IDSymbol) then
IDtype = MsgBox.ListAsString({"num_nat","number","station"},"National number, full number or station?","Maucha: Choose labelling")
end
theTable = _InorganicData
theVTab = theTable.getVTab
numrecs = theVTab.GetNumRecords
av.ShowStopButton
av.ShowMsg("Starting process...")
for each record in theVTab
'for each record in 1 .. 10
theField = theVTab.FindField("WMA")
wma = theVtab.ReturnValue(theField,record)
theGraphicGroup = GraphicGroup.Make
theField = theVtab.FindField("Latitude")
lat = theVtab.ReturnValue(theField,record)
theField = theVtab.FindField("Longitude")
lon = theVtab.ReturnValue(theField,record)
theField = theVtab.FindField("Latitude")
latitude = theVtab.ReturnValue(theField,record)
theField = theVtab.FindField("Longitude")
longitude = theVtab.ReturnValue(theField,record)
theField = theVtab.FindField("Station")
feat_name = theVtab.ReturnValue(theField,record)
theField = theVtab.FindField("Medtalkali")
theField = theVtab.FindField("Medtalkalinity")
theField = theVtab.FindField("MedTAL")
TotalAlk = theVtab.ReturnValue(theField,record)
theField = theVtab.FindField("Medsodium")
theField = theVtab.FindField("MedNa")
Sodium = theVtab.ReturnValue(theField,record)
theField = theVtab.FindField("Medmagnesi")
theField = theVtab.FindField("Medmagnesium")
theField = theVtab.FindField("MedMg")
Magnesium = theVtab.ReturnValue(theField,record)
theField = theVtab.FindField("Medsulphat")
theField = theVtab.FindField("Medsulphate")
theField = theVtab.FindField("MedSO4")
Sulphate = theVtab.ReturnValue(theField,record)
theField = theVtab.FindField("Medchlorid")
theField = theVtab.FindField("Medchloride")
theField = theVtab.FindField("MedCl")
Chloride = theVtab.ReturnValue(theField,record)
theField = theVtab.FindField("Medpotassi")
theField = theVtab.FindField("Medpotassium")
theField = theVtab.FindField("MedK")
Potassium = theVtab.ReturnValue(theField,record)
theField = theVtab.FindField("Medcalcium")
theField = theVtab.FindField("MedCa")
Calcium = theVtab.ReturnValue(theField,record)
theField = theVtab.FindField("MedTDS")
TotalDisSalt = theVtab.ReturnValue(theField,record)
theField = theVtab.FindField(IDtype)
_station = theVtab.ReturnValue(theField,record).AsString
'xSymbol = longitude
'ySymbol = latitude
geopoint = Point.Make(lon,lat)
prjpoint = geopoint.ReturnProjected(thePrj)
xSymbol = prjpoint.GetX
ySymbol = prjpoint.GetY
' Convert milligrams per litre to equivalents per litre (equivalents =
' molarity times charge on ion) altering the sequence to allow for the
' way in which AML works out its degrees...
Ion2 = Potassium*1.0/39.10
Ion1 = Sodium*1.0/22.99
Ion8 = Calcium*2.0/40.08
Ion7 = Magnesium*2.0/24.31
Ion6 = Sulphate*2.0/96.07
Ion5 = Chloride*1.0/35.45
A = TotalAlk*2.0/100.00
pH = logH
TDS = TotalDisSalt
Ion4 = A
Ion3 = 0
theIonList = {Ion1,Ion2,Ion3,Ion4,Ion5,Ion6,Ion7,Ion8}
' Allocate colours to the various ions:
theIonColourList = List.Make
theIonColourList.Add(Color.GetYellow)
theIonColourList.Add(Color.GetMagenta)
theIonColourList.Add(Color.GetGray)
theIonColourList.Add(Color.GetCyan)
theIonColourList.Add(Color.GetGreen)
theIonColourList.Add(Color.GetBlue)
theIonColourList.Add(Color.GetGray)
theIonColourList.Add(Color.GetRed)
theIonColourList.Add(Color.GetBlack)
theIonNameList = List.Make
theIonNameList.Add("Na+")
theIonNameList.Add("K+")
theIonNameList.Add("CO3=")
theIonNameList.Add("HCO3-")
theIonNameList.Add("Cl-")
theIonNameList.Add("SO4=")
theIonNameList.Add("Mg++")
theIonNameList.Add("Ca++")
'&if [translate %Key%] = KEY &then
'&do
' &sv Posn1 = 'CL'
' &sv Posn2 = 'LL'
' &sv Posn3 = 'LR'
' &sv Posn4 = 'CR'
' &sv Posn5 = 'CR'
' &sv Posn6 = 'UR'
' &sv Posn7 = 'UL'
' &sv Posn8 = 'CL'
' &sv Posn9 = 'LR'
'&end
' If pH <=1 the calling routine is signalling to skip CO2 calculations:
' (in this version the CO2 step is ignored)
'&if %pH% <= 1 &then
'&do
Ion4 = A
Ion3 = 0
'&end
AnionSum = Ion3+Ion4+Ion5+Ion6
CationSum = Ion1+Ion2+Ion7+Ion8
EquivSum = CationSum+AnionSum
message = "Sum cations ="++CationSum.SetFormat("d.dd").asString++"meq"
message = message++"Sum anions ="++AnionSum.SetFormat("d.dd").asString++"meq"
'MsgBox.Info ( message, "Ions" )
' Calculate the dimensions of the Maucha diagram:
pi = Number.GetPi
angrad = pi/180
XCentre = 0.0
YCentre = 0.0
TotalArea = (EquivSum+1).ln
TotalRadius = ((0.125*TotalArea)/(angrad*22.5).sin).sqrt
AnnoRadius = 1.5*TotalRadius
xcirc = (XCentre*Scale)+xSymbol
ycirc = (YCentre*Scale)+ySymbol
RadiusCirc = TotalRadius*Scale
aCircle = Circle.Make(xcirc@ycirc,RadiusCirc)
gCircle = GraphicShape.Make(aCircle)
aSymbol = RasterFill.Make
aSymbol.SetColor(Color.GetWhite)
aSymbol.SetOLColor(Color.GetBlack)
gCircle.SetSymbol(aSymbol)
gCircle.SetObjectTag("Maucha")
theGraphicGroup.Add(gCircle)
theVtxOffList = {1}
for each Angle in 0..7
RayNumber = Angle+1
Equivalents = theIonList.Get(Angle)
IonArea = Equivalents/EquivSum*TotalArea
IonRadius = IonArea/(TotalRadius*((angrad*22.5).sin))
StartAngle = Angle*45.0
IonAngle = StartAngle+22.5
EndAngle = RayNumber*45.0
x1 = XCentre
y1 = YCentre
x2 = TotalRadius*(angrad*StartAngle).cos
y2 = TotalRadius*(angrad*StartAngle).sin
x3 = IonRadius*(angrad*IonAngle).cos
y3 = IonRadius*(angrad*IonAngle).sin
x4 = TotalRadius*(angrad*EndAngle).cos
y4 = TotalRadius*(angrad*EndAngle).sin
x5 = XCentre
y5 = YCentre
theVtxList = {x1,y1,x2,y2,x3,y3,x4,y4,x5,y5}
'
theVtxOffList = List.Make
for each vtx in 0..4
xv = (theVtxList.Get(2*vtx)*Scale)+xSymbol
yv = (theVtxList.Get((2*vtx)+1)*Scale)+ySymbol
theVtxOffList.Add(xv@yv)
end
mPol = polygon.Make({theVtxOffList})
gmPol = GraphicShape.Make(mPol)
aSym = RasterFill.Make
aSym.SetColor(theIonColourList.Get(Angle))
gmPol.SetSymbol(aSym)
gmPol.SetObjectTag("Maucha")
'theGraphics.Add(gmPol)
theGraphicGroup.Add(gmPol)
av.SetStatus((record/numrecs)*100)
nowDate = Date.Now
nowDate.SetFormat("yyyy-MM-dd")
nowTime = Date.Now
nowTime.SetFormat("hhhh:m")
dtStamp = nowDate.asString++" "++nowTime.AsString
av.ShowMsg(dtStamp++"("++record.AsString++")" ++feat_name.asString++message++lat.SetFormat("d.ddddd").asString++lon.SetFormat("d.ddddd").asString)
'&if [translate %Key%] = KEY &then
'&do
' TEXTJUSTIFICATION [unquote [value Posn%RayNumber%]]
' TEXTSCALE 2
' MOVE %.x6% %.y6%
' TEXT [value IoName%RayNumber%]
' TEXTSCALE 1
'&end
end
if (_IDSymbol) then
nTs = 5 '(number of halo text copies)
HaloColour = Color.GetWhite
TextColour = Color.GetBlack
FontSize = 10
pi = Number.GetPi
angrad = pi/180
hAngle = 360/nTs
Tx = xcirc-(RadiusCirc*0.9)
Ty = ycirc
IDText = GraphicText.Make(_station.AsString.Left(6),Tx@Ty)
IDTextSymbol = IDText.ReturnSymbols.Get(0)
IDTextSymbol.SetSize(FontSize)
IDText.GetSymbol.SetColor(TextColour)
HaloDist = RadiusCirc/20
for each nT in 0..nTs ' draw text several times, slightly offset, for "halo"
HaloText = GraphicText.Make(_station.AsString.Left(6),Tx@Ty)
angvalue = hAngle*nT
theOffX = HaloDist*(angrad*angvalue).cos
theOffY = HaloDist*(angrad*angvalue).sin
HaloText.Offset(Point.Make(theOffX,theOffY))
HaloTextSymbol = HaloText.ReturnSymbols.Get(0)
HaloTextSymbol.SetSize(FontSize)
HaloText.GetSymbol.SetColor(HaloColour)
theGraphicGroup.Add(HaloText)
end
theGraphicGroup.Add(IDText)
end
theGraphicGroup.SetObjectTag("Maucha")
theGraphicGroup.SetUniformScaling(true)
theGraphics.Add(theGraphicGroup)
theGraphics.Invalidate
end
av.ClearMsg
av.ClearStatus
'-----------------------------------------------------------------------------
' Comments from the original AML:
' Michael Silberbauer - August to October 1993
' mauchsym.aml takes a set of chemical concentrations in mg/L
' (K+, Na+, Ca++, Mg++, SO4=, Cl-, HCO3-, CO3=) and produces the
' coordinates required to plot them in the format used by
' Maucha (1932) Hydrochemische Metoden in der Limnologie.
' Binnengewasser 12, 173p,
' and modified by
' Broch E S & Yake W (1969) A modification of Maucha's
' ionic diagram to include ionic concentrations.
' Limnology and Oceanography 14, 933-935.
'
' Please refer to
' "Silberbauer M J & King J M (1991) Geographical trends in the
' water chemistry of wetlands in the south-western Cape Province, South
' Africa. Southern African Journal of Aquatic Sciences 17 (1/2) 82 - 88."
' if you wish to cite this programme or use it for any publication.
' Other useful documents are:
' Day J A & King J M (1995) Geographical patterns, and their origins, in the
' dominance of major ions in South African rivers.
' South African Journal of Science 91, 299-306.
' Day J A (1993) The major ion chemistry of some southern African saline systems.
' Hydrobiologia 267, 37-59.
'
' The Maucha diagram is the specific case of a radial plot in which the major
' cations are plotted on the right-hand side of an eight-point star and the
' major anions are plotted on the left-hand side.
'
' The aqueous CO2 should be split between H2CO3, HCO3- and CO3= using the approach of
' Mackereth FJH, Heron J & Talling JF (1978) Water Analysis: Some Revised
' Methods for Limnologists. Freshwater Biological Association Scientific
' Publication No 36. pp39-42. (but this option not implemented)
'
' Broch and Yake's modification was to use the TDS value to scale the whole
' diagram up or down, so that an idea of the relative concentration of various
' samples can be obtained. A logarithmic scale is necessary when widely
' differing concentrations are plotted on the same map.
'-----------------------------------------------------------------------------