Commit eab8091d authored by Ronnie Leeper's avatar Ronnie Leeper
Browse files

Included MMTS flagging for out-of-bound observations of temps; -73 and 124 C

parent ac057a98
auto_roxygenize_for_build_and_reload="0"
auto_roxygenize_for_build_package="1"
auto_roxygenize_for_check="1"
makefile_args=""
{
}
\ No newline at end of file
~%2FDocuments%2FResearch%2FMMTS-BiasBudget%2FScripts%2FprocessTowerData.R="19B8543D"
{
"contents" : "#Process Tower Data\n\n#Created: Ronald D. Leeper\n#Date: 11.02.2015\n#Update: 11.13.2015\n\n#Purpose:\n#Read in tower station (A-D) data and create daily aggregates from the 16 second observations. From the tower data, generate\n#datasets for the respective biases:\n# 1). Sensor bias\n# 2). Sheilding bias\n# 3). Sitting bias (land cover)\n# 4). Time-of-observation (TOB) bias\n\n#11.13.2015 Update:\n# a). Included hourly moving averages of wind speed and solar radiation to evaluate atmospheric conditions at the time of observation\n# b). Set outlier observations from MMTS, -73.278 and 124.170 C, to NA\n\n#This script is archieved via gitlab at this URL: git@k3.cicsnc.org:ronnieleeper/MMTS-BiasBudget.git\n\n#Directory Setup\nSystemDir = \"/Users/ronnieleeper/Documents/\"\nProjectDir = \"Research/MMTS-BiasBudget/\"\nDataDir = \"Data/ATDD-TowerData/\"\nOutputDir = \"Data/\"\n\n#Tower IDs\ntowers = c(\"A\",\"B\",\"C\",\"D\")\n\n#Loop per tower (A-D)\nfor (tower in towers){\n print(paste(\"Processing tower: \",tower,sep=\"\"))\n \n #Readin current tower dataset and header file\n filename = paste(SystemDir,ProjectDir,DataDir,\"data\",tower,\".csv\",sep=\"\")\n data = read.table(file=filename, header=FALSE, sep=\",\")\n \n #Tower Header file\n filename = paste(SystemDir,ProjectDir,DataDir,\"hdrsTower\",tower,\".csv\",sep=\"\")\n hdrs = read.table(file=filename, sep=\",\", stringsAsFactors=FALSE)\n \n #Invoke column headers using the header file\n colnames(data)<-hdrs\n \n #Create unique date field (combine YYYYMMDD) from which to aggragate over\n #Covert the day field to a two-digit ID\n data[data[,\"Day\"] < 10,\"Day\"] = paste(\"0\",data[data[,\"Day\"] < 10,\"Day\"],sep=\"\")\n #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))\n data[,length(data)+1] = paste(data[,\"Year\"],data[,\"Month\"],data[,\"Day\"],sep=\"\")\n hdrs[length(hdrs)+1] = \"DailyID (12am Observer)\"\n colnames(data)<-hdrs\n \n #Identify the relative position (column number) of interested variables in the dataset\n if (tower == \"D\"){\n for (x in 1:length(hdrs)){\n if (hdrs[x] == \"NimbusTemp\") mmtsNimbus = x\n if (hdrs[x] == \"PRT1 (First CRN temperature)\") met1PRT1 = x\n if (hdrs[x] == \"PRT2 (Second CRN Temperature)\") met1PRT2 = x\n if (hdrs[x] == \"PRT3 (Third CRN temperature)\") met1PRT3 = x\n if (hdrs[x] == \"PRT_Nim (PRT inside of MMTS shelter)\") mmtsPRT = x\n if (hdrs[x] == \"Rgin (Incoming shortwave radiation)\") inShtR = x\n if (hdrs[x] == \"Rgout (Outgoing shortwave radiation)\") outShtR = x\n if (hdrs[x] == \"Lwin (Incoming Longwave Radiation Wm-2) \") inLngR = x\n if (hdrs[x] == \"Lwout (Outgoing Longwave Radiation Wm-2)\") outLngR = x\n if (hdrs[x] == \"KZTemp (4-ban net radiometer body temperature)\") KZTemp = x\n if (hdrs[x] == \"U (Sonic anem ws)\") uwind = x\n if (hdrs[x] == \"V (Sonic ws)\") vwind = x\n if (hdrs[x] == \"W (Sonic vertical ws)\") wwind = x\n if (hdrs[x] == \"batt_volt_Min\") battVol = x\n if (hdrs[x] == \"DailyID (12am Observer)\") dailyID = x\n }\n \n #Correct radiation values using body temperature (pulled from scar_view.m file)\n data[,length(data)+1] = data[,inLngR] + 5.67e-8*(data[,KZTemp]+273.15)^4\n hdrs[,length(data)] = \"correctedLwin\"\n cortedInLng = length(data)\n data[,length(data)+1] = data[,outLngR] + 5.67e-8*(data[,KZTemp]+273.15)^4\n hdrs[,length(data)] = \"correctedLwout\"\n cortedOutLng = length(data)\n colnames(data)<-hdrs\n \n #Set all observations to NA if battery voltages are less-than 11.2 volts\n minBatVolt = 11.2\n data[data[,battVol] < minBatVolt,c(mmtsNimbus,mmtsPRT,met1PRT1,met1PRT2,met1PRT3,inShtR,outShtR,inLngR,outLngR,uwind,vwind,wwind,cortedOutLng,cortedInLng)] = NA\n \n #Set outlier observations to NA\n minOutlier = -73.278\n maxOutlier = 124.170\n #Set Minimum outlier to NA\n data[data[,mmtsNimbus] == minOutlier & !is.na(data[,mmtsNimbus]), mmtsNimbus] = NA\n #Set Maximum outlier to NA\n data[data[,mmtsNimbus] == maxOutlier & !is.na(data[,mmtsNimbus]), mmtsNimbus] = NA\n \n #Calculate moving averages of windspeed and solar radition based on lag-times (60 minutes)\n sixtyMinuteMovingWindowCount = 226 #There are 225 16 second periods within a 60 minute window\n #Wind moving averages\n #Uwind\n #data[,length(data)+1] = SMA(data[,uwind], n=sixtyMinuteMovingWindowCount) #Calculate uwind moving average\n data[,length(data)+1] = rollapply(data[,uwind], width=sixtyMinuteMovingWindowCount, align=\"right\", FUN=mean, na.rm=T, fill=NA) #Calculate uwind moving average\n hdrs[length(data)] = \"uwind 60minMavg\" #Add variable to file header\n uwindmavg = length(hdrs) #Create a referece variable for calculated average\n #Vwind\n data[,length(data)+1] = rollapply(data[,vwind], width=sixtyMinuteMovingWindowCount, align=\"right\", FUN=mean, na.rm=T, fill=NA) #Calculate vwind moving average\n hdrs[length(data)] = \"vwind 60minMavg\" #Add variable to file header\n vwindmavg = length(hdrs) #Create a referece variable for calculated average\n \n #Solar moving averages\n #Incoming Shortwave Radiation\n data[,length(data)+1] = rollapply(data[,inShtR], width=sixtyMinuteMovingWindowCount, align=\"right\", FUN=mean, na.rm=T, fill=NA) #Calculate inShtR moving average\n hdrs[length(data)] = \"Rgin 60minMavg\" #Add variable to file header\n inShtRmavg = length(hdrs) #Create a referece variable for calculated average\n #Incoming Longwave Radiation\n data[,length(data)+1] = rollapply(data[,cortedInLng], width=sixtyMinuteMovingWindowCount, align=\"right\", FUN=mean, na.rm=T, fill=NA) \n hdrs[length(data)] = \"cLwin 60minMavg\" #Add variable to file header\n inCLngRmavg = length(hdrs) #Create a referece variable for calculated average\n #Outgoing Shortwave Radiation\n data[,length(data)+1] = rollapply(data[,outShtR], width=sixtyMinuteMovingWindowCount, align=\"right\", FUN=mean, na.rm=T, fill=NA)\n hdrs[length(data)] = \"Rgout 60minMavg\" #Add variable to file header\n outShtRmavg = length(hdrs) #Create a referece variable for calculated average\n #Outgoing Longwave Radiation\n data[,length(data)+1] = rollapply(data[,cortedOutLng], width=sixtyMinuteMovingWindowCount, align=\"right\", FUN=mean, na.rm=T, fill=NA)\n hdrs[length(data)] = \"cLout 60minMavg\" #Add variable to file header\n outCLngRmavg = length(hdrs) #Create a referece variable for calculated average\n \n #Compute PRT 5-minute running means\n fiveMinuteMovingWindowCount = 20 #There are 20 16-second periods within a five minute window\n #PRT1\n data[,length(data)+1] = SMA(data[,met1PRT1], n=fiveMinuteMovingWindowCount) #Calculate moving average\n hdrs[length(data)] = \"PRT1 5minMavg\" #Add variable to file header\n met1PRT1mavg = length(hdrs) #Create a referece variable for calculated average\n #PRT2\n data[,length(data)+1] = SMA(data[,met1PRT2], n=fiveMinuteMovingWindowCount) #Calculate moving average\n hdrs[length(data)] = \"PRT2 5minMavg\" #Add variable to file header\n met1PRT2mavg = length(hdrs) #Create a referece variable for calculated average\n #PRT3\n data[,length(data)+1] = SMA(data[,met1PRT3], n=fiveMinuteMovingWindowCount) #Calculate moving average\n hdrs[length(data)] = \"PRT3 5minMavg\" #Add variable to file header\n met1PRT3mavg = length(hdrs) #Create a referece variable for calculated average\n \n #Aggregate to daily values using dailyID\n #Maximum Aggregates\n tDdailyMx = aggregate(data[,c(mmtsNimbus,mmtsPRT,met1PRT1mavg,met1PRT2mavg,met1PRT3mavg)], by=list(data[,dailyID]), FUN=max, na.rm=T) #towerD Daily Maximum aggregates\n #Add prefix \"max\" to dataset headers\n dataHdrs = colnames(tDdailyMx)\n dataHdrs[2:length(dataHdrs)] = paste(\"max\",dataHdrs[2:length(dataHdrs)],sep=\"\")\n colnames(tDdailyMx)<-dataHdrs\n #Minimum Aggregates\n tDdailyMn = aggregate(data[,c(mmtsNimbus,mmtsPRT,met1PRT1mavg,met1PRT2mavg,met1PRT3mavg)], by=list(data[,dailyID]), FUN=min, na.rm=T) #towerD Daily Minimum aggregates\n #Add prefix \"min\" to dataset headers\n dataHdrs = colnames(tDdailyMn)\n dataHdrs[2:length(dataHdrs)] = paste(\"min\",dataHdrs[2:length(dataHdrs)],sep=\"\")\n colnames(tDdailyMn)<-dataHdrs\n }\n}\n",
"created" : 1447164802351.000,
"dirty" : false,
"encoding" : "UTF-8",
"folds" : "",
"hash" : "2945194069",
"id" : "10ADB5A9",
"lastKnownWriteTime" : 1447440117,
"path" : "~/Documents/Research/MMTS-BiasBudget/Scripts/processTowerData.R",
"project_path" : "Scripts/processTowerData.R",
"properties" : {
},
"relative_order" : 1,
"source_on_save" : false,
"type" : "r_source"
}
\ No newline at end of file
{
"contents" : "",
"created" : 1447419772710.000,
"dirty" : false,
"encoding" : "",
"folds" : "",
"hash" : "0",
"id" : "1D3D84CA",
"lastKnownWriteTime" : 140735261503232,
"path" : null,
"project_path" : null,
"properties" : {
"cacheKey" : "ksayulz9dr",
"caption" : "data",
"contentUrl" : "grid_resource/gridviewer.html?env=&obj=data&cache_key=ksayulz9dr",
"displayedObservations" : "2839964",
"environment" : "",
"object" : "data",
"totalObservations" : "2839964",
"variables" : "35"
},
"relative_order" : 5,
"source_on_save" : false,
"type" : "r_dataframe"
}
\ No newline at end of file
{
"contents" : "",
"created" : 1447426773235.000,
"dirty" : false,
"encoding" : "",
"folds" : "",
"hash" : "0",
"id" : "219652B9",
"lastKnownWriteTime" : 4864344320,
"path" : null,
"project_path" : null,
"properties" : {
"cacheKey" : "mfopsv6dwn",
"caption" : "data[1:300, c(1:6, uwind, 34)]",
"contentUrl" : "grid_resource/gridviewer.html?env=_rs_no_env&obj=&cache_key=mfopsv6dwn",
"displayedObservations" : 300,
"environment" : "_rs_no_env",
"object" : "",
"totalObservations" : 300,
"variables" : 8
},
"relative_order" : 6,
"source_on_save" : false,
"type" : "r_dataframe"
}
\ No newline at end of file
{
"contents" : "",
"created" : 1447178699036.000,
"dirty" : false,
"encoding" : "",
"folds" : "",
"hash" : "0",
"id" : "4A867CA0",
"lastKnownWriteTime" : 0,
"path" : null,
"project_path" : null,
"properties" : {
"cacheKey" : "w82ddm7zm8",
"caption" : "data[7400:7410, ]",
"contentUrl" : "grid_resource/gridviewer.html?env=_rs_no_env&obj=&cache_key=w82ddm7zm8",
"displayedObservations" : 11,
"environment" : "_rs_no_env",
"object" : "",
"totalObservations" : 11,
"variables" : 31
},
"relative_order" : 5,
"source_on_save" : false,
"type" : "r_dataframe"
}
\ No newline at end of file
{
"contents" : "",
"created" : 1447169040883.000,
"dirty" : false,
"encoding" : "",
"folds" : "",
"hash" : "0",
"id" : "63A78CB5",
"lastKnownWriteTime" : 7310030925004284416,
"path" : null,
"project_path" : null,
"properties" : {
"cacheKey" : "zstb4ggtap",
"caption" : "tDdailyMn",
"contentUrl" : "grid_resource/gridviewer.html?env=&obj=tDdailyMn&cache_key=zstb4ggtap",
"displayedObservations" : "528",
"environment" : "",
"object" : "tDdailyMn",
"totalObservations" : "528",
"variables" : "6"
},
"relative_order" : 2,
"source_on_save" : false,
"type" : "r_dataframe"
}
\ No newline at end of file
{
"contents" : "",
"created" : 1447169084618.000,
"dirty" : false,
"encoding" : "",
"folds" : "",
"hash" : "0",
"id" : "86572549",
"lastKnownWriteTime" : 4341388032,
"path" : null,
"project_path" : null,
"properties" : {
"cacheKey" : "chx32xmrz4",
"caption" : "tDdailyMx",
"contentUrl" : "grid_resource/gridviewer.html?env=&obj=tDdailyMx&cache_key=chx32xmrz4",
"displayedObservations" : 510,
"environment" : "",
"object" : "tDdailyMx",
"totalObservations" : 510,
"variables" : 6
},
"relative_order" : 3,
"source_on_save" : false,
"type" : "r_dataframe"
}
\ No newline at end of file
Version: 1.0
RestoreWorkspace: Default
SaveWorkspace: Default
AlwaysSaveHistory: Default
EnableCodeIndexing: Yes
UseSpacesForTab: Yes
NumSpacesForTab: 2
Encoding: UTF-8
RnwWeave: Sweave
LaTeX: pdfLaTeX
......@@ -2,7 +2,7 @@
#Created: Ronald D. Leeper
#Date: 11.02.2015
#Update: 11.10.2015
#Update: 11.13.2015
#Purpose:
#Read in tower station (A-D) data and create daily aggregates from the 16 second observations. From the tower data, generate
......@@ -12,6 +12,10 @@
# 3). Sitting bias (land cover)
# 4). Time-of-observation (TOB) bias
#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
#This script is archieved via gitlab at this URL: git@k3.cicsnc.org:ronnieleeper/MMTS-BiasBudget.git
#Directory Setup
......@@ -39,6 +43,8 @@ for (tower in towers){
colnames(data)<-hdrs
#Create unique date field (combine YYYYMMDD) from which to aggragate over
#Covert the day field to a two-digit ID
data[data[,"Day"] < 10,"Day"] = paste("0",data[data[,"Day"] < 10,"Day"],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="")
hdrs[length(hdrs)+1] = "DailyID (12am Observer)"
......@@ -61,15 +67,88 @@ for (tower in towers){
if (hdrs[x] == "V (Sonic ws)") vwind = x
if (hdrs[x] == "W (Sonic vertical ws)") wwind = x
if (hdrs[x] == "batt_volt_Min") battVol = x
if (hdrs[x] == "midObserDailyID") dailyID = x
if (hdrs[x] == "DailyID (12am Observer)") dailyID = x
}
#TODO:
#1). Correct longwave radiation values
#Aggregate by dailyID using max, min, and mean mathmatical methods
#Daily Maximum Values
dailyData = aggregate(data[,c()]
#Correct radiation values using body temperature (pulled from scar_view.m file)
data[,length(data)+1] = data[,inLngR] + 5.67e-8*(data[,KZTemp]+273.15)^4
hdrs[,length(data)] = "correctedLwin"
cortedInLng = length(data)
data[,length(data)+1] = data[,outLngR] + 5.67e-8*(data[,KZTemp]+273.15)^4
hdrs[,length(data)] = "correctedLwout"
cortedOutLng = length(data)
colnames(data)<-hdrs
#Set all observations to NA if battery voltages are less-than 11.2 volts
minBatVolt = 11.2
data[data[,battVol] < minBatVolt,c(mmtsNimbus,mmtsPRT,met1PRT1,met1PRT2,met1PRT3,inShtR,outShtR,inLngR,outLngR,uwind,vwind,wwind,cortedOutLng,cortedInLng)] = NA
#Set outlier observations to NA
minOutlier = -73.278
maxOutlier = 124.170
#Set Minimum outlier to NA
data[data[,mmtsNimbus] == minOutlier & !is.na(data[,mmtsNimbus]), mmtsNimbus] = NA
#Set Maximum outlier to NA
data[data[,mmtsNimbus] == maxOutlier & !is.na(data[,mmtsNimbus]), mmtsNimbus] = NA
#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
#Wind moving averages
#Uwind
#data[,length(data)+1] = SMA(data[,uwind], n=sixtyMinuteMovingWindowCount) #Calculate uwind moving average
data[,length(data)+1] = rollapply(data[,uwind], width=sixtyMinuteMovingWindowCount, align="right", FUN=mean, na.rm=T, fill=NA) #Calculate uwind moving average
hdrs[length(data)] = "uwind 60minMavg" #Add variable to file header
uwindmavg = length(hdrs) #Create a referece variable for calculated average
#Vwind
data[,length(data)+1] = rollapply(data[,vwind], width=sixtyMinuteMovingWindowCount, align="right", FUN=mean, na.rm=T, fill=NA) #Calculate vwind moving average
hdrs[length(data)] = "vwind 60minMavg" #Add variable to file header
vwindmavg = length(hdrs) #Create a referece variable for calculated average
#Solar moving averages
#Incoming Shortwave Radiation
data[,length(data)+1] = rollapply(data[,inShtR], width=sixtyMinuteMovingWindowCount, align="right", FUN=mean, na.rm=T, fill=NA) #Calculate inShtR moving average
hdrs[length(data)] = "Rgin 60minMavg" #Add variable to file header
inShtRmavg = length(hdrs) #Create a referece variable for calculated average
#Incoming Longwave Radiation
data[,length(data)+1] = rollapply(data[,cortedInLng], width=sixtyMinuteMovingWindowCount, align="right", FUN=mean, na.rm=T, fill=NA)
hdrs[length(data)] = "cLwin 60minMavg" #Add variable to file header
inCLngRmavg = length(hdrs) #Create a referece variable for calculated average
#Outgoing Shortwave Radiation
data[,length(data)+1] = rollapply(data[,outShtR], width=sixtyMinuteMovingWindowCount, align="right", FUN=mean, na.rm=T, fill=NA)
hdrs[length(data)] = "Rgout 60minMavg" #Add variable to file header
outShtRmavg = length(hdrs) #Create a referece variable for calculated average
#Outgoing Longwave Radiation
data[,length(data)+1] = rollapply(data[,cortedOutLng], width=sixtyMinuteMovingWindowCount, align="right", FUN=mean, na.rm=T, fill=NA)
hdrs[length(data)] = "cLout 60minMavg" #Add variable to file header
outCLngRmavg = length(hdrs) #Create a referece variable for calculated average
#Compute PRT 5-minute running means
fiveMinuteMovingWindowCount = 20 #There are 20 16-second periods within a five minute window
#PRT1
data[,length(data)+1] = SMA(data[,met1PRT1], n=fiveMinuteMovingWindowCount) #Calculate moving average
hdrs[length(data)] = "PRT1 5minMavg" #Add variable to file header
met1PRT1mavg = length(hdrs) #Create a referece variable for calculated average
#PRT2
data[,length(data)+1] = SMA(data[,met1PRT2], n=fiveMinuteMovingWindowCount) #Calculate moving average
hdrs[length(data)] = "PRT2 5minMavg" #Add variable to file header
met1PRT2mavg = length(hdrs) #Create a referece variable for calculated average
#PRT3
data[,length(data)+1] = SMA(data[,met1PRT3], n=fiveMinuteMovingWindowCount) #Calculate moving average
hdrs[length(data)] = "PRT3 5minMavg" #Add variable to file header
met1PRT3mavg = length(hdrs) #Create a referece variable for calculated average
#Aggregate to daily values using dailyID
#Maximum Aggregates
tDdailyMx = aggregate(data[,c(mmtsNimbus,mmtsPRT,met1PRT1mavg,met1PRT2mavg,met1PRT3mavg)], by=list(data[,dailyID]), FUN=max, na.rm=T) #towerD Daily Maximum aggregates
#Add prefix "max" to dataset headers
dataHdrs = colnames(tDdailyMx)
dataHdrs[2:length(dataHdrs)] = paste("max",dataHdrs[2:length(dataHdrs)],sep="")
colnames(tDdailyMx)<-dataHdrs
#Minimum Aggregates
tDdailyMn = aggregate(data[,c(mmtsNimbus,mmtsPRT,met1PRT1mavg,met1PRT2mavg,met1PRT3mavg)], by=list(data[,dailyID]), FUN=min, na.rm=T) #towerD Daily Minimum aggregates
#Add prefix "min" to dataset headers
dataHdrs = colnames(tDdailyMn)
dataHdrs[2:length(dataHdrs)] = paste("min",dataHdrs[2:length(dataHdrs)],sep="")
colnames(tDdailyMn)<-dataHdrs
}
}
Supports Markdown
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