forked from NEONScience/EC_Validity
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathLoadAnalyzeECTEFiles.R
131 lines (105 loc) · 5.2 KB
/
LoadAnalyzeECTEFiles.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
##install required packages
#install.packages("neonUtilities")
#install.packages("tidyverse")
#install.packages("BiocManager")
#BiocManager::install('rhdf5')
##strings not automatically assigned to factor and load packages
options(stringsAsFactors=F)
library(neonUtilities)
library(ggplot2)
library(tidyverse)
library(lubridate)
##Enter NEON site code and date range
sitename <- "STEI"
startdate <- "2019-01"
enddate <- "2019-12"
dirtemp <- tempdir()
if (!dir.exists(paste0(dirtemp, "/", sitename))) {
dirfile <- dir.create(paste0(dirtemp, "/", sitename))
}
##Download eddy co bundled data package
zipsByProduct(dpID="DP4.00200.001", package="expanded",
site=c(sitename),
startdate=startdate, enddate=enddate,
savepath=paste0(dirtemp, "/", sitename),
check.size=F)
##Stack and save to list of data frames
list_flux <- stackEddy(paste0(dirtemp, "/", sitename, "/filesToStack00200"),
level = "dp04")
##save to data frame
flux <- list_flux[[1]]
##Save flux data frame locally
#saveRDS(list_flux$UNDE, "C:/Users/cslemmons/Documents/ExploreNEON/EDDY/filesToStack00200/UNDEflux.rds")
##Load flux data
#flux <- readRDS(file="C:/Users/cslemmons/Documents/ExploreNEON/EDDY/filesToStack00200/UNDEflux.rds")
##Convert time stamp to R date - time format, add to data frame
flux$TimeB <- as.POSIXct(flux$timeBgn,
format="%Y-%m-%dT%H:%M:%S",
tz="GMT")
##Add year and month column for summary
flux$YearMonth <- substr(flux$TimeB,0,7)
##Add year, month, day column for summary
flux$YearMonthDay = substr(flux$TimeB,0,10)
##Evaluate fluxCor and fluxRaw column for NAs, calculate % NAs by month
NAsPerMonth <- flux %>%
mutate(FluxCorNA=is.na(data.fluxCo2.turb.fluxCor)) %>%
mutate(FluxRawNA=is.na(data.fluxCo2.turb.fluxRaw)) %>%
group_by(YearMonth) %>%
summarise(FluxCorNAPercent = mean(FluxCorNA*100), FluxRawNAPercent = mean(FluxRawNA*100))
##Evaluate fluxCor and fluxRaw column for NAs, calculate % NAs by day
NAsPerDay <- flux %>%
mutate(FluxCorNA=is.na(data.fluxCo2.turb.fluxCor)) %>%
mutate(FluxRawNA=is.na(data.fluxCo2.turb.fluxRaw)) %>%
group_by(YearMonthDay) %>%
summarise(FluxCorNAPercent = mean(FluxCorNA*100), FluxRawNAPercent = mean(FluxRawNA*100))
##calculate difference between raw and corrected flux, save to data frame. Calculate mean difference by month
flux$DiffCor <- flux$data.fluxCo2.turb.fluxRaw - flux$data.fluxCo2.turb.fluxCor
fluxDiffbyMonth <- flux %>%
mutate(DiffCor=data.fluxCo2.turb.fluxRaw-data.fluxCo2.turb.fluxCor) %>%
group_by(YearMonth) %>%
summarise(DiffCor = mean(DiffCor, na.rm=TRUE))
##calculate difference between raw and corrected flux and calculate mean difference grouped by day
fluxDiffbyDay <- flux %>%
mutate(DiffCor=data.fluxCo2.turb.fluxRaw-data.fluxCo2.turb.fluxCor) %>%
group_by(YearMonthDay) %>%
summarise(DiffCor = mean(DiffCor, na.rm=TRUE))
##calculate number of final Turb QF by Month and Day
finQFbyMonth<-flux %>%
group_by(YearMonth) %>%
summarise(QFPercent=mean(qfqm.fluxCo2.turb.qfFinl*100))
finQFbyDay<-flux %>%
group_by(YearMonthDay) %>%
summarise(QFPercent=mean(qfqm.fluxCo2.turb.qfFinl*100))
##Plot raw and corrected turbulent flux data for whole dataset
ggplot(flux) +geom_point(aes(x=TimeB, y=data.fluxCo2.turb.fluxRaw, color="blue"))+
geom_point(aes(x=TimeB, y=data.fluxCo2.turb.fluxCor, color="green"))+
labs(title=paste0(sitename, " 2019 Turbulent Flux"))+
xlab("Date")+
ylab("CO2 umol m2/s-1")+
scale_color_identity(name=" ",
breaks = c("blue", "green"),
labels = c("Raw Flux","Corrected Flux"),
guide = "legend")+ theme(plot.title = element_text(hjust=0.5),legend.position = c(.9,.1), legend.background = element_rect(color=NULL, fill="transparent"))
##Plot mean difference between raw and corrected flux per Month
ggplot(fluxDiffbyMonth) +geom_point(aes(x=YearMonth, y=DiffCor))+
theme(axis.text.x=element_text(angle=90))+
labs(title=paste0(sitename, " 2019 Difference Between Raw and Corrected Flux"), x="Date", y="CO2 umol m2/s-1")
##Plot mean difference between raw and corrected flux per day
ggplot(fluxDiffbyDay) +geom_point(aes(x=YearMonthDay, y=DiffCor))+
theme(axis.text.x=element_text(angle=90))+
scale_x_discrete(breaks=c("2019-01-02","2019-02-02","2019-03-01","2019-04-01","2019-05-01", "2019-06-01", "2019-07-01", "2019-08-01", "2019-09-01", "2019-10-01", "2019-11-01", "2019-12-01"))+
labs(title=paste0(sitename," 2019 Mean Daily Difference Between Raw and Corrected Flux"))+
xlab("Date")+
ylab("CO2 umol m2/s-1")
##Plot Percent NA and Percent QF by Month
ggplot(NAsPerMonth) + aes(x=YearMonth, y=FluxCorNAPercent)+
geom_col()+
ylim(0,100)+
labs(x = "Date", y = "Percent", title = paste0(sitename, " 2019 - Mean Percent of Missing Corrected Turbulent Flux"))
ggplot(finQFbyMonth) + aes(x=YearMonth, y=QFPercent)+
geom_col()+
ylim(0,100)+
labs(x = "Date", y = "Percent", title = paste0(sitename, " 2019 - Mean Percent of Quality Flags Turbulent Flux"))
##save to CSV
#savefilename <- "C:/Users/cslemmons/Documents/ExploreNEON/filesToStack00200/UNDEFlux.csv"
#write.table(flux, "C:/Users/cslemmons/Documents/ExploreNEON/filesToStack00200/UNDEFlux.csv", sep=",")