diff divandfull.jl @ 0:484930fdc002 draft

planemo upload for repository https://github.com/galaxyecology/tools-ecology/tree/master/tools/ocean commit e395cfee9cab90bbed58ac52fb8467c896f51824

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/divandfull.jl	Thu Aug 01 09:46:44 2024 +0000
@@ -0,0 +1,211 @@
+#Julia script
+
+###############################
+##    DIVAndrun analsysis    ##
+###############################
+import Pkg; 
+using Pkg
+Pkg.status()
+
+### Import packages
+using DIVAnd
+using Dates
+using Printf
+# Getting the arguments from the command line
+args = ARGS
+
+# Import data
+if length(args) < 4
+    error("This tool needs at least 4 arguments")
+else
+    netcdf_data = args[1]
+    longmin = parse(Float64, args[2])
+    longmax = parse(Float64, args[3])
+    latmin = parse(Float64, args[4])
+    latmax = parse(Float64, args[5])
+    startdate = args[6] # yyyy,mm,dd
+    enddate = args[7]
+    varname = args[8]
+    selmin = parse(Float64, args[9])
+    selmax = parse(Float64, args[10])
+    bathname = args[11]
+end
+
+## This script will create a climatology:
+# 1. ODV data reading.
+# 2. Extraction of bathymetry and creation of mask
+# 3. Data download from other sources and duplicate removal.
+# 4. Quality control.
+# 5. Parameter optimisation.
+# 6. Spatio-temporal interpolation with DIVAnd.
+
+
+### Configuration
+# Define the horizontal, vertical (depth levels) and temporal resolutions.
+# Select the variable of interest
+
+dx, dy = 0.125, 0.125
+lonr = longmin:dx:longmax
+latr = latmin:dy:latmax
+
+# Convert string in date
+startdate = Date(startdate, "yyyy-mm-dd")
+
+# extract year month day
+startyear = year(startdate)
+startmonth = month(startdate)
+startday = day(startdate)
+
+# Convert string in date
+enddate = Date(enddate, "yyyy-mm-dd")
+
+# extract year month day
+endyear = year(enddate)
+endmonth = month(enddate)
+endday = day(enddate)
+
+timerange = [Date(startyear, startmonth, startday),Date(endyear, endmonth, endday)];
+
+depthr = [0.,5., 10., 15., 20., 25., 30., 40., 50., 66, 
+    75, 85, 100, 112, 125, 135, 150, 175, 200, 225, 250, 
+    275, 300, 350, 400, 450, 500, 550, 600, 650, 700, 750, 
+    800, 850, 900, 950, 1000, 1050, 1100, 1150, 1200, 1250, 
+    1300, 1350, 1400, 1450, 1500, 1600, 1750, 1850, 2000];
+depthr = [0.,10.,20.];
+
+varname = varname
+yearlist = [1900:2023];
+monthlist = [[1,2,3],[4,5,6],[7,8,9],[10,11,12]];
+
+# 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.
+
+TS = DIVAnd.TimeSelectorYearListMonthList(yearlist,monthlist);
+@show TS;
+
+figdir = "outputs/"
+if ~(isdir(figdir))
+    mkdir(figdir)
+else
+    @info("Figure directory already exists")
+end
+### 1. Read your ODV file
+# Adapt the datadir and datafile values.
+# The example is based on a sub-setting of the Mediterranean Sea aggregated dataset.
+# The dataset has been extracted around the Adriatic Sea and exported to a netCDF using Ocean Data 
+datadir = "../data"
+
+datafile = netcdf_data
+
+# Then you can read the full file:
+@time obsval,obslon,obslat,obsdepth,obstime,obsid = NCODV.load(Float64, datafile, 
+    "Water body $(varname)");
+
+# Check the extremal values of the observations
+checkobs((obslon,obslat,obsdepth,obstime),obsval,obsid)
+
+### 2. Extract the bathymetry
+
+# It is used to delimit the domain where the interpolation is performed.
+## 2.1 Choice of bathymetry
+
+# Modify bathname according to the resolution required.
+
+@time bx,by,b = load_bath(bathname,true,lonr,latr);
+
+## 2.2 Create mask
+# False for sea
+# True for land
+
+mask = falses(size(b,1),size(b,2),length(depthr))
+for k = 1:length(depthr)
+    for j = 1:size(b,2)
+        for i = 1:size(b,1)
+            mask[i,j,k] = b[i,j] >= depthr[k]
+        end
+    end
+end
+@show size(mask)
+
+### 3. Quality control
+# We check the salinity value.
+# Adapt the criteria to your region and variable.
+
+sel = (obsval .<= selmax) .& (obsval .>= selmin);
+
+obsval = obsval[sel]
+obslon = obslon[sel]
+obslat = obslat[sel]
+obsdepth = obsdepth[sel]
+obstime = obstime[sel]
+obsid = obsid[sel];
+
+### 4. Analysis parameters
+# Correlation lengths and noise-to-signal ratio
+
+# We will use the function diva3D for the calculations.
+# With this function, the correlation length has to be defined in meters, not in degrees.
+
+sz = (length(lonr),length(latr),length(depthr));
+lenx = fill(100_000.,sz)   # 100 km
+leny = fill(100_000.,sz)   # 100 km
+lenz = fill(25.,sz);      # 25 m 
+len = (lenx, leny, lenz);
+epsilon2 = 0.1;
+
+### Output file name
+outputdir = "outputs_netcdf/"
+if !isdir(outputdir)
+    mkpath(outputdir)
+end
+filename = joinpath(outputdir, "Water_body_$(replace(varname," "=>"_")).nc")
+
+### 7. Analysis
+# Remove the result file before running the analysis, otherwise you'll get the message
+if isfile(filename)
+    rm(filename) # delete the previous analysis
+    @info "Removing file $filename"
+end
+
+## 7.1 Plotting function
+# Define a plotting function that will be applied for each time index and depth level.
+# All the figures will be saved in a selected directory.
+     
+function plotres(timeindex,sel,fit,erri)
+    tmp = copy(fit)
+    nx,ny,nz = size(tmp)
+    for i in 1:nz
+        figure("Additional-Data")
+        ax = subplot(1,1,1)
+        ax.tick_params("both",labelsize=6)
+        ylim(39.0, 46.0);
+        xlim(11.5, 20.0);
+        title("Depth: (timeindex)", fontsize=6)
+        pcolor(lonr.-dx/2.,latr.-dy/2, permutedims(tmp[:,:,i], [2,1]);
+               vmin = 33, vmax = 40)
+        colorbar(extend="both", orientation="vertical", shrink=0.8).ax.tick_params(labelsize=8)
+
+        contourf(bx,by,permutedims(b,[2,1]), levels = [-1e5,0],colors = [[.5,.5,.5]])
+        aspectratio = 1/cos(mean(latr) * pi/180)
+        gca().set_aspect(aspectratio)
+        
+        figname = varname * @sprintf("_%02d",i) * @sprintf("_%03d.png",timeindex)
+        plt.savefig(joinpath(figdir, figname), dpi=600, bbox_inches="tight");
+        plt.close_figs()
+    end
+end
+
+## 7.2 Create the gridded fields using diva3d
+# Here only the noise-to-signal ratio is estimated.
+# Set fitcorrlen to true to also optimise the correlation length.
+@time dbinfo = DIVAnd.diva3d((lonr,latr,depthr,TS),
+    (obslon,obslat,obsdepth,obstime), obsval,
+    len, epsilon2,
+    filename,varname,
+    bathname=bathname,
+    fitcorrlen = false,
+    niter_e = 2,
+    surfextend = true
+    );
+
+# Save the observation metadata in the NetCDF file.
+DIVAnd.saveobs(filename,(obslon,obslat,obsdepth,obstime),obsid);