##==== HEADER ========================================================================
## 
## Purpose:
## This is a script to plot temperatur data from British Columbia
## Ministry of Forest, Land and Natural Resources Operation
##
## Modification History:
## code generation -- May 2013, Hasi
## seasonal N-factors with fdd and tdd
##
## Autors: 
## Hasi (Andreas Hasler), andreas_hasler@gmx.ch, BC Forest Dept.
## 
## Description:		Plots of SO and nFactor over time at one location.
##
## Dependencies:		Ttools.R (own toolbox for temperature timeseries treatment,
## 				(contained in STATws.R workspace)
##
## I/O, variables and units:
## input:		STATws.RData workspace (generated by procedure "stat_TTSMAT.R")
##				contains toolbox and different aggregates of time series
##				for details see header of "stat_TTSMAT.R"
##				SO_TO_proflist.csv = list of variables in each profile
##		
## output:		plots Tair, Ts, dT and nFactor time series
##
## units:		t is in inported as "mm-dd-yyyy HH:MM:SS AM/PM"
##				t is internaly in R POSIXct
##				temperature is in ?C
##						
##======================================================================================

## clear current workspace to avoid env. dependency
#rm(list=ls())

##======================================================================================
##   control variables
##======================================================================================

site <- "MID"

loc <- "wx"

## set plot mode: "png","pdf" or "scr" (for screen)
pl <- "png"

## define smoothing window (days)
win <- 15

## make bar for winter and summer N-fact
Nfactbar <- T

## define plot window (use 3 full years!!!)
beginyear <- 2007
tbeg <- paste (beginyear,"-01-01",sep="")
tend <- paste (beginyear+2,"-12-31",sep="")

## set limits in degC for plots
ylim1 <- c(-30,20)
ylim2 <- c(-20,30)

##======================================================================================
##   preparation
##======================================================================================


## load workspace with all data and tools 
if (!exists("STAWS")) load("RData/STATws.RData")

## load toolbox Ttools
source("R_source/Ttools.R")

## get list of sites:
sites <- levels(factor(slist$site))
lsl <- length(sites) ## length of short list
siteloc <- paste(site,loc,sep="_")

## getting list with profil definitions for field site
proflist <- read.csv("meta/SO_TO_proflist.csv", as.is=TRUE)
if (dim(proflist)[2]==1) proflist <- read.csv("meta/SO_TO_proflist.csv", as.is=TRUE, sep=";")
proflist <- proflist[proflist$site==site,]

##get variable names and surface tiype and lapse rate corr:
if (exists("vnam")) rm("vnam")
for (i in 1:(dim(proflist)[1])){
    lloc <- unlist(strsplit(proflist$ground[i],split="_",fixed=T))[1]
    if (!is.na(lloc)) if (lloc == loc) vnam <- proflist[i,2:4]
    if (!is.na(lloc)) if (lloc == loc) typ <- proflist$surftype[i]
    if (!is.na(lloc)) if (lloc == loc) lps <- proflist$lapsrc[i]
}

## col for Ts (ussually = 2; if dT then use 3)
if (typ=="r") Ts <- 3 else Ts <- 2


##======================================================================================
##======================================================================================
## get data:
d <- get("DMT_"&site)
SO <- get("SOTO_"&site)[loc]
if (typ=="r") SO <- SO[2,1] else SO <- SO[1,1]
SO <- round(as.numeric(SO),1)
d <- ts.df.runmean(d,win=win,qlim=0.95)
d <- ts.cut(d,posix(tbeg),posix(tend))

t <- d[,1]
Tair <- d[,as.character(vnam[1])] + lps
Ts <- d[,as.character(vnam[Ts])] 
N <- Ts/Tair

## calculate seasonal N-factors
DMT <- data.frame(t=t, Tair=Tair, Ts=Ts)
DMT[is.na(DMT$Tair)|is.na(DMT$Ts),2] <- NA
DMT[is.na(DMT$Tair)|is.na(DMT$Ts),3] <- NA
tdd <- deg.day(DMT)$tdd
fdd <- deg.day(DMT)$fdd
Nth <- tdd$Ts / tdd$Tair
Nfr <- fdd$Ts / fdd$Tair

## erease first and last freezing N factor (2 complete winters)
Nfr <- Nfr[2:3]
if (siteloc == "MID_wx") Nth[1] <- NA

print("mean SO over time series of " & loc)
print(mean(Tair-Ts,na.rm=T))

##======================================================================================
##   ploting
##======================================================================================

## define grafic device size
plfn <- paste("plots/SO_N_fact_",site,loc,sep='')
pldev(x=4,y=4,pld=pl,plfn=plfn)

## make rows and cols on sheet c(rows,cols):
par(mfrow=c(3,1),xpd=FALSE)
## set margins c(bottom, left, top, right)
par(mar = c(1.5, 3, 1.5, 1) + 0.2, mgp=c(2,0.7,0))   # Leave space for z axis
## title size:
par(cex.main=0.9,cex=0.75)

## first plot *************************
## title 
tlab <- "15d running mean temperature"
title <- paste(siteloc,tlab,sep=' ')
## ploting variables
ylab <- "T (°C)"
plot(t,Tair,main=title,ylab=ylab,col="darkblue",ylim=ylim1,type="l")
lines(t,Ts,col="red")
nice.abline(iy=10)

legend("topright",c("Tair","GST"),text.col=c("darkblue","red"),bg="white",cex=0.7)

## second plot *************************
## plots surface offset over time  
dT <- Ts-Tair

## parameters for barplot
l <- length(dT)

## color for barplot
#if (max(!is.na(dT))==0) yaxt="n" else yaxt="s"
col <- dT >= 0
col[col] <- "orange"
col[dT < 0] <- "blue"
ylab <- "GST - Tair (°C)"

barplot(dT,ylim=ylim2*1.08,space=0,col=col,axes=F,border=NA)

##make lineplot
par (new=TRUE)
## title 
tlab <- "15d temperature difference"
title <- paste(siteloc,tlab,"   (SO=",SO,"°C)",sep=' ')
plot(t,dT,main=title,ylab=ylab,ylim=ylim2,type="l")
nice.abline(iy=10)
legend("topright",c("GST > Tair","GST < Tair"),fill=c("orange","blue"),bg="white",cex=0.7)


## third plot *************************
col <- Tair >= 0
col1 <- rgb(1,0.85,0)
col2 <- "lightblue"

if(Nfactbar){
  ylim<-c(-0.1,1.4)*1.08 - 0.05
  barplot(Nth,xlim=c(0.2,2.82),ylim=ylim,width=0.35,space=1.5,col=col1,border=NA,axes=F)
  par (new=TRUE)  
  barplot(Nfr,xlim=c(-0.25,2.6),ylim=ylim,width=0.35,space=1.7,col=col2,border=NA,axes=F)
  par (new=TRUE) 
}
## colors for pointplot
col <- Tair >= 0
col[col] <- "orange"
col[Tair < 0] <- "blue"
## title 
tlab <- "15d and seasonal N-factor"
title <- paste(siteloc,tlab,sep=' ')
#plot(t,N,main=title,xaxs="i",ylim=c(-0.1,1.4),col=col,cex=0.3,pch=19)
plot(t,N,main=title,ylim=c(-0.1,1.4),col=col,cex=0.3,pch=19)
nice.abline(iy=0.2)
abline(h=1)
legend("topright",c("Nth","Nfr"),col=c("orange","blue"),pch=19,bg="white",cex=0.7)


## close ploting dev for png or pdf
if (pl!="scr") dev.off()

##======================================================================================
##   end
##======================================================================================
