##==== 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 -- April 2013, Hasi
##
## Autors: 
## Hasi (Andreas Hasler), andreas_hasler@gmx.ch, BC Forest Dept.
## 
## Description:		Plots of SO and TO profiles for each field site.
##
## 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_MAT.R" and "stat_SOTO.R")
##				and WXws.RData workspace
##						
## output:		plot profiles
##
## units:		t is in inported as "mm-dd-yyyy HH:MM:SS AM/PM"
##				t is internaly in R POSIXct
##				temperature is in ?C
##				 
## variables:	SOTO_XXX data frames with surface offsets
##				
## 
## units:		t is internaly in R POSIXct
##				T is temperature in ?C
##				RAIN is in mm
##				RH is relative humidity in %
##				SWR is in W / m'2
##				WS in m / s
##				P is in mBar
##				DP is in ?C
##				WSG in m / s
##				WD in ?		
##						
##======================================================================================

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

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

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

## save ANN_PRE file and dataframe?
save <- TRUE

##======================================================================================
##   preparation and import
##======================================================================================

## load workspaces with all data and tools 
if (!exists("STAWS")) load("RData/STATws.RData")
source("R_source/Ttools.R")
load("RData/WXws.RData")

## get list of sites:
sites <- c("HUD","GUN","PIN","NON","MID","TET","POP")
lsl <- length(sites) ## length of short list

## get prec values from map for each site sites:
##		c("HUD","GUN","PIN","NON","MID","TET","POP")
MAP <- c(  1300, 770,  710,  630,  550,  710,  540)
MSP <- c(  430,  540,  495,  395,  360,  475,  345)

## prepare ploting
## define grafic device size
plfn <- "plots/SOvsPRE"
pldev(x=5,y=5,pld=pl,plfn=plfn)
## make rows and cols on sheet c(rows,cols):
par(mfrow=c(1,1),xpd=FALSE)

##======================================================================================
##======================================================================================
##   main loop through sites

for (k in 1:lsl){	
	## get site name	
	site <- sites[k]
	print("ploting "&site)

 	## getting dataframe SOTO of respective sensor assigned to "soto"
	dfn <- "SOTO_" & site 
	soto <- get(dfn)
	cn <- dim(soto)[2] 
	typ <- soto[soto$var=="typ",2:cn]
	offset <- as.numeric(soto[1,2:cn])
	offset[is.na(offset)] <- as.numeric(soto[2,2:cn])[is.na(offset)]
	mso <- mean(offset,na.rm=TRUE)
	
	## get annual prec from map: (in mapprec)
	precip <- MAP[k]

	
	## get number of offsets per site
	nprof <- length(offset)
	

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

	## plot 
	xlim <- c(370,1310)	
	ylim <- c(0,6)
  ylab <- "SO (°C)"
  xlab <- "annual precipitation (mm)"
	
		## point type
	pty <- c(1,19,2,6,17,0,5,15)
	
	## points for different surface types
	pt2 <- as.numeric(typ)
	pt2[typ=="r"] <- 18
	pt2[typ=="s"] <- 21
	pt2[typ=="f"] <- 17
	pt2[typ=="d"] <- 22
	
	## plot
	if(k==1) plot(precip,mso,type="p",xlim=xlim,ylim=ylim,xlab=xlab,ylab=ylab,cex=1.4)
	points(precip,mso,col=1,pch=pty[k],cex=1.4)
	for (i in 1:nprof) points(precip,offset[i],col="grey",pch=pt2[i],cex=0.8)
	
		
	## plot variation arrows:
	#arrows(mx-v, my, mx+v, my, length = 0.03, angle = 90, code = 3, lwd = 1)
	
}#   END of main loop through the sites

##======================================================================================
##======================================================================================
	
legend("topleft",sites,pch=pty[1:lsl],pt.cex=1.3,bty="n")
legend("bottomleft",c("rock","debris","soil","forest"),pch=c(18,22,21,17),cex=1,pt.cex=0.8,bty="n",col="grey",text.col="darkgrey")

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

if (save) {
	
	ANN_PRE <- data.frame(site = sites, map = MAP, msp = MSP)
	
	## remove unused variables
	rm(list=c("precip","dfn","k","i","lsl","site","mso","soto","cn","pl","MSP","MAP","sites"))
	
 	## save all dataframes for later use (for faster data access)
	save(list=ls(), file = "RData/STATws.RData") 
	
	## save ANN_AMP csv file
	write.csv (ANN_PRE,"CSV_pro/ANN_PRE.csv",row.names=FALSE)
	
}


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