Mercurial > repos > ecology > argo_getdata

comparison divandfull.jl @ 1:28c8f085d1b4 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
planemo upload for repository https://github.com/galaxyecology/tools-ecology/tree/master/tools/ocean commit e395cfee9cab90bbed58ac52fb8467c896f51824
author ecology
date Thu, 01 Aug 2024 09:46:52 +0000
parents
children
comparison
equal deleted inserted replaced
0:055a934a380f 1:28c8f085d1b4
1 #Julia script
2
3 ###############################
4 ## DIVAndrun analsysis ##
5 ###############################
6 import Pkg;
7 using Pkg
8 Pkg.status()
9
10 ### Import packages
11 using DIVAnd
12 using Dates
13 using Printf
14 # Getting the arguments from the command line
15 args = ARGS
16
17 # Import data
18 if length(args) < 4
19 error("This tool needs at least 4 arguments")
20 else
21 netcdf_data = args[1]
22 longmin = parse(Float64, args[2])
23 longmax = parse(Float64, args[3])
24 latmin = parse(Float64, args[4])
25 latmax = parse(Float64, args[5])
26 startdate = args[6] # yyyy,mm,dd
27 enddate = args[7]
28 varname = args[8]
29 selmin = parse(Float64, args[9])
30 selmax = parse(Float64, args[10])
31 bathname = args[11]
32 end
33
34 ## This script will create a climatology:
35 # 1. ODV data reading.
36 # 2. Extraction of bathymetry and creation of mask
37 # 3. Data download from other sources and duplicate removal.
38 # 4. Quality control.
39 # 5. Parameter optimisation.
40 # 6. Spatio-temporal interpolation with DIVAnd.
41
42
43 ### Configuration
44 # Define the horizontal, vertical (depth levels) and temporal resolutions.
45 # Select the variable of interest
46
47 dx, dy = 0.125, 0.125
48 lonr = longmin:dx:longmax
49 latr = latmin:dy:latmax
50
51 # Convert string in date
52 startdate = Date(startdate, "yyyy-mm-dd")
53
54 # extract year month day
55 startyear = year(startdate)
56 startmonth = month(startdate)
57 startday = day(startdate)
58
59 # Convert string in date
60 enddate = Date(enddate, "yyyy-mm-dd")
61
62 # extract year month day
63 endyear = year(enddate)
64 endmonth = month(enddate)
65 endday = day(enddate)
66
67 timerange = [Date(startyear, startmonth, startday),Date(endyear, endmonth, endday)];
68
69 depthr = [0.,5., 10., 15., 20., 25., 30., 40., 50., 66,
70 75, 85, 100, 112, 125, 135, 150, 175, 200, 225, 250,
71 275, 300, 350, 400, 450, 500, 550, 600, 650, 700, 750,
72 800, 850, 900, 950, 1000, 1050, 1100, 1150, 1200, 1250,
73 1300, 1350, 1400, 1450, 1500, 1600, 1750, 1850, 2000];
74 depthr = [0.,10.,20.];
75
76 varname = varname
77 yearlist = [1900:2023];
78 monthlist = [[1,2,3],[4,5,6],[7,8,9],[10,11,12]];
79
80 # We create here the variable TS (for "tDataset(netcdf_data,"r")ime selector"), which allows us to work with the observations corresponding to each period of interest.
81
82 TS = DIVAnd.TimeSelectorYearListMonthList(yearlist,monthlist);
83 @show TS;
84
85 figdir = "outputs/"
86 if ~(isdir(figdir))
87 mkdir(figdir)
88 else
89 @info("Figure directory already exists")
90 end
91 ### 1. Read your ODV file
92 # Adapt the datadir and datafile values.
93 # The example is based on a sub-setting of the Mediterranean Sea aggregated dataset.
94 # The dataset has been extracted around the Adriatic Sea and exported to a netCDF using Ocean Data
95 datadir = "../data"
96
97 datafile = netcdf_data
98
99 # Then you can read the full file:
100 @time obsval,obslon,obslat,obsdepth,obstime,obsid = NCODV.load(Float64, datafile,
101 "Water body $(varname)");
102
103 # Check the extremal values of the observations
104 checkobs((obslon,obslat,obsdepth,obstime),obsval,obsid)
105
106 ### 2. Extract the bathymetry
107
108 # It is used to delimit the domain where the interpolation is performed.
109 ## 2.1 Choice of bathymetry
110
111 # Modify bathname according to the resolution required.
112
113 @time bx,by,b = load_bath(bathname,true,lonr,latr);
114
115 ## 2.2 Create mask
116 # False for sea
117 # True for land
118
119 mask = falses(size(b,1),size(b,2),length(depthr))
120 for k = 1:length(depthr)
121 for j = 1:size(b,2)
122 for i = 1:size(b,1)
123 mask[i,j,k] = b[i,j] >= depthr[k]
124 end
125 end
126 end
127 @show size(mask)
128
129 ### 3. Quality control
130 # We check the salinity value.
131 # Adapt the criteria to your region and variable.
132
133 sel = (obsval .<= selmax) .& (obsval .>= selmin);
134
135 obsval = obsval[sel]
136 obslon = obslon[sel]
137 obslat = obslat[sel]
138 obsdepth = obsdepth[sel]
139 obstime = obstime[sel]
140 obsid = obsid[sel];
141
142 ### 4. Analysis parameters
143 # Correlation lengths and noise-to-signal ratio
144
145 # We will use the function diva3D for the calculations.
146 # With this function, the correlation length has to be defined in meters, not in degrees.
147
148 sz = (length(lonr),length(latr),length(depthr));
149 lenx = fill(100_000.,sz) # 100 km
150 leny = fill(100_000.,sz) # 100 km
151 lenz = fill(25.,sz); # 25 m
152 len = (lenx, leny, lenz);
153 epsilon2 = 0.1;
154
155 ### Output file name
156 outputdir = "outputs_netcdf/"
157 if !isdir(outputdir)
158 mkpath(outputdir)
159 end
160 filename = joinpath(outputdir, "Water_body_$(replace(varname," "=>"_")).nc")
161
162 ### 7. Analysis
163 # Remove the result file before running the analysis, otherwise you'll get the message
164 if isfile(filename)
165 rm(filename) # delete the previous analysis
166 @info "Removing file $filename"
167 end
168
169 ## 7.1 Plotting function
170 # Define a plotting function that will be applied for each time index and depth level.
171 # All the figures will be saved in a selected directory.
172
173 function plotres(timeindex,sel,fit,erri)
174 tmp = copy(fit)
175 nx,ny,nz = size(tmp)
176 for i in 1:nz
177 figure("Additional-Data")
178 ax = subplot(1,1,1)
179 ax.tick_params("both",labelsize=6)
180 ylim(39.0, 46.0);
181 xlim(11.5, 20.0);
182 title("Depth: (timeindex)", fontsize=6)
183 pcolor(lonr.-dx/2.,latr.-dy/2, permutedims(tmp[:,:,i], [2,1]);
184 vmin = 33, vmax = 40)
185 colorbar(extend="both", orientation="vertical", shrink=0.8).ax.tick_params(labelsize=8)
186
187 contourf(bx,by,permutedims(b,[2,1]), levels = [-1e5,0],colors = [[.5,.5,.5]])
188 aspectratio = 1/cos(mean(latr) * pi/180)
189 gca().set_aspect(aspectratio)
190
191 figname = varname * @sprintf("_%02d",i) * @sprintf("_%03d.png",timeindex)
192 plt.savefig(joinpath(figdir, figname), dpi=600, bbox_inches="tight");
193 plt.close_figs()
194 end
195 end
196
197 ## 7.2 Create the gridded fields using diva3d
198 # Here only the noise-to-signal ratio is estimated.
199 # Set fitcorrlen to true to also optimise the correlation length.
200 @time dbinfo = DIVAnd.diva3d((lonr,latr,depthr,TS),
201 (obslon,obslat,obsdepth,obstime), obsval,
202 len, epsilon2,
203 filename,varname,
204 bathname=bathname,
205 fitcorrlen = false,
206 niter_e = 2,
207 surfextend = true
208 );
209
210 # Save the observation metadata in the NetCDF file.
211 DIVAnd.saveobs(filename,(obslon,obslat,obsdepth,obstime),obsid);