Mercurial > repos > ecology > argo_getdata
changeset 1:28c8f085d1b4 draft default tip
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 | 055a934a380f |
| children | |
| files | divandfull.jl test-data/data_from_Eutrophication_Med_profiles_2022_unrestricted.nc |
| diffstat | 2 files changed, 211 insertions(+), 0 deletions(-) [+] |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/divandfull.jl Thu Aug 01 09:46:52 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);