Vera's R code:
Scripts Directory:
General:
NOTE: this is
simply a copy of the input commands without formatting, comments, or
hierarchy :-).
Starting display moniter:
xhost + stratus
/sbin/ifconfig (to get the ip adress of work station)
export DISPLAY=anzak.hpc.unbc.ca:0
export DISPLAY=142.207.117.35:0.0
OR
xhost + status
ssh -X stratus
mat<-(matrix(1:2,1,2))
layout(mat)
layout.show(2)
mat<-(matrix(1:4,2,2))
layout(mat)
layout.show(4)
mat<-(matrix(1:6,2,3))
layout(mat)
layout.show(6)
hotd<-read.table(file="hotday3.out",header=F,na.strings="-999",col.name=c("dat","jd","yr","mo","dy","n","tmp","tmp3","ppt","spd","dir"))
hotJA<-read.table(file="hotday3JA.out",na.strings="-999",header=F,col.name=c("dat","jd","yr","mo","dy","n","tmp","tmp3","ppt","spd","dir"))
peak<-read.table(file="peak.out",header=F,na.strings="-999",col.name=c("dat","jd","yr","mo","dy","n","tmp","tmp3","ppt","spd","dir"))
peakJA<-read.table(file="peakJA.out",header=F,na.strings="-999",col.name=c("dat","jd","yr","mo","dy","n","tmp","tmp3","ppt","spd","dir"))
brja<-read.table(file="in.dat.HC2",header=F,col.name=c("dat","yr","mo","dy","tmp3","tmp","uk1","uk2","uk3","day3"),na.string="-999",skip=1)
plot(peak$mo,peak$n,xlab="Month",ylab="Cycle Length in Days",main="PE
For All Months",xlim=c(4,10),ylim=c(0,50))
text(6,48, as.expression(substitute(max==12)))
text(6,46,as.expression(substitute(min==4)))
text(6,44,as.expression(substitute(mean==5.2)))
text(6,42,as.expression(substitute(cycles==92)))
dev.copy2eps(file="peak.eps")
plot(hotd$mo,hotd$n,xlab="Month",ylab="Cycle Lenght in Days",main="FW
For All Months",xlim=c(4,10))
text(5,48, as.expression(substitute(max==48)))
text(5,46,as.expression(substitute(min==4)))
text(5,44,as.expression(substitute(mean==8.1)))
text(5,42,as.expression(substitute(cycles==388)))
dev.copy2eps(file="Fair_Weather.eps")
#dev.copy2eps(file="peak_and_FW.eps")
hist(peak$n,breaks=c(4,5,6,7,8,9,10,11,12,13,14,15,16,17,18),main="PE
Cycle Lenghts for All Months",xlab="length of Cycle in
Days",ylim=c(0,70))
text(14,60,as.expression(substitute(mean==5.2)))
#dev.copy2eps(file="CycleFreqPeakAll.eps")
hist(peakJA$n,breaks=c(4,5,6,7,8,9,10,11,12,13,14,15,16,17,18),main="PE
Cycle Lenghts for Jul. & Aug.",xlab="Length of Cycle in
Days",ylim=c(0,70))
text(14,60,as.expression(substitute(mean==5.3)))
#dev.copy2eps(file="CycleFreqPeakJA.eps")
hist(hotd$n,breaks=seq(from=0,to=50,by=1),main="FW Cycle Lenghts for
All Months",xlab="Length of Cycle in Days",freq=T)
text(30,80,as.expression(substitute(mean==8.1)))
#dev.copy2eps(file="CycleFreqFWAll.eps")
hist(hotJA$n,breaks=seq(from=0,to=50,by=1),main="FW Cycle Lenghts for
Jul. & Aug.",xlab="Length of Cycle in Days",freq=T,ylim=c(0,100))
text(30,80,as.expression(substitute(mean==9.1)))
#dev.copy2eps(file="CycleFreqFWJA.eps")
dev.copy2eps(file="histograms.eps")
mat<-(matrix(1:2,1,2))
layout(mat)
layout.show(2)
hist(peak$mo,breaks=c(3,4,5,6,7,8,9,10),ylab="Number of
cycles",bg="blue",main="PE Cycles per
Month",xlab="Month",ylim=c(0,130),axes=F)
axis(2)
axis(1,labels=c(3,4,5,6,7,8,9,10),adj="0")
text(6.5,30,"July")
text(7.5,30,"August")
text(6.5,25,"30%")
text(7.5,25,"25%")
dev.copy2eps(file="MonthFreqPeak.eps")
hist(hotd$mo,main="FW Cycles per Month",ylab="Number of
Cycles",xlab="Month",freq=T,breaks=c(3,4,5,6,7,8,9,10),ylim=c(0,130),axes=FALSE)
axis(2)
axis(1,at=c(3,3.5,4,4.5,5,5.5,6,6.5,7,7.5,8,8.5,9,9.5,10,10.5),pch=c(3,3.5,4,4.5,5,5.5,6,6.5,7,7.5,8,8.5,9,9.5,10,10.5),labels=c(NA,4,NA,5,NA,6,NA,7,NA,8,NA,9,NA,10,NA,NA))
text(6.5,80,"July")
text(7.5,80,"August")
text(6.5,75,"33%")
text(7.5,75,"23%")
#dev.copy2eps(file="MonthFreqFW2.eps")
#dev.copy2eps(file="MonthFreqBoth.eps")
boxplot(split(hotd$n,hotd$mo),main="FW Cycle Length by Month",
xlab="Month", ylab="Length of Cycle in Days")
#dev.copy2eps(file="BPCycleByMonthFW.eps")
boxplot(split(peak$n,peak$mo),main="PE Cycle Length by Month",
xlab="Month", ylab="Length of Cycle in Days",ylim=c(0,50))
#dev.copy2eps(file="BPCycleByMonthPeak.eps")
boxplot(split(peak$yr,peak$mo),main="Annual Distribution of PE Cycles
by Month", xlab="Month")
#dev.copy2eps(file="BPAnnualPeak.eps")
boxplot(split(hotd$yr,hotd$mo),main="Annual Distribution of FW Cycles
by Month", xlab="Month")
#dev.copy2eps(file="BPAnualFW.eps")
dev.copy2eps(file="BoxPlots.eps")
Climatology comparison:
all<-read.table(file="date.out",header=F,na.strings="-999",col.name=c("dat","jd","yr","mo","dy","tmp","tmp3","ppt","spd","dir"))
clm<-read.table(file="dateJA.out",na.strings="-999",header=F,col.name=c("dat","jd","yr","mo","dy","tmp3","tmp","ppt","spd","dir"))
yr,mo,dy,hr,tmp,tmp2,prs,spd,dir,sol
ALL THE FOLLOWING HLY DATA IS FOR JULY AND AUGUST ONLY:
hly<-read.table(file="hlyJA.out",na.strings="-999",header=F,col.name=c("yr","mo","dy","hr","tmp","tmp2","prs","spd","dir","sol"))
hlyp<-read.table(file="hlypeakJA.out",na.strings="-999",header=F,col.name=c("yr","mo","dy","hr","tmp","tmp2","prs","spd","dir","sol"))
hlyfw<-read.table(file="hlyhotday3JA.out",na.strings="-999",header=F,col.name=c("yr","mo","dy","hr","tmp","tmp2","prs","spd","dir","sol"))
hly15<-read.table(file="hlyJA15.out",na.strings="-999",header=F,col.name=c("yr","mo","dy","hr","tmp","tmp2","prs","spd","dir","sol"))
hlyp15<-read.table(file="hlypeakJA15.out",na.strings="-999",header=F,col.name=c("yr","mo","dy","hr","tmp","tmp2","prs","spd","dir","sol"))
hlyfw15<-read.table(file="hlyhotday3JA15.out",na.strings="-999",header=F,col.name=c("yr","mo","dy","hr","tmp","tmp2","prs","spd","dir","sol"))
wind<-read.table(file="hourly/windJA.out",na.strings="-999",header=F,col.name=c("spd","dir"))
plot(peakJA$jd,peakJA$n,main="Jul15-Aug.15 vs all Jul-Aug PE Cycle
Friquency",ylab="Cycle Length in Days",xlab="Julian Date",sub=("Blue
lines deliniate Jul.15 to Aug.15 and include 71% of cycles and 1/2 total
time."))
abline(v=227,col="blue",untf=F)
abline(v=196,untf=F,col="blue")
text(188,11,"Cycles=11")
text(235,11,"Cycles=9")
text(205,11,"Cycles=51")
text(205,11.5,"Frequency=1.7")
text(188,11.5,"Friquency=0.7")
text(235,11.5,"Friquency=0.6")
#dev.copy2eps(file="peakJAvsJA15.eps")
Loop for counting Months:
i<-1;no<-0
while(i<=length(hotd$mo)){
if(hotd$mo[i]==8){
no<-no+1
};
i<-i+1}
mat<-(matrix(1:4,2,2))
layout(mat)
layout.show(4)
hist(clm$dir,breaks=35,main="Wind Direction",xlab="wind
direction",xlim=c(1,36),ylim=c(0,200))
hist(peakJA$dir,add=T,border="green",breaks=35)
dev.copy2eps(file="WindDirClmVsPeak.eps")
hist(clm$spd,breaks=120,main="Extreem Ghust Speed",xlab="Wind
Speed",ylim=c(0,110))
hist(peakJA$spd,breaks=120,add=T,border="green")
dev.copy2eps(file="WindSpdClmVsPeak.eps")
hist(clm$ppt,breaks=seq(from=0,to=500,by=5),freq=T,main="Precipitation",xlab="Precipitation
in 0.1mm",ylim=c(0,130))
hist(peakJA$ppt,breaks=seq(from=0,to=500,by=5),freq=T,add=T,border="green")
dev.copy2eps(file="PptClmVsPeak.eps")
hist(peakJA$ppt,breaks=seq(from=0,to=10,by=1),ylim=c(0,10),freq=T,main="Peak
Precipitation",xlab="Precipitation in 0.1mm")
dev.copy2eps(file="PptPeak.eps")
hist(hly$prs)
hist(hlyp$prs,add=T,border="green")
boxplot(split(hly$prs,hly$hr),main="Pressure",ylab="pressure in 0.01
kPa",xlab="Hour",sub="Means: Peak Emergance=9368,Climatology=9358")
boxplot(split(hlyp$prs,hlyp$hr),add=T,border="blue")
boxplot(split(hlyfw$prs,hlyfw$hr),add=T,border="purple")
#dev.copy2eps(file="PressureClmVsPeak.eps")
boxplot(split(hly$spd,hly$hr),main="Wind Speed",ylab="Speed in
kph",xlab="Hour",sub="Means: Peak Emergance=,Climatology=")
boxplot(split(hlyp$spd,hlyp$hr),add=T,border="blue")
#dev.copy2eps(file="PressureClmVsPeak.eps")
peakJA
boxplot(split(hly$dir,hly$hr),main="Wind
Direction",xlab="Hour",sub="Means: Peak Emergance=,Climatology=")
boxplot(split(hlyp$spd,hlyp$hr),add=T,border="blue")
#dev.copy2eps(file="PressureClmVsPeak.eps")
boxplot(split(hly$tmp,hly$hr),plot=F)
boxplot(split(hlyp$prs,hly$hr),plot=F)
boxplot(split(hlyp$tmp,hly$hr),plot=F)
boxplot(split(hly$tmp,hly$hr),plot=F)
value<-split(hly$prs,hly$hr)
hlym<-sapply(value,mean,na.rm=T)
value<-split(hlyp$prs,hlyp$hr)
hlypm<-sapply(value,mean,na.rm=T)
value<-split(hlyfw$prs,hlyfw$hr)
hlyfwm<-sapply(value,mean,na.rm=T)
value<-split(hly$tmp,hly$hr)
hlymt<-sapply(value,mean,na.rm=T)
value<-split(hlyp$tmp,hlyp$hr)
hlypmt<-sapply(value,mean,na.rm=T)
value<-split(hlyfw$tmp,hlyfw$hr)
hlyfwmt<-sapply(value,mean,na.rm=T)
LINE PLOT
plot(hlym,pch="-",ylim=c(9300,9410),type="l",xlab="Hours",ylab="Pressure",main="Diurnal
Mean Temperature and Pressure",axes=F)
axis(2)
axis(1,c(1,24:1))
add=t
points(hlypm,pch="-",col="blue",type="l")
points(hlyfwm,pch="-",col="purple",type="l")
par(new=T)
plot(hlymt,pch="-",col="orange",ylim=c(90,300),ylab="",xlab="",type="l",axes=F)
axis(side=4)
add=T
points(hlypmt,pch="-",col="red",type="l")
points(hlyfwmt,pch="-",col="dark red",type="l")
text(2,300,"temp.all",col="orange")
text(2,290,"temp.peak",col="red")
text(2,280,"tem.Fair Weather",col="dark red")
text(2,270,"pressure.all")
text(2,260,"pressure.peak",col="blue")
text(2,250,"pressure.Fair Weahter",col="purple")
#dev.copy2eps(file="PrsTmpMeanAll.eps")
mtext(4,Temperature,srt=90,"Temperature")
Point PLOT
plot(hlym,pch=15,ylim=c(9300,9410),xlab="Hours",ylab="Pressure",main="Diurnal
Mean Temperature and Pressure",axes=F)
axis(2)
axis(1,c(1,24:1))
add=t
points(hlypm,pch=17,type="p")
points(hlyfwm,pch=20)
par(new=T)
plot(hlymt,pch=22,ylim=c(90,300),ylab="",xlab="",axes=F)
axis(side=4)
add=T
points(hlypmt,pch=24)
points(hlyfwmt,pch=21)
text(2,300,"tem.Fair Weather=")
text(4,300,pch=15)
text(2,290,"temp.all",col="orange")
text(2,280,"temp.peak",col="red")
text(2,270,"pressure.all")
text(2,260,"pressure.peak",col="blue")
text(2,250,"pressure.Fair Weahter",col="purple")
#dev.copy2eps(file="PrsTmpMeanAll.eps")
ADDING JULY15 TO AUG15
value<-split(hly15$prs,hly15$hr)
hlym15<-sapply(value,mean,na.rm=T)
value<-split(hlyp15$prs,hlyp15$hr)
hlypm15<-sapply(value,mean,na.rm=T)
value<-split(hlyfw15$prs,hlyfw15$hr)
hlyfwm15<-sapply(value,mean,na.rm=T)
value<-split(hly15$tmp,hly15$hr)
hlymt15<-sapply(value,mean,na.rm=T)
value<-split(hlyp15$tmp,hlyp15$hr)
hlypmt15<-sapply(value,mean,na.rm=T)
value<-split(hlyfw15$tmp,hlyfw15$hr)
hlyfwmt15<-sapply(value,mean,na.rm=T)
plot(hlym,pch="-",ylim=c(9300,9410),type="l",xlab="Hours",ylab="Pressure",sub="dashed
lines = Jul15 to Aug15 data only",main="Diurnal Mean Temperature and
Pressure",axes=F)
axis(2)
axis(1,c(1,24:1))
add=t
points(hlypm,pch="-",col="blue",type="l")
points(hlyfwm,pch="-",col="purple",type="l")
par(new=T)
plot(hlymt,pch="-",col="orange",ylim=c(90,300),ylab="",xlab="",type="l",axes=F)
axis(side=4)
add=T
points(hlypmt,pch="-",col="red",type="l")
points(hlyfwmt,pch="-",col="dark red",type="l")
text(2,300,"temp.all",col="orange")
text(2,290,"temp.peak",col="red")
text(2,280,"tem.Fair Weather",col="dark red")
text(2,270,"pressure.all")
text(2,260,"pressure.peak",col="blue")
text(2,250,"pressure.Fair Weahter",col="purple")
#dev.copy2eps(file="PrsTmpMeanAll.eps")
par(new=T)
plot(hlym15,pch="-",ylim=c(9300,9410),type="l",lty=2,xlab="Hours",ylab="Pressure",main="Diurnal
Mean Temperature and Pressure",axes=F)
axis(2)
axis(1,c(1,24:1))
add=t
points(hlypm15,pch="-",col="blue",type="l",lty=2)
points(hlyfwm15,pch="-",col="purple",type="l",lty=2)
par(new=T)
plot(hlymt15,pch="-",col="orange",ylim=c(90,300),ylab="",xlab="",type="l",lty=2,axes=F)
axis(side=4)
add=T
points(hlypmt15,pch="-",col="red",type="l",lty=2)
points(hlyfwmt15,pch="-",col="dark red",type="l",lty=2)
#dev.copy2eps(file="MeanPrsAll.eps")
SAME PLOT FOR SOL:
value<-split(hly$sol,hly$hr)
hlyms<-sapply(value,mean,na.rm=T)
value<-split(hlyp$sol,hlyp$hr)
hlypms<-sapply(value,mean,na.rm=T)
value<-split(hlyfw$sol,hlyfw$hr)
hlyfwms<-sapply(value,mean,na.rm=T)
value<-split(hly15$sol,hly15$hr)
hlym15s<-sapply(value,mean,na.rm=T)
value<-split(hlyp15$sol,hlyp15$hr)
hlypm15s<-sapply(value,mean,na.rm=T)
value<-split(hlyfw15$sol,hlyfw15$hr)
hlyfwm15s<-sapply(value,mean,na.rm=T)
LINE PLOT
plot(hlyms,pch="-",ylim=c(0,3000),type="l",xlab="Hours",ylab="Global
Solar Radiation",sub="dashed lines = Jul15 to Aug15 data
only",main="Diurnal Mean Temperature and Solar Radiation",axes=F)
axis(2)
axis(1,c(1,24:1))
add=t
points(hlypms,pch="-",col="blue",type="l")
points(hlyfwms,pch="-",col="purple",type="l")
points(hlypm15s,pch="*",col="blue",type="l",lty=2)
points(hlyfwm15s,pch="*",col="purple",type="l",lty=2)
points(hlym15s,pch="*",type="l",lty=2)
par(new=T)
plot(hlymt,pch="-",col="orange",ylim=c(90,300),ylab="",xlab="",type="l",axes=F)
axis(side=4)
add=T
points(hlypmt,pch="-",col="red",type="l")
points(hlyfwmt,pch="-",col="dark red",type="l")
points(hlypmt15,pch="-",col="red",type="l",lty=2)
points(hlyfwmt15,pch="-",col="dark red",type="l",lty=2)
points(hlymt15,pch="*",type="l",lty=2)
text(2,300,"temp.all",col="orange")
text(2,290,"temp.peak",col="red")
text(2,280,"tem.Fair Weather",col="dark red")
text(2,270,"sol.all")
text(2,260,"sol.peak",col="blue")
text(2,250,"sol.Fair Weahter",col="purple")
#dev.copy2eps(file="MeanSolAll.eps")
SAME PLOT FOR SPD:
value<-split(hly$spd,hly$hr)
hlymspd<-sapply(value,mean,na.rm=T)
value<-split(hlyp$spd,hlyp$hr)
hlypmspd<-sapply(value,mean,na.rm=T)
value<-split(hlyfw$spd,hlyfw$hr)
hlyfwmspd<-sapply(value,mean,na.rm=T)
value<-split(hly15$spd,hly15$hr)
hlym15spd<-sapply(value,mean,na.rm=T)
value<-split(hlyp15$spd,hlyp15$hr)
hlypm15spd<-sapply(value,mean,na.rm=T)
value<-split(hlyfw15$spd,hlyfw15$hr)
hlyfwm15spd<-sapply(value,mean,na.rm=T)
LINE PLOT
plot(hlymspd,pch="-",type="l",xlab="Hours",ylab="Wind Speed
km/hr",sub="dashed lines = Jul15 to Aug15 data only",main="Diurnal Mean
Temperature and Wind Speed",axes=F)
axis(2)
axis(1,c(1,24:1))
add=t
points(hlypmspd,pch="-",col="blue",type="l")
points(hlyfwmspd,pch="-",col="purple",type="l")
points(hlypm15spd,pch="*",col="blue",type="l",lty=2)
points(hlyfwm15spd,pch="*",col="purple",type="l",lty=2)
points(hlym15spd,pch="*",type="l",lty=2)
par(new=T)
plot(hlymt,pch="-",col="orange",ylim=c(90,300),ylab="",xlab="",type="l",axes=F)
axis(side=4)
add=T
points(hlypmt,pch="-",col="red",type="l")
points(hlyfwmt,pch="-",col="dark red",type="l")
points(hlypmt15,pch="-",col="red",type="l",lty=2)
points(hlyfwmt15,pch="-",col="dark red",type="l",lty=2)
points(hlymt15,pch="*",type="l",lty=2)
text(4,300,"temp.all",col="orange")
text(4,290,"temp.peak",col="red")
text(4,280,"tem.Fair Weather",col="dark red")
text(4,270,"spd.all")
text(4,260,"spd.peak",col="blue")
text(4,250,"spd.Fair Weahter",col="purple")
dev.copy2eps(file="MeanSpdAll.eps")
SAME PLOT FOR SPD amd DIR:
value<-split(hly$dir,hly$hr)
hlymdir<-sapply(value,mean,na.rm=T)
value<-split(hlyp$dir,hlyp$hr)
hlypmdir<-sapply(value,mean,na.rm=T)
value<-split(hlyfw$dir,hlyfw$hr)
hlyfwmdir<-sapply(value,mean,na.rm=T)
value<-split(hly15$dir,hly15$hr)
hlym15dir<-sapply(value,mean,na.rm=T)
value<-split(hlyp15$dir,hlyp15$hr)
hlypm15dir<-sapply(value,mean,na.rm=T)
value<-split(hlyfw15$dir,hlyfw15$hr)
hlyfwm15dir<-sapply(value,mean,na.rm=T)
LINE PLOT
plot(hlymspd,pch="-",type="l",xlab="Hours",ylab="Wind Speed
km/hr",sub="dashed lines = Jul15 to Aug15 data only",main="Diurnal Mean
Temperature and Wind Speed",axes=F)
axis(2)
axis(1,c(1,24:1))
add=t
points(hlypmspd,pch="-",col="blue",type="l")
points(hlyfwmspd,pch="-",col="purple",type="l")
points(hlypm15spd,pch="*",col="blue",type="l",lty=2)
points(hlyfwm15spd,pch="*",col="purple",type="l",lty=2)
points(hlym15spd,pch="*",type="l",lty=2)
par(new=T)
plot(hlymdir,pch="-",col="orange",ylab="",xlab="",type="l",axes=F)
axis(side=4)
add=T
points(hlypmdir,pch="-",col="red",type="l")
points(hlyfwmdir,pch="-",col="dark red",type="l")
points(hlypm15dir,pch="-",col="red",type="l",lty=2)
points(hlyfwm15dir,pch="-",col="dark red",type="l",lty=2)
points(hlym15dir,pch="*",type="l",lty=2)
text(4,20,"dir.all",col="orange")
text(4,18,"dir.peak",col="red")
text(4,16,"dir.Fair Weather",col="dark red")
text(4,14,"spd.all")
text(4,28,"spd.peak",col="blue")
text(4,30,"spd.Fair Weahter",col="purple")
#dev.copy2eps(file="MeanSpdAll.eps")