/
derived.jl
154 lines (126 loc) · 5.1 KB
/
derived.jl
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
function deepest(array::Array{Union{Missing,T},3}) where T
#function deepest(array::AbstractArray{Union{Missing,T},3}) where T
sz = size(array)
array2d = Array{Union{Missing,T},3-1}(undef,sz[1],sz[2])
array2d .= missing
for k = 1:sz[3]
for j = 1:sz[2]
for i = 1:sz[1]
if !ismissing(array[i,j,k])
# possibly overwrite previous value
array2d[i,j] = array[i,j,k]
end
end
end
end
return array2d
end
function ncdeepest(nc,varname,suffix,thresholds_value)
T = Float32
ncvar = nc[varname * suffix]
sz = size(ncvar)
@assert length(sz) == 4
dims = ("lon", "lat", "time")
newvarname = varname * "_deepest" * suffix
longname = ncvar.attrib["long_name"]
valex = ncvar.attrib["_FillValue"]
longname_deepest = "Deepest values of $(longname)"
if thresholds_value !== nothing
longname_deepest = longname_deepest * " masked using relative error threshold $(thresholds_value)"
end
# value of deepest layer
@info("Creating new variable $(newvarname)")
ncvar_deepest = defVar(nc, newvarname, T, dims, attrib = OrderedDict(
"_FillValue" => T(valex),
"missing_value" => T(valex),
"units" => ncvar.attrib["units"],
"long_name" => "Deepest values of $(longname)",
))
for n = 1:sz[4]
field3D = nc[varname][:,:,:,n]
ncvar_deepest[:,:,n] = deepest(field3D)
end
end
"""
DIVAnd.derived(filename,varname,new_filename;
error_thresholds = [("L1", 0.3), ("L2", 0.5)],
)
Compute derived quantities from a DIVAnd analyse (in particular the deepest
value of an analysis) using the NetCDF file `filename` with the variable
`varname`. The result will be written in `new_filename`.
See `DIVAnd.diva3d` for the optional parameter `error_thresholds`.
Example:
```julia
using DIVAnd
filename = "Water_body_dissolved_oxygen_concentration_monthly.nc"
new_filename = "Water_body_dissolved_oxygen_concentration_monthly2.nc"
varname = "Water body dissolved oxygen concentration"
DIVAnd.derived(filename,varname,new_filename)
```
"""
function derived(source_filename,varname,datafile;
error_thresholds = [("L1", 0.3), ("L2", 0.5)],
)
if source_filename !== datafile
cp(source_filename,datafile, force=true)
end
NCDatasets.Dataset(datafile, "a") do nc
ncvar = nc[varname]
ncvar_relerr = nc[varname * "_relerr"]
longname = ncvar.attrib["long_name"]
T = Float32
sz = size(ncvar)
ncattrib = OrderedDict(ncvar.attrib)
ncdeepest(nc,varname,"",nothing)
dims = ("lon", "lat", "time")
#@show ncvar[10,3,:,1]
#@show nc[varname * "_deepest"][10,3,1]
for (thresholds_name, thresholds_value) in error_thresholds
newvarname = "$(varname)_deepest_$(thresholds_name)"
ncattrib["long_name"] = "Deepest values of $(longname) masked using relative error threshold $(thresholds_value)"
ncvar_Lx = defVar(nc, newvarname, T, dims, attrib = ncattrib)
for n = 1:sz[4]
field = deepest(ncvar[:,:,:,n])
relerr = deepest(ncvar_relerr[:,:,:,n])
for ind in eachindex(relerr)
if !ismissing(relerr[ind])
if relerr[ind] .> thresholds_value
field[ind] = missing
end
end
end
ncvar_Lx[:,:,n] = field
end
end
# add depth of deepest layer
# Load the depth variable
depth = nc["depth"][:]
@info("Working on $(length(depth)) depth layers")
# Load the 4D variable (we can load the first time instance,
# as depth don't change with time)
field3D = nc[varname][:,:,:,1]
@info("size = $(size(field3D))")
# Get missing value from field
valex = nc[varname].attrib["missing_value"]
newvarname = varname * "_deepest_depth"
@info("Creating new variable $(newvarname)")
ncvardeepestdepth = defVar(nc, newvarname, T, ("lon", "lat"))
ncvardeepestdepth.attrib["long_name"] = "Deepest depth for $(varname)"
ncvardeepestdepth.attrib["_FillValue"] = T(valex)
ncvardeepestdepth.attrib["missing_value"] = T(valex)
ncvardeepestdepth.attrib["units"] = "meters"
ncvardeepestdepth.attrib["positive"] = "down"
# Loop on depth: start from surface and go to the bottom
# (I also add "abs" in case depth are negative, but not probable)
depthindex = sortperm(abs.(depth))
for idepth in depthindex
@info("Depth index: $(idepth)")
# Look for non-missing values at the considered depth
nonmissing = (.!ismissing.(field3D[:,:,idepth]))
@info("Found $(sum(nonmissing)) non missing values for depth $(depth[idepth])")
# Write the variable
ncvardeepestdepth[nonmissing] .= depth[idepth]
end
@info("Written new variable deepest depth")
end
end