Commit 83e9a04b authored by Ronnie Leeper's avatar Ronnie Leeper
Browse files

Fixed processing bug and created script to process sensor biases within the...

Fixed processing bug and created script to process sensor biases within the Nimbus shield at tower D
parent bf767202
# Analyzed tower-D for sensor related baises
# Created: Ronald D. Leeper
# Date: 03.28.2016
# Updated:
# Directory Setup
SystemDir = "/Users/ronnieleeper/Documents/"
ProjectDir = "Research/MMTS-BiasBudget/"
DataDir = "Data/towerD/"
OutputDir = "Analysis/SensorBiases/"
# Read Sensor data
filename = paste(SystemDir,ProjectDir,DataDir,"DailyObs-towerD.csv",sep="")
d = read.table(file=filename, sep=",", header=T, na.strings="NA")
hdrs = colnames(d)
# Discover data columns of interest
for(x in 1:length(hdrs)){
if(hdrs[x] == "mmtsNimMx") mmtsMx = x
if(hdrs[x] == "mmtsNimMn") mmtsMn = x
if(hdrs[x] == "prtNimMx") prtNimMx = x
if(hdrs[x] == "prtNimMn") prtNimMn = x
if(hdrs[x] == "Date") dt = x
if(hdrs[x] == "maxHr") mxHr = x
if(hdrs[x] == "minHr") mnHr = x
if(hdrs[x] == "mxUWind") mxU = x
if(hdrs[x] == "mnUWind") mnU = x
if(hdrs[x] == "mxVWind") mxV = x
if(hdrs[x] == "mnVWind") mnV = x
if(hdrs[x] == "mxRgIn") mxRin = x
if(hdrs[x] == "mnRgIn") mnRin = x
if(hdrs[x] == "mxRgOut") mxRout = x
if(hdrs[x] == "mnRgOut") mnRout = x
if(hdrs[x] == "mxLWin") mxLin = x
if(hdrs[x] == "mnLWin") mnLin = x
if(hdrs[x] == "mxLWout") mxLout = x
if(hdrs[x] == "mnLWout") mnLout = x
}
#Add Year month column
d[,length(d)+1] = paste(substr(d[,dt],6,7),substr(d[,dt],1,4),sep="-")
hdrs[length(d)] = "dateLabel"
colnames(d)<-hdrs
dtlab=length(hdrs)
# Preliminary Figures
plot(y=d[,mmtsMx],x=d[,prtNimMx],pch=16,col="red",
xlab="NIMBUS PRT Max Temp. (C)",
ylab="NIMBUS MMTS Max Temp. (C)",
main="TowerD Nimbus Sensors")
abline(0,1)
plot(y=d[,mmtsMn],x=d[,prtNimMn],pch=16,col="blue",
xlab="NIMBUS PRT Min Temp. (C)",
ylab="NIMBUS MMTS Min Temp. (C)",
main="TowerD Nimbus Sensors")
abline(0,1)
# Timeseries of max
mar.default <- c(5,4,4,2) + 0.1
par(mar = mar.default + c(0, 1, 0, 0))
tseq = seq(from=1, to=length(d[,1]),by=15)
mxV = max(max(d[,mmtsMx],na.rm=T), max(d[,prtNimMx], na.rm=T))*10
mnV = min(min(d[,mmtsMx],na.rm=T), min(d[,prtNimMx], na.rm=T))*10
devisor = max((1*10^(nchar(trunc(abs(mxV)))-1)),(1*10^(nchar(trunc(abs(mnV)))-1)))
mx = ((floor(mxV/devisor)+1)*devisor)/10
if(abs(mx - mxV) < devisor) mx = ((floor(mxV/devisor)+2)*devisor)/10
mn = (floor(mnV/devisor)-1)*devisor/10
plot(d[,prtNimMx],col="red",type='l',axes=F, xaxs="i", yaxs="i",
xlab="Datetime",ylim=c(mn,mx), lwd=3, main="TowerD Nimbus Sensors",
ylab="Daily Maximum Temperature (C)")
axis(1, at=tseq, label=d[tseq,dtlab])
axis(2)
lines(d[,mmtsMx], col="blue", lwd=3)
legend("topleft", legend=c("PRT-Nimbus","MMTS-Nimbus"), ncol=2, fill, lty=1, lwd=3, col=c("red","blue"))
# Timeseries of min
mar.default <- c(5,4,4,2) + 0.1
par(mar = mar.default + c(0, 1, 0, 0))
tseq = seq(from=1, to=length(d[,1]),by=15)
mxV = max(max(d[,mmtsMn],na.rm=T), max(d[,prtNimMn], na.rm=T))*10
mnV = min(min(d[,mmtsMn],na.rm=T), min(d[,prtNimMn], na.rm=T))*10
devisor = max((1*10^(nchar(trunc(abs(mxV)))-1)),(1*10^(nchar(trunc(abs(mnV)))-1)))
mx = ((floor(mxV/devisor)+1)*devisor)/10
if(abs(mx - mxV) < devisor) mx = ((floor(mxV/devisor)+2)*devisor)/10
mn = (floor(mnV/devisor)-1)*devisor/10
plot(d[,prtNimMn],col="red",type='l',axes=F, xaxs="i", yaxs="i",
xlab="Datetime",ylim=c(mn,mx), lwd=3, main="TowerD Nimbus Sensors",
ylab="Daily Minimum Temperature (C)")
axis(1, at=tseq, label=d[tseq,dtlab])
axis(2)
lines(d[,mmtsMn], col="blue", lwd=3)
legend("topleft", legend=c("PRT-Nimbus","MMTS-Nimbus"), ncol=2, fill, lty=1, lwd=3, col=c("red","blue"))
# Maximum Diffs
mar.default <- c(5,4,4,2) + 0.1
par(mar = mar.default + c(0, 1, 0, 0))
tseq = seq(from=1, to=length(d[,1]),by=15)
mxV = max(d[,mmtsMx] - d[,prtNimMx],na.rm=T)*10
mnV = min(d[,mmtsMx] - d[,prtNimMx],na.rm=T)*10
devisor = max((1*10^(nchar(trunc(abs(mxV)))-1)),(1*10^(nchar(trunc(abs(mnV)))-1)))
mx = ((floor(mxV/devisor)+1)*devisor)/10
if(abs(mx - mxV) < devisor) mx = ((floor(mxV/devisor)+2)*devisor)/10
mn = (floor(mnV/devisor)-1)*devisor/10
plot(d[,mmtsMx] - d[,prtNimMx], col="red",type='l',axes=F, xaxs="i", yaxs="i",
xlab="Datetime",ylim=c(-4,4), lwd=3, main="TowerD Nimbus Sensors",
ylab="Daily Maximum Temperature Difference (C)")
axis(1, at=tseq, label=d[tseq,dtlab])
axis(2)
legend("topleft", legend=c("MMTS - PRT"), ncol=1, fill, lty=1, lwd=3, col=c("red"))
abline(0,0)
# Maximum Diffs
mar.default <- c(5,4,4,2) + 0.1
par(mar = mar.default + c(0, 1, 0, 0))
tseq = seq(from=1, to=length(d[,1]),by=15)
mxV = max(d[,mmtsMn] - d[,prtNimMn],na.rm=T)*10
mnV = min(d[,mmtsMn] - d[,prtNimMn],na.rm=T)*10
devisor = max((1*10^(nchar(trunc(abs(mxV)))-1)),(1*10^(nchar(trunc(abs(mnV)))-1)))
mx = ((floor(mxV/devisor)+1)*devisor)/10
if(abs(mx - mxV) < devisor) mx = ((floor(mxV/devisor)+2)*devisor)/10
mn = (floor(mnV/devisor)-1)*devisor/10
plot(d[,mmtsMn] - d[,prtNimMn], col="blue",type='l',axes=F, xaxs="i", yaxs="i",
xlab="Datetime",ylim=c(-4,4), lwd=3, main="TowerD Nimbus Sensors",
ylab="Daily Minimum Temperature Difference (C)")
axis(1, at=tseq, label=d[tseq,dtlab])
axis(2)
legend("topleft", legend=c("MMTS - PRT"), ncol=1, fill, lty=1, lwd=3, col=c("blue"))
abline(0,0)
mean(d[,mmtsMn] - d[,prtNimMn],na.rm=T)
mean(d[,mmtsMx] - d[,prtNimMx],na.rm=T)
......@@ -44,7 +44,7 @@
SystemDir = "/Users/ronnieleeper/Documents/"
ProjectDir = "Research/MMTS-BiasBudget/"
DataDir = "Data/ATDD-TowerData/"
OutputDir = "Data/"
OutputDir = "Data/towerD/"
# Start and End Dates
endDate = as.Date("20130815", format="%Y%m%d") #Several days before the Aug 23rd end date
......@@ -66,7 +66,7 @@ for (tower in towers){
data = read.table(file=filename, header=FALSE, sep=",")
#Tower Header file
filename = paste(SystemDir,ProjectDir,DataDir,"hdrsTower",tower,".csv",sep="")
filename = paste(SystemDir,ProjectDir,DataDir,"headers/hdrsTower",tower,".csv",sep="")
hdrs = read.table(file=filename, sep=",", stringsAsFactors=FALSE)
#Invoke column headers using the header file
......@@ -206,7 +206,7 @@ for (tower in towers){
# Remove duplicate columns
dailyD = dailyD[,-c(which(duplicated(dailyHdrs)==T))]
dailyHdrs = colnames(dailyHdrs)
dailyHdrs = colnames(dailyD)
# Find wind and solar radiation values at time (last if more than one) daily extremes
# Create unique id consisting of date & temperature (a seach column)
......@@ -221,6 +221,12 @@ for (tower in towers){
dailyHdrs[length(dailyD)] = "minSearch"
colnames(dailyD)<-dailyHdrs
# Add time of min and max
dailyD[,length(dailyD)+1] = NA
dailyHdrs[length(dailyD)] = "maxHr"
dailyD[,length(dailyD)+1] = NA
dailyHdrs[length(dailyD)] = "minHr"
# Create columns to store time of max wind and solar conditions
dailyD[,length(dailyD)+1] = NA
dailyHdrs[length(dailyD)] = "mxUWind"
......@@ -247,6 +253,7 @@ for (tower in towers){
dailyD[,length(dailyD)+1] = NA
dailyHdrs[length(dailyD)] = "mnLWout"
colnames(dailyD)<-dailyHdrs
dHdrs = colnames(data)
#Loop per day
for(x in 1:length(dailyD[,1])){
......@@ -257,6 +264,10 @@ for (tower in towers){
mnInd = which(data[,length(data)] == dailyD[x,"minSearch"])
mnInd = mnInd[length(mnInd)]
#Extract time of max and min
dailyD[x,"maxHr"] = data[mxInd,"Hour"]
dailyD[x,"minHr"] = data[mnInd,"Hour"]
#Extract solar and wind observations at these points
#At time of maximum
dailyD[x,"mxRgIn"] = data[mxInd,inShtRmavg] # Incoming Shortwave Radiation
......@@ -277,109 +288,122 @@ for (tower in towers){
}
}
#Relable column titles
dailyHdrs[which(dailyHdrs=="maxNimbusTemp")] = "mmtsNimMx"
dailyHdrs[which(dailyHdrs=="minNimbusTemp")] = "mmtsNimMn"
dailyHdrs[which(dailyHdrs=="maxPRT_Nim (PRT inside of MMTS shelter)")] = "prtNimMx"
dailyHdrs[which(dailyHdrs=="minPRT_Nim (PRT inside of MMTS shelter)")] = "prtNimMn"
dailyHdrs[which(dailyHdrs=="maxPRT1 5minMavg")] = "prt1_MetMvAvgMx"
dailyHdrs[which(dailyHdrs=="maxPRT2 5minMavg")] = "prt2_MetMvAvgMx"
dailyHdrs[which(dailyHdrs=="maxPRT3 5minMavg")] = "prt3_MetMvAvgMx"
dailyHdrs[which(dailyHdrs=="minPRT1 5minMavg")] = "prt1_MetMvAvgMn"
dailyHdrs[which(dailyHdrs=="minPRT2 5minMavg")] = "prt2_MetMvAvgMn"
dailyHdrs[which(dailyHdrs=="minPRT3 5minMavg")] = "prt3_MetMvAvgMn"
colnames(dailyD)<-dailyHdrs
#Output dailyD results
filename = paste(SystemDir,ProjectDir,DataDir,"DailyObs-tower",tower,".csv",sep="")
filename = paste(SystemDir,ProjectDir,OutputDir,"DailyObs-tower",tower,".csv",sep="")
write.table(dailyD, file=filename, row.names=F, na="NA",sep=",")
#TODO: REMOVE
print("Completed Tower processing!!")
## Preliminary anlaysis ##
tDdailyMx[,length(tDdailyMx)+1] = tDdailyMx[,3] - tDdailyMx[,2]
data[,length(data)+1] = abs(data[,12] - data[,8])
daily16s10Ccount = aggregate(data[data[,43] <=10 ,c(mmtsNimbus)], by=list(data[data[,43] <=10,dailyID]), FUN=length)
daily16s5Ccount = aggregate(data[data[,43] <=5 ,c(mmtsNimbus)], by=list(data[data[,43] <=5,dailyID]), FUN=length)
tseq = seq(from=1, to=length(data[,1]), by=5)
plot(daily16s10Ccount[,2]/5400*100,col="red",pch=16, ylab=expression(bold("Percent of Available Nimbus Observations (%)")),
main=expression(bold("Tower D; 10"*degree*"C Threshold")),
xlab=expression(bold("Date")),axes=F, xaxs="i", yaxs="i", ylim=c(0,120))
axis(1, at=tseq, label=daily16s10Ccount[tseq,1])
axis(2)
tseq = seq(from=1, to=length(data[,1]), by=5)
plot(daily16s5Ccount[,2]/5400*100,col="red",pch=16, ylab=expression(bold("Percent of Available Nimbus Observations (%)")),
main=expression(bold("Tower D; 5"*degree*"C Threshold")),
xlab=expression(bold("Date")),axes=F, xaxs="i", yaxs="i", ylim=c(0,120))
axis(1, at=tseq, label=daily16s5Ccount[tseq,1])
axis(2)
hist(tDdailyMx[,7],xlab="Daily Maximum Difference MMTS-Nimbus minus MMTS-PRT", main="Tower D", xlim=c(-5,5), breaks=seq(from=-30,to=30,by=0.25),xaxs="i", yaxs="i")
plot(x=tDdailyMx[,2],y=tDdailyMx[,3],xlab=expression(bold("MMTS-Nimbus Daily Maximum ("*degree*"C)")),ylab=expression(bold("MMTS-PRT Daily Maximum ("*degree*"C)")),main="Tower D", xlim=c(-20,40), ylim=c(-20,40))
abline(1,1)
minV=-70
maxV=120
filename = paste(SystemDir,ProjectDir,OutputDir,"tower",tower,"16Secplot.png",sep="")
tiff(file=filename,width=700,height=700)
plot(x=data[data[,43] <= 5, 8],y=data[data[,43] <= 5,12],xlim=c(minV,maxV), ylim=c(minV,maxV),
xlab=expression(bold("MMTS-Nimbus Temperature ("*degree*"C)")),
ylab=expression(bold("MMTS-PRT Temperature ("*degree*"C)")),
main="Tower D", col="gray",pch=16)
points(x=data[data[,43] > 5 & data[,43] <= 10, 8],y=data[data[,43] > 5 & data[,43] <= 10,12],xlim=c(minV,maxV), ylim=c(minV,maxV),col="lightblue",pch=16)
points(x=data[data[,43] > 10 & data[,43] <= 15, 8],y=data[data[,43] > 10 & data[,43] <= 15,12],xlim=c(minV,maxV), ylim=c(minV,maxV),col="green",pch=16)
points(x=data[data[,43] > 15 & data[,43] <= 20, 8],y=data[data[,43] > 15 & data[,43] <= 20,12],xlim=c(minV,maxV), ylim=c(minV,maxV),col="yellow",pch=16)
points(x=data[data[,43] > 20 & data[,43] <= 40, 8],y=data[data[,43] > 20 & data[,43] <= 40,12],xlim=c(minV,maxV), ylim=c(minV,maxV),col="orange",pch=16)
points(x=data[data[,43] > 40 , 8],y=data[data[,43] > 40 ,12],xlim=c(minV,maxV), ylim=c(minV,maxV),col="red",pch=16)
legend("topleft", legend=c("<=5","10","15","20","40",">40"),ncol=1,pch=16,col=c("gray","lightblue","green","yellow","orange","red"),title=expression(bold("Difference Magnitude "*degree*"C")))
abline(1,1)
dev.off()
#par(mar = mar.default + c(0, 1, 0, 0))
tseq = seq(from=1, to=length(tDdailyMx[,1]), by=5)
plot(tDdailyMx[,3], type='l', col="blue", axes=F, lwd=3, xaxs="i", yaxs="i", xlab=expression(bold("Date (YYYYMM)")),
ylab=expression(bold("Daily Maximum Temperature ("*degree*"C)")), main="Tower D", ylim=c(-80,150))
axis(1, at=tseq, label=tDdailyMx[tseq,1])
axis(2)
lines(tDdailyMx[,2], col="red", lwd=3)
legend("bottomleft", legend=c("MMTS-PRT","MMTS-Nimbus"), ncol=2, fill, lty=1, lwd=3, col=c("blue","red"))
#16s frequecy timeseries
#All Obs.
filename = paste(SystemDir,ProjectDir,OutputDir,"tower",tower,"16SecTSplotAllObs.png",sep="")
tiff(file=filename,width=700,height=700)
tseq = seq(from=1, to=length(data[,1]), by=3000)
plot(data[,8], ylim=c(-100,150), ylab=expression(bold("MMTS Nimbus Temperature ("*degree*"C)")), main="Tower D; All Obs",
xlab=expression(bold("Date")),col="red", lwd=3, xaxs="i", yaxs="i", axes=F)
lines(data[,12], col="blue", lwd=3)
axis(1, at=tseq, label=data[tseq,31])
axis(2)
legend("topright", legend=c("Nimbus","PRT"), ncol=2, fill, lwd=3, col=c("red","blue"))
dev.off()
#Obs within a 40C check
filename = paste(SystemDir,ProjectDir,OutputDir,"tower",tower,"16SecTSplot40C.png",sep="")
tiff(file=filename,width=700,height=700)
tseq = seq(from=1, to=length(data[,1]), by=3000)
plot(data[data[,43] <= 40,8], ylim=c(-100,150), ylab=expression(bold("MMTS Nimbus Temperature ("*degree*"C)")),
xlab=expression(bold("Date")), main="Tower D; Obs within 40C",
col="red", lwd=3, xaxs="i", yaxs="i", axes=F)
lines(data[data[,43] <= 40 ,12], col="blue", lwd=3)
axis(1, at=tseq, label=data[tseq,31])
axis(2)
legend("topright", legend=c("Nimbus","PRT"), ncol=2, fill, lwd=3, col=c("red","blue"))
dev.off()
#Obs within a 10C check
filename = paste(SystemDir,ProjectDir,OutputDir,"tower",tower,"16SecTSplot10C.png",sep="")
tiff(file=filename,width=700,height=700)
tseq = seq(from=1, to=length(data[,1]), by=3000)
plot(data[data[,43] <= 10,8], ylim=c(-50,50), ylab=expression(bold("MMTS Nimbus Temperature ("*degree*"C)")),
xlab=expression(bold("Date")), main="Tower D; Obs within 10C",
col="red", lwd=3, xaxs="i", yaxs="i", axes=F)
lines(data[data[,43] <= 10 ,12], col="blue", lwd=3)
axis(1, at=tseq, label=data[tseq,31])
axis(2)
legend("topright", legend=c("Nimbus","PRT"), ncol=2, fill, lwd=3, col=c("red","blue"))
dev.off()
#Obs within a 5C check
filename = paste(SystemDir,ProjectDir,OutputDir,"tower",tower,"16SecTSplot5C.png",sep="")
tiff(file=filename,width=700,height=700)
tseq = seq(from=1, to=length(data[,1]), by=3000)
plot(data[data[,43] <= 5,8], ylim=c(-50,50), ylab=expression(bold("MMTS Nimbus Temperature ("*degree*"C)")),
xlab=expression(bold("Date")), main="Tower D; Obs within 5C",
col="red", lwd=3, xaxs="i", yaxs="i", axes=F)
lines(data[data[,43] <= 5 ,12], col="blue", lwd=3)
axis(1, at=tseq, label=data[tseq,31])
axis(2)
legend("topright", legend=c("Nimbus","PRT"), ncol=2, fill, lwd=3, col=c("red","blue"))
dev.off()
\ No newline at end of file
# tDdailyMx[,length(tDdailyMx)+1] = tDdailyMx[,3] - tDdailyMx[,2]
# data[,length(data)+1] = abs(data[,12] - data[,8])
# daily16s10Ccount = aggregate(data[data[,43] <=10 ,c(mmtsNimbus)], by=list(data[data[,43] <=10,dailyID]), FUN=length)
# daily16s5Ccount = aggregate(data[data[,43] <=5 ,c(mmtsNimbus)], by=list(data[data[,43] <=5,dailyID]), FUN=length)
#
# tseq = seq(from=1, to=length(data[,1]), by=5)
# plot(daily16s10Ccount[,2]/5400*100,col="red",pch=16, ylab=expression(bold("Percent of Available Nimbus Observations (%)")),
# main=expression(bold("Tower D; 10"*degree*"C Threshold")),
# xlab=expression(bold("Date")),axes=F, xaxs="i", yaxs="i", ylim=c(0,120))
# axis(1, at=tseq, label=daily16s10Ccount[tseq,1])
# axis(2)
#
# tseq = seq(from=1, to=length(data[,1]), by=5)
# plot(daily16s5Ccount[,2]/5400*100,col="red",pch=16, ylab=expression(bold("Percent of Available Nimbus Observations (%)")),
# main=expression(bold("Tower D; 5"*degree*"C Threshold")),
# xlab=expression(bold("Date")),axes=F, xaxs="i", yaxs="i", ylim=c(0,120))
# axis(1, at=tseq, label=daily16s5Ccount[tseq,1])
# axis(2)
#
# hist(tDdailyMx[,7],xlab="Daily Maximum Difference MMTS-Nimbus minus MMTS-PRT", main="Tower D", xlim=c(-5,5), breaks=seq(from=-30,to=30,by=0.25),xaxs="i", yaxs="i")
# plot(x=tDdailyMx[,2],y=tDdailyMx[,3],xlab=expression(bold("MMTS-Nimbus Daily Maximum ("*degree*"C)")),ylab=expression(bold("MMTS-PRT Daily Maximum ("*degree*"C)")),main="Tower D", xlim=c(-20,40), ylim=c(-20,40))
# abline(1,1)
# minV=-70
# maxV=120
# filename = paste(SystemDir,ProjectDir,OutputDir,"tower",tower,"16Secplot.png",sep="")
# tiff(file=filename,width=700,height=700)
# plot(x=data[data[,43] <= 5, 8],y=data[data[,43] <= 5,12],xlim=c(minV,maxV), ylim=c(minV,maxV),
# xlab=expression(bold("MMTS-Nimbus Temperature ("*degree*"C)")),
# ylab=expression(bold("MMTS-PRT Temperature ("*degree*"C)")),
# main="Tower D", col="gray",pch=16)
# points(x=data[data[,43] > 5 & data[,43] <= 10, 8],y=data[data[,43] > 5 & data[,43] <= 10,12],xlim=c(minV,maxV), ylim=c(minV,maxV),col="lightblue",pch=16)
# points(x=data[data[,43] > 10 & data[,43] <= 15, 8],y=data[data[,43] > 10 & data[,43] <= 15,12],xlim=c(minV,maxV), ylim=c(minV,maxV),col="green",pch=16)
# points(x=data[data[,43] > 15 & data[,43] <= 20, 8],y=data[data[,43] > 15 & data[,43] <= 20,12],xlim=c(minV,maxV), ylim=c(minV,maxV),col="yellow",pch=16)
# points(x=data[data[,43] > 20 & data[,43] <= 40, 8],y=data[data[,43] > 20 & data[,43] <= 40,12],xlim=c(minV,maxV), ylim=c(minV,maxV),col="orange",pch=16)
# points(x=data[data[,43] > 40 , 8],y=data[data[,43] > 40 ,12],xlim=c(minV,maxV), ylim=c(minV,maxV),col="red",pch=16)
# legend("topleft", legend=c("<=5","10","15","20","40",">40"),ncol=1,pch=16,col=c("gray","lightblue","green","yellow","orange","red"),title=expression(bold("Difference Magnitude "*degree*"C")))
# abline(1,1)
# dev.off()
# #par(mar = mar.default + c(0, 1, 0, 0))
# tseq = seq(from=1, to=length(tDdailyMx[,1]), by=5)
# plot(tDdailyMx[,3], type='l', col="blue", axes=F, lwd=3, xaxs="i", yaxs="i", xlab=expression(bold("Date (YYYYMM)")),
# ylab=expression(bold("Daily Maximum Temperature ("*degree*"C)")), main="Tower D", ylim=c(-80,150))
# axis(1, at=tseq, label=tDdailyMx[tseq,1])
# axis(2)
# lines(tDdailyMx[,2], col="red", lwd=3)
# legend("bottomleft", legend=c("MMTS-PRT","MMTS-Nimbus"), ncol=2, fill, lty=1, lwd=3, col=c("blue","red"))
#
# #16s frequecy timeseries
# #All Obs.
# filename = paste(SystemDir,ProjectDir,OutputDir,"tower",tower,"16SecTSplotAllObs.png",sep="")
# tiff(file=filename,width=700,height=700)
# tseq = seq(from=1, to=length(data[,1]), by=3000)
# plot(data[,8], ylim=c(-100,150), ylab=expression(bold("MMTS Nimbus Temperature ("*degree*"C)")), main="Tower D; All Obs",
# xlab=expression(bold("Date")),col="red", lwd=3, xaxs="i", yaxs="i", axes=F)
# lines(data[,12], col="blue", lwd=3)
# axis(1, at=tseq, label=data[tseq,31])
# axis(2)
# legend("topright", legend=c("Nimbus","PRT"), ncol=2, fill, lwd=3, col=c("red","blue"))
# dev.off()
#
# #Obs within a 40C check
# filename = paste(SystemDir,ProjectDir,OutputDir,"tower",tower,"16SecTSplot40C.png",sep="")
# tiff(file=filename,width=700,height=700)
# tseq = seq(from=1, to=length(data[,1]), by=3000)
# plot(data[data[,43] <= 40,8], ylim=c(-100,150), ylab=expression(bold("MMTS Nimbus Temperature ("*degree*"C)")),
# xlab=expression(bold("Date")), main="Tower D; Obs within 40C",
# col="red", lwd=3, xaxs="i", yaxs="i", axes=F)
# lines(data[data[,43] <= 40 ,12], col="blue", lwd=3)
# axis(1, at=tseq, label=data[tseq,31])
# axis(2)
# legend("topright", legend=c("Nimbus","PRT"), ncol=2, fill, lwd=3, col=c("red","blue"))
# dev.off()
#
# #Obs within a 10C check
# filename = paste(SystemDir,ProjectDir,OutputDir,"tower",tower,"16SecTSplot10C.png",sep="")
# tiff(file=filename,width=700,height=700)
# tseq = seq(from=1, to=length(data[,1]), by=3000)
# plot(data[data[,43] <= 10,8], ylim=c(-50,50), ylab=expression(bold("MMTS Nimbus Temperature ("*degree*"C)")),
# xlab=expression(bold("Date")), main="Tower D; Obs within 10C",
# col="red", lwd=3, xaxs="i", yaxs="i", axes=F)
# lines(data[data[,43] <= 10 ,12], col="blue", lwd=3)
# axis(1, at=tseq, label=data[tseq,31])
# axis(2)
# legend("topright", legend=c("Nimbus","PRT"), ncol=2, fill, lwd=3, col=c("red","blue"))
# dev.off()
#
# #Obs within a 5C check
# filename = paste(SystemDir,ProjectDir,OutputDir,"tower",tower,"16SecTSplot5C.png",sep="")
# tiff(file=filename,width=700,height=700)
# tseq = seq(from=1, to=length(data[,1]), by=3000)
# plot(data[data[,43] <= 5,8], ylim=c(-50,50), ylab=expression(bold("MMTS Nimbus Temperature ("*degree*"C)")),
# xlab=expression(bold("Date")), main="Tower D; Obs within 5C",
# col="red", lwd=3, xaxs="i", yaxs="i", axes=F)
# lines(data[data[,43] <= 5 ,12], col="blue", lwd=3)
# axis(1, at=tseq, label=data[tseq,31])
# axis(2)
# legend("topright", legend=c("Nimbus","PRT"), ncol=2, fill, lwd=3, col=c("red","blue"))
# dev.off()
\ No newline at end of file
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment