-
Notifications
You must be signed in to change notification settings - Fork 0
/
plot_TS.R
105 lines (81 loc) · 2.75 KB
/
plot_TS.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
#############################################################################
# This script does the time series plot
#############################################################################
plot_TS <- function(M, PARAM_NAME, plot_name, zoom_pres=NULL, zoom_param=NULL, zoom_x=NULL, date_axis=FALSE, logscale=FALSE, max_qc=2) {
### find array dimensions
n_prof = dim(M)[2]
N_DEPTH = rep(NA, n_prof)
for (i in seq(1,n_prof)) {
N_DEPTH[i] = length(M[,i]$PRES)
}
n_depth = max(N_DEPTH)
### build arrays from M
pres = array(NA, c(n_depth, n_prof))
param = array(NA, c(n_depth, n_prof))
param_qc = array(NA, c(n_depth, n_prof))
juld = array(NA, c(n_depth, n_prof))
param_units = rep(NA, n_prof)
profile_id = array(NA, c(n_depth, n_prof))
for (i in seq(1,n_prof)) {
n_pres_i = length(M[,i]$PRES)
pres[1:n_pres_i, i] = M[,i]$PRES
param[1:n_pres_i, i] = M[,i]$PARAM
param_qc[1:n_pres_i, i] = M[,i]$PARAM_QC
juld[1:n_pres_i, i] = rep(M[,i]$JULD, n_pres_i)
profile_id[1:n_pres_i, i] = rep(M[,i]$profile_id, n_pres_i)
param_units[i] = M[,i]$param_units
}
# transform to vector
pres = as.vector(pres)
param = as.vector(param)
param_qc = as.vector(param_qc)
juld = as.vector(juld)
profile_id = as.vector(profile_id)
if (!is.null(zoom_pres)) {
pres[which( pres<zoom_pres[1] | pres>zoom_pres[2] )] = NA
}
if (!is.null(zoom_param)) {
param[which( param<zoom_param[1] )] = zoom_param[1]
param[which( param>zoom_param[2] )] = zoom_param[2]
}
if (logscale) {
param = log10(param)
}
# remove all data with QC worse than max_qc (but not data with QC>=5)
param[which( param_qc>max_qc & param_qc<=4 )] = NA
not_na_axis = which( !is.na(pres) & is.finite(param) )
profile_id = profile_id[not_na_axis]
juld = juld[not_na_axis]
param_qc = param_qc[not_na_axis]
param = param[not_na_axis]
pres = pres[not_na_axis]
# define labeling parameters
dates = as.Date(juld, origin='1950-01-01')
if (date_axis) {
Xaxis = dates
Xlabel = "Time"
} else {
Xaxis = profile_id
Xlabel = "Profile number"
}
param_units = unique(param_units[which(!is.na(param_units))])
if (logscale) {
param_units = paste("log10( ", param_units, " )")
}
print(range(param))
colors = colormap(param, zlim=c(min(param), max(param)))
# create plot
png(plot_name, width = 800, height = 400)
layout(t(1:2), widths=c(10,1))
par(mar=c(5,4,4,0.5))
plot(Xaxis, pres, col=colors$zcol, pch=16, cex=0.5, xlim=zoom_x, ylim=rev(range(pres, na.rm=T)), xlab=Xlabel, ylab="Pressure (decibar)")
title(main=PARAM_NAME)
#image.plot(juld,pres,param)
# colorbar
par(mar=c(5,0,4,3.5))
image(y=colors$breaks, z=t(colors$breaks), col=colors$col, axes=FALSE)
axis(4, cex.axis=0.8, mgp=c(0,0.5,0))
mtext(param_units, side=4, line=2)
dev.off()
return(0)
}