-
Notifications
You must be signed in to change notification settings - Fork 2
/
Plotting_diagnostic_plots_for_site_selection.R
188 lines (108 loc) · 4.64 KB
/
Plotting_diagnostic_plots_for_site_selection.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
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
library(ncdf4)
#clear R environment
rm(list=ls(all=TRUE))
#Set path
path <- "/srv/ccrc/data04/z3509830/Fluxnet_data/All_flux_sites_processed_PLUMBER2/"
#Source function
source(paste0(path, "/scripts/functions/plot_timeseries.R"))
#Variables to plot
vars <- c("SWdown", "LWdown", "Precip", "Tair", "Qair", "Wind", "CO2air",
"Rnet", "Qle", "Qh", "Qg", "Ebal", "LAI_MODIS", "LAI_Copernicus")
#Variable types
var_type <- c("Met", "Met", "Met", "Met", "Met", "Met", "Met",
"Flux", "Flux", "Flux", "Flux", "Flux", "Met", "Met")
#Find sites
site_files <- list.files(paste0(path, "/all_sites_no_duplicates/Nc_files/Met_with_LAI/"),
full.names=TRUE)
#Open file handles
site_nc_met <- lapply(site_files, nc_open)
#Get site codes
site_codes <- sapply(site_nc_met, function(x) ncatt_get(x, varid=0, "site_code")$value)
#Set output directory
outdir <- paste0(path, "/PLUMBER2_site_selection_figs/")
dir.create(outdir)
#Loop through sites
for (s in 1:length(site_codes)) {
#Get matching flux file
nc_flux <- nc_open(list.files(paste0(path, "/all_sites_no_duplicates/Nc_files/Flux"),
pattern=site_codes[s], full.names=TRUE))
#Get time (used to var lengths etc)
time <- ncvar_get(nc_flux, "time")
#Get timing information
timing <- GetTimingNcfile(nc_flux)
#QC information
QC_measured <- 0
if(substr(site_codes[s], 1, 3) == "AU-") QC_measured <- append(QC_measured, 10)
#Set up file
outfile <- paste0(outdir, "/", site_codes[s], "_timeseries_plot.png")
png(outfile, height=1300, width=3500, res=50, pointsize=20)
par(mfcol=c(3,5))
#Loop through variables
for (v in 1:length(vars)) {
#Get data
#Energy balance
if (vars[v] == "Ebal") {
rnet <- tryCatch(ncvar_get(nc_flux, "Rnet"), error=function(e) rep(0, length(time)))
qle <- tryCatch(ncvar_get(nc_flux, "Qle"), error=function(e) rep(0, length(time)))
qh <- tryCatch(ncvar_get(nc_flux, "Qh"), error=function(e) rep(0, length(time)))
#Look for ground heat flux, set to zero if not found
qg <- tryCatch(ncvar_get(nc_flux, "Qg"), error=function(e) rep(0, length(time)))
#Calculate ratio (Rnet-G) / (Qle + Qh)
var_data <- (rnet - qg) / (qle + qh)
qc_data <- NA
data_unit <- "-"
#Cap values to 0-5
var_data[var_data > 1.5] <- 1.5
var_data[var_data < 0] <- 0
var_data[which(is.infinite(var_data))] <- 1.5
#All other variables
} else {
#Set file for met variable
if (var_type[v] == "Met") {
nc <- site_nc_met[[s]]
#Set file for flux variable
} else if (var_type[v] == "Flux") {
nc <- nc_flux
}
#Get variable data
var_data <- tryCatch(ncvar_get(nc, vars[v]), error=function(e) rep(NA, length(time)))
#Get QC data
qc_data <- tryCatch(ncvar_get(nc, paste0(vars[v], "_qc")), error=function(e) NA)
#Get data units
data_unit <- tryCatch(ncatt_get(nc, varid=vars[v], "units")$value, error=function(e) NA)
}
#Convert precip and Tair units
if (vars[v] == "Precip") {
var_data <- var_data * (time[2] - time[1])
data_unit <- paste("mm/", (time[2] - time[1])/60, "min", sep="")
}
if (vars[v] == "Tair") {
var_data <- var_data - 273.15
data_unit <- "deg C"
}
#Extract QC data and replace all gap-filled values with 0
# and measured with 1 (opposite to Fluxnet but what PALS expects)
if (any(!is.na(qc_data))) {
qc_data[!(qc_data %in% QC_measured)] <- 2 #replace gap-filled values with a temporary value
qc_data[qc_data %in% QC_measured] <- 1 #set measured to 1
qc_data[qc_data == 2] <- 0 #set gap-filled to 0
#If first value missing, set to measured (to avoid an error when PALS
#checks if first value -1)
if(is.na(qc_data[1])){ qc_data[1] <- 0}
#Else set to PALS option corresponding to no QC data
} else {
qc_data <- matrix(-1, nrow = 1, ncol = 1)
}
#Plot file
Timeseries(obslabel=site_codes[s], tsdata=as.matrix(var_data),
varname=vars[v],
ytext=paste(vars[v], " (", data_unit, ")", sep=""),
legendtext=vars[v],
plotcex=2 , timing=timing,
smoothed = FALSE, winsize = 1,
plotcolours="black",
vqcdata = as.matrix(qc_data),
na.rm=TRUE)
} #vars
dev.off()
} #sites