Commit 238dabf0 authored by Ronnie Leeper's avatar Ronnie Leeper
Browse files

Extract wind and solar observations at time of Nimbus max and min.

parent 5771f3f2
......@@ -4,28 +4,40 @@
#Date: 11.02.2015
#Update: 11.17.2015
#Purpose:
#Read in tower station (A-D) data and create daily aggregates from the 16 second observations. From the tower data, generate
#datasets for the respective biases:
# Purpose:
# Read in tower station (A-D) data and create daily aggregates from the 16 second observations. From the tower data, generate
# datasets for the respective biases:
# 1). Sensor bias
# 2). Sheilding bias
# 3). Sitting bias (land cover)
# 4). Time-of-observation (TOB) bias
#11.13.2015 Update:
# 11.13.2015 Update:
# a). Included hourly moving averages of wind speed and solar radiation to evaluate atmospheric conditions at the time of observation
# b). Set outlier observations from MMTS, -73.278 and 124.170 C, to NA
#11.17.2015 Update:
# IS THIS REASONABLE??
# 11.17.2015 Update:
# IS THIS REASONABLE?? YES
# Refined the outlier QC checks to set MMTS-Nimbus obs to NA if they differed from MMTS-PRT obs by more than 5 C
# Added preliminary analysis code to evalute the effectiveness of MMTS-NIMBUS QC checks
#11.23.2015 Update:
# 11.23.2015 Update:
# The refined outlier QC check seems remove all the outliers. According to emails, the study ended (for MMTS systems) on Aug 23 2013.
# The datasets should be truncated to that period. Should we have a common start date for all four towers? It could be immedately
# following the relocation of towers from the period of colocation: Nov. 1 2012
#12.04.2015 Update:
# Determined the 5 Degree C threshold check is reasonable on the 16s frequency data
# TODO:
# 1). Evaluate wind and solar radiation at time of daily extremes -- Completed
# 2). Store daily counts along with daily extremes to QC out later
# 3). Set both Nimbus and PRT observations to missing if either are also missing. -- Completed
# 12.10.2015 Update:
# Implenmented a routine to find the index (rowID) of the last daily max and min occurance and extract wind and solar
# moving averages at that specific row. This will give us wind and solar observations (hourly average) at time of maximum and minimum.
# These values were only extracted at time of Nimbus max and min. Also set PRT and Nimbus observations to NA if either were no available.
#This script is archieved via gitlab at this URL: git@k3.cicsnc.org:ronnieleeper/MMTS-BiasBudget.git
#Directory Setup
......@@ -35,8 +47,8 @@ DataDir = "Data/ATDD-TowerData/"
OutputDir = "Data/"
# Start and End Dates
endDate = 20130815 #Several days before the Aug 23rd end date
startDate = 20121101 #Date towers were positioned at locations from A to D.
endDate = as.Date("20130815", format="%Y%m%d") #Several days before the Aug 23rd end date
startDate = as.Date("20121101", format="%Y%m%d") #Date towers were positioned at locations from A to D.
#Tower IDs
towers = c("A","B","C","D")
......@@ -65,12 +77,12 @@ for (tower in towers){
data[data[,"Day"] < 10,"Day"] = paste("0",data[data[,"Day"] < 10,"Day"],sep="")
data[data[,"Month"] < 10,"Month"] = paste("0",data[data[,"Month"] < 10,"Month"],sep="")
#NOTE: This creates a midnight-to-midnight observer (we may want to apply various day definitions such as 7am-to-7am or 5pm-to-5pm))
data[,length(data)+1] = paste(data[,"Year"],data[,"Month"],data[,"Day"],sep="")
data[,length(data)+1] = as.Date(paste(data[,"Year"],data[,"Month"],data[,"Day"],sep=""),format="%Y%m%d")
hdrs[length(hdrs)+1] = "DailyID (12am Observer)"
colnames(data)<-hdrs
#Truncate data by start and end dates
data = data[data[,length(data)] >= startDate & data[,length(data)] <= endDate,]
#Identify the relative position (column number) of interested variables in the dataset
if (tower == "D"){
......@@ -111,6 +123,15 @@ for (tower in towers){
outLierCriteria = 5
data[abs(data[,mmtsNimbus]-data[,mmtsPRT])>outLierCriteria & !is.na(data[,mmtsNimbus]),c(1:6,mmtsNimbus,mmtsPRT)] = NA
#Set PRT and MMTS obs to missing if either are missing
#Set Nimbus to missing if any of the PRTs are missing
data[is.na(data[,mmtsPRT]) | is.na(data[,met1PRT1]) | is.na(data[,met1PRT2]) | is.na(data[,met1PRT3]), mmtsNimbus] = NA
#Set PRTs to missing with Nimbus is NA
data[is.na(data[,mmtsNimbus]),mmtsPRT] = NA #MMTS-PRT
data[is.na(data[,mmtsNimbus]),met1PRT1] = NA #met1PRT1
data[is.na(data[,mmtsNimbus]),met1PRT2] = NA #met1PRT2
data[is.na(data[,mmtsNimbus]),met1PRT3] = NA #met1PRT3
#Calculate Moving Averages
#Calculate moving averages of windspeed and solar radition based on lag-times (60 minutes)
sixtyMinuteMovingWindowCount = 226 #There are 225 16 second periods within a 60 minute window
......@@ -175,24 +196,182 @@ for (tower in towers){
dataHdrs[2:length(dataHdrs)] = paste("min",dataHdrs[2:length(dataHdrs)],sep="")
colnames(tDdailyMn)<-dataHdrs
#TODO; Evaluate wind and solar radiation at time of maximum (loop for each daily extreme)
#Days with all missing data report infinity; set those observations to missing
# Days with all missing data report infinity; set those observations to missing
tDdailyMx[is.infinite(tDdailyMx[,2]),c(2,3)] = NA
tDdailyMn[is.infinite(tDdailyMn[,2]),c(2,3)] = NA
# Merge daily extreme (max and min files)
dailyD = cbind(tDdailyMx,tDdailyMn)
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)
data[,length(data)+1] = paste(data[,dailyID],data[,mmtsNimbus],sep="")
hdrs[length(data)] = "tempSearch"
colnames(data)<-hdrs
# Create max and min search columns for the daily file
dailyD[,length(dailyD)+1] = paste(dailyD[,"Date"],dailyD[,"maxNimbusTemp"],sep="")
dailyHdrs[length(dailyD)] = "maxSearch"
dailyD[,length(dailyD)+1] = paste(dailyD[,"Date"],dailyD[,"minNimbusTemp"],sep="")
dailyHdrs[length(dailyD)] = "minSearch"
colnames(dailyD)<-dailyHdrs
# Create columns to store time of max wind and solar conditions
dailyD[,length(dailyD)+1] = NA
dailyHdrs[length(dailyD)] = "mxUWind"
dailyD[,length(dailyD)+1] = NA
dailyHdrs[length(dailyD)] = "mnUWind"
dailyD[,length(dailyD)+1] = NA
dailyHdrs[length(dailyD)] = "mxVWind"
dailyD[,length(dailyD)+1] = NA
dailyHdrs[length(dailyD)] = "mnVWind"
dailyD[,length(dailyD)+1] = NA
dailyHdrs[length(dailyD)] = "mxRgIn"
dailyD[,length(dailyD)+1] = NA
dailyHdrs[length(dailyD)] = "mnRgIn"
dailyD[,length(dailyD)+1] = NA
dailyHdrs[length(dailyD)] = "mxRgOut"
dailyD[,length(dailyD)+1] = NA
dailyHdrs[length(dailyD)] = "mnRgOut"
dailyD[,length(dailyD)+1] = NA
dailyHdrs[length(dailyD)] = "mxLWin"
dailyD[,length(dailyD)+1] = NA
dailyHdrs[length(dailyD)] = "mxLWout"
dailyD[,length(dailyD)+1] = NA
dailyHdrs[length(dailyD)] = "mnLWin"
dailyD[,length(dailyD)+1] = NA
dailyHdrs[length(dailyD)] = "mnLWout"
colnames(dailyD)<-dailyHdrs
#Loop per day
for(x in 1:length(dailyD[,1])){
print(paste("Processing date: ",dailyD[x,"Date"]))
#Determine which rows have the last observed maximum and minimum temperature
mxInd = which(data[,length(data)] == dailyD[x,"maxSearch"])
mxInd = mxInd[length(mxInd)]
mnInd = which(data[,length(data)] == dailyD[x,"minSearch"])
mnInd = mnInd[length(mnInd)]
#Extract solar and wind observations at these points
#At time of maximum
dailyD[x,"mxRgIn"] = data[mxInd,inShtRmavg] # Incoming Shortwave Radiation
dailyD[x,"mxRgOut"] = data[mxInd,outShtRmavg] # Outgoing Shortwave Radiation
dailyD[x,"mxLWin"] = data[mxInd,inCLngRmavg] # Incoming Longwave Radiation
dailyD[x,"mxLWout"] = data[mxInd,outCLngRmavg] # Incoming Longwave Radiation
dailyD[x,"mxUWind"] = data[mxInd,uwindmavg] # U component of the wind
dailyD[x,"mxVWind"] = data[mxInd,vwindmavg] # V component of the wind
#At time of minimum
dailyD[x,"mnRgIn"] = data[mnInd,inShtRmavg] # Incoming Shortwave Radiation
dailyD[x,"mnRgOut"] = data[mnInd,outShtRmavg] # Outgoing Shortwave Radiation
dailyD[x,"mnLWin"] = data[mnInd,inCLngRmavg] # Incoming Longwave Radiation
dailyD[x,"mnLWout"] = data[mnInd,outCLngRmavg] # Incoming Longwave Radiation
dailyD[x,"mnUWind"] = data[mnInd,uwindmavg] # U component of the wind
dailyD[x,"mnVWind"] = data[mnInd,vwindmavg] # V component of the wind
}
}
}
#TODO: REMOVE
print("Completed Tower processing!!")
## Preliminary anlaysis ##
hist(tDdailyMx[,7],xlab="Daily Maximum Difference MMTS-Nimbus minus MMTS-PRT", main="Tower D", xlim=c(-5,5), breaks=seq(from=-6,to=6,by=0.25),xaxs="i", yaxs="i")
plot(x=tDdailyMx[,2],y=tDdailyMx[,3],xlab="MMTS-Nimbus Daily Maximum (C)",ylab="MMTS-PRT Daily Maximum (C)",main="Tower D", xlim=c(-20,40), ylim=c(-20,40))
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="Date (YYYYMM)",
ylab=expression(bold("Daily Maximum Temperature")), main=pltTitles, ylim=c(-20,40))
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("topleft", legend=c("MMTS-PRT","MMTS-Nimbus"), ncol=2, fill, lty=1, lwd=3, col=c("blue","red"))
\ No newline at end of file
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