# Maucha plotting on a map, with acknowledgement to Susanne Walford's R trig skills
# Michael Silberbauer 2011-02-03
# Please cite as
# Silberbauer, M and Walford, S (2011) R script to plot Maucha ionic diagrams spatially.

Maucha<-function(scale, xscale, xm, ym, K, Na, Ca, Mg, SO4, Cl, HCO3, CO3, label){
  colours=c("yellow3","purple","black","cyan","green","blue","grey29","red")
  mlabels=c("Na","K","-","TAL","Cl","SO4","Mg","Ca")
  #ion <- c(K, Na, Ca, Mg, SO4, Cl, HCO3, CO3)
  ion <-c (Na, K, CO3, HCO3, Cl, SO4, Mg, Ca)
  #print(ion)
  anions  <- SO4 + Cl + HCO3 + CO3
  cations <- K + Na + Ca + Mg
  total.ions <- cations + anions
  ps <- total.ions * scale
  sectorAngle <- 360 / 16
  A <- log10(ps+1) # for widely-varying TDS, a log transformation helps to keep the symbols visible
  #A <- ps
  #A <- sqrt(ps)
  sector <- sectorAngle * ( pi / 180 )
  R <- ( (A / 16) *2 / sin(sector) ) ^ 0.5
  circleSector <- sector / 4
  # Plot circle
  x.circ<-array(1:64)
  y.circ<-array(1:64)
  for (i in 1:64 ) {
    x.circ[i] <- R * cos(circleSector*i) * xscale
    y.circ[i] <- R * sin(circleSector*i)
  }
  x.circ <- x.circ + xm
  y.circ <- y.circ + ym
  polygon( x.circ, y.circ, col="white", border="grey80", lwd=0.7 )
  # Plot Maucha symbol
  count <- 1
  for (i in 1:8 ) {
    ion.conc <- ion[i]
    #print(paste(i,": ",ion[i])) # debug statement
    Ai <- ion.conc / total.ions * A
    h <- count - 1
    j <- count + 1
    x.c <- R * cos(sector*h) * xscale
    y.c <- R * sin(sector*h)
    x.d <- R * cos(sector*j) * xscale
    y.d <- R * sin(sector*j)
    a <- Ai / ( R * sin(sector) )
    x.b <- a * cos(sector*count) * xscale
    y.b <- a * sin(sector*count)
    x.t <- 1.1 * max(a,R) * cos(sector*count) * xscale
    y.t <- 1.1 * max(a,R) * sin(sector*count)
    x <- c(0, x.d, x.b, x.c, 0) + xm
    y <- c(0, y.d, y.b, y.c, 0) + ym
    polygon( x, y, col=colours[[i]], border=NA )
    x.t<-x.t+xm
    y.t<-y.t+ym
    if(label) print(paste(x.t,y.t,label,mlabels[i]))
    if(label) text(x.t,y.t,mlabels[i],cex=0.3)
    count = count + 2
  }
}

MauchaData<-function(chem){
    md<-NULL
    md$K  <- chem$K_Diss_Water / 39.0983
    md$Na <- chem$Na_Diss_Water / 22.98977
    md$Ca <- 2 * ( chem$Ca_Diss_Water / 40.08 )
    md$Mg <- 2 * ( chem$Mg_Diss_Water / 24.305 )
    md$SO4<- 2 * ( chem$SO4_Diss_Water / ( 32.06 + (4 * 15.9994) ) )
    md$Cl <- chem$Cl_Diss_Water / 35.453
    md$SO4<- 2 * ( chem$SO4_Diss_Water / ( 32.06 + (4 * 15.9994) ) )
    md$TAL<- chem$TAL_Diss_Water / ( 1.0079 + 12.011 + (3 * 15.9994) )
    md$day<-chem$Jday
    md$year<-chem$Year
    md$d<-chem$Date
    md<-data.frame(md)
    md<-na.omit(md)
}

require(RSQLite)
require(maptools)
debugging=FALSE
hydroyear=TRUE # select whether to use calendar years or hydrological years

pris<-"C:/data_large/av/drainage/hca_1_simple.shp"
pri<-readShapeLines(pris, proj4string=CRS("+proj=longlat +datum=WGS84"))
ters<-"C:/data_large/av/drainage/hca_3.shp"
ter<-readShapeLines(ters, proj4string=CRS("+proj=longlat +datum=WGS84"))
rivr<-"C:/data_large/av/drainage/wri_4_7_SimplifyLine1.shp"
riv<-readShapeLines(rivr, proj4string=CRS("+proj=longlat +datum=WGS84"))
rivr3<-"C:/data_large/av/drainage/wri_3_SimplifyLine.shp"
riv3<-readShapeLines(rivr3, proj4string=CRS("+proj=longlat +datum=WGS84"))
rivr1_2<-"C:/data_large/av/drainage/wri_1_2_SimplifyLine.shp"
riv1_2<-readShapeLines(rivr1_2, proj4string=CRS("+proj=longlat +datum=WGS84"))

yr1<-2001
yr2<-2010
#yr1=1973
#yr2=2011
#yr1=2007
#yr2=2007

drv<-dbDriver("SQLite")
con<-dbConnect(drv, dbname="C:/data_large/av/WMS/wmrq2.db")
dbListTables(con)
pwide<-8.27
phigh<-11.69
pdftitle<-paste("Maucha map","plotted by RQS",format(Sys.time(),"%Y-%m-%d %H:%M:%S"),R.version.string)
pdfname<-paste("C:/tmp/","Maucha_map_",yr1,"_",yr2,".pdf",sep="")
if (hydroyear) {
  pdftitle<-paste("Maucha map plotted by RQS using hydrological years (October to September)",format(Sys.time(),"%Y-%m-%d %H:%M:%S"),R.version.string)
  pdfname<-paste("C:/tmp/","Maucha_map_hydroyears",yr1,"_",yr2,".pdf",sep="")
  }
if(!debugging) pdf(file=pdfname, paper="a4", width=pwide, height=phigh, title=pdftitle)

#w<-read.csv("c:/data_large/av/wms/NCMP_sites_2010_11_16.csv") # choose site list
w<-read.csv("C:/data/AV/National_water_quality/GEMS/Sites_TDS_EC_paper.csv")
for (year in yr1:yr2) {
  date1<-paste(year,"-01-01",sep="")
  date2<-paste(year,"-01-01",sep="")
  dates<-year
  if(hydroyear) {
    date1<-paste(year-1,"-10-01",sep="")
    date2<-paste(year,"-09-30",sep="")
    dates<-paste("- Hydrological year ",year-1,"-",year,sep="")
  }
  or<-par(mar=0.1+c(1.2,1.2,1.2,1.2)) # set up layout
  # Main map
  xr<-c(16.5,33)
  yr<-c(-34,-22)
  plot(ter,col="grey85",xlim=xr,ylim=yr,lwd=0.7)
  plot(pri,col="grey70",xlim=xr,ylim=yr,add=TRUE)
  plot(riv3,col="lightskyblue1",xlim=xr,ylim=yr,add=TRUE,lwd=0.5)
  plot(riv,col="lightskyblue",xlim=xr,ylim=yr,add=TRUE)
  rect(par("usr")[1],par("usr")[3],par("usr")[2],par("usr")[4],border="grey")
  for (nAxis in 1:4) {
    axis(nAxis,col="grey",tcl=-0.2,labels=FALSE)
  }
  for (nAxis in c(1,2,4)) {
    axis(nAxis,col="grey",line=-1,cex.axis=0.5,col.axis="grey",las=1,lwd=0,lwd.ticks=0)
  }
  title(main=paste("NCMP Maucha diagrams", dates),cex.main=1,font.main=2)
  title(sub=paste("Data from WMS plotted by RQS on",format(Sys.time(),"%Y-%m-%d %H:%M:%S"),R.version.string),cex.sub=0.5,line=-1,col.sub="grey")
  
  nid=0
  
  for(featid in c(w$ID)) {
  #-#for(featid in c(w$ID[12:16])) { # select data range
  #for(featid in c(w$ID[1])) {
    nid<-nid+1
    q <- paste ("SELECT PointID, Date, K_Diss_Water, Na_Diss_Water, Ca_Diss_Water, Mg_Diss_Water, SO4_Diss_Water, Cl_Diss_Water, TAL_Diss_Water FROM chem WHERE PointID  = \"",featid,"\" AND Date >= \"",date1,"\" AND Date <= \"",date2,"\"", sep="") 
    c<-dbGetQuery(con, q)
    q<-paste("SELECT PointID, RefCode, Name, Region, Latitude, Longitude, FirstDate, LastDate FROM inv2 WHERE PointID = \"",featid,"\"", sep="")
    s<-dbGetQuery(con, q)
    #print (paste(featid,s$RefCode[1],s$Latitude[1],s$Longitude[1],s$Name[1])) # debug statement
    points(s$Longitude[1], s$Latitude[1], pch=16, cex=0.1, col="grey")
    if( nrow(c) > 1 ) {
      c$Year<-as.numeric(substr(c$Date,1,4))
      c$Date<-as.Date(c$Date, "%Y-%m-%d")
      mc<-MauchaData(c)
      if ( nrow(mc) > 1 ) {
        xscale<-(phigh / pwide) * (par("usr")[2]-par("usr")[1])/(par("usr")[4]-par("usr")[3])
        CO3 <- 0
        mscale<-0.005
        text(s$Longitude[1], s$Latitude[1], substr(s$RefCode[1],1,6), col="grey", cex=0.15)
        print(paste("(",nid,") Maucha",s$RefCode[1],featid,s$Longitude[1], s$Latitude[1]))
        Maucha(mscale, xscale, s$Longitude[1], s$Latitude[1], median(mc$K), median(mc$Na), median(mc$Ca), median(mc$Mg), median(mc$SO4), median(mc$Cl), median(mc$TAL), CO3, FALSE)
        }
        else print (paste(featid,"nrows for Maucha =",nrow(mc) ))
      }
      else print(paste("(",nid,") Maucha: *no data for",year,s$RefCode[1],featid,s$Longitude[1], s$Latitude[1]))
    }
}
if(!debugging) dev.off()

# to execute in R:
# source ("c:/data/program/r/Maucha_map.R")