# HG changeset patch # User ecology # Date 1722505604 0 # Node ID 484930fdc002e3c0cd84b960067cf7ddfe73d6e4 planemo upload for repository https://github.com/galaxyecology/tools-ecology/tree/master/tools/ocean commit e395cfee9cab90bbed58ac52fb8467c896f51824 diff -r 000000000000 -r 484930fdc002 argo_getdata.py --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/argo_getdata.py Thu Aug 01 09:46:44 2024 +0000 @@ -0,0 +1,191 @@ +# author: Marie Jossé + +# Python script + +############################# +# Argo data access # +############################# + +# Packages : argopy + + +# Load arguments +import argparse +import sys + +import argopy + +command_line_args = sys.argv[1:] + + +parser = argparse.ArgumentParser(description="Retrieve argo Data") +# Add arguments + +parser.add_argument("--user", type=str, + help="User mode : standard, expert or research") +parser.add_argument("--cardinal_1", type=float, help="Longitude min") +parser.add_argument("--cardinal_2", type=float, help="Longitude max") +parser.add_argument("--cardinal_3", type=float, help="Latitude min") +parser.add_argument("--cardinal_4", type=float, help="Latitude max") +parser.add_argument("--pressure_1", type=float, help="Pressure min") +parser.add_argument("--pressure_2", type=float, help="Pressure max") +parser.add_argument("--date_1", type=str, help="Starting date") +parser.add_argument("--date_2", type=str, help="Ending date.") +parser.add_argument("--wmo", type=str, help="WMO: argo's identifier") +parser.add_argument("--profile", type=str, help="Number of profiles") +parser.add_argument("--params", type=str, help="List of bgc parameters") +parser.add_argument("--measured", type=str, help="List of bgc parameters") +parser.add_argument("--output_argo", type=str, help="Output data from argo") + +args = parser.parse_args(command_line_args) + + +# Parse the command line arguments + +print(args) +# Import data + +user = args.user +cardinal_1 = args.cardinal_1 +cardinal_2 = args.cardinal_2 +cardinal_3 = args.cardinal_3 +cardinal_4 = args.cardinal_4 +pressure_1 = args.pressure_1 +pressure_2 = args.pressure_2 +date_1 = args.date_1 +date_2 = args.date_2 +wmo = args.wmo +if wmo is not None: + wmo = list(map(int, wmo.split(","))) +profile = args.profile +if profile is not None: + profile = list(map(int, profile.split(","))) +params = args.params +if params is not None: + params = params.split(",") + if len(params) == 83: + params = "all" +measured = args.measured +if measured is not None: + measured = measured.split(",") + +# Let’s import the argopy data fetcher: + +###################### +# User mode # +###################### +# By default, +# all argopy data fetchers are set to work with a standard user mode. +# To change that + +argopy.set_options(mode=user) + +###################### +# Fetching Argo data # +###################### +# Data selection # + +# To access Argo data with a DataFetcher, +# you need to define how to select your data of interest. +# argopy provides 3 different data selection methods: + +argo_data = argopy.DataFetcher() + +# 🗺 For a space/time domain # + +if (cardinal_1 is not None): + mode = "region" + argo_data = argo_data.region([cardinal_1, cardinal_2, + cardinal_3, cardinal_4, + pressure_1, pressure_2, + date_1, date_2]) + +# ⚓ For one or more profiles # +# Use the fetcher access point argopy.DataFetcher.profile() +# to specify the float WMO platform number +# and the profile cycle number(s) to retrieve profiles for. +elif (wmo is not None and profile is not None): + argo_data = argo_data.profile(wmo, profile) + # can also be argo_data = argo_data.profile(6902755, [3, 12]) + mode = "profile" + +# 🤖 For one or more floats # +# If you know the Argo float unique identifier number called a WMO number +# you can use the fetcher access point DataFetcher.float() +# to specify one or more float WMO platform numbers to select. +else: + argo_data = argo_data.float(wmo) + # can also be argo_data = argo_data.float([6902746, 6902755]) + mode = "float" + +# Data sources # +# Let’s start with standard import: +# argopy.reset_options() +# Specify data source erddap, gdac or argovis + +# if (ftp != "") : + # argopy.set_options(src = "gdac", ftp = ftp) +# else : + # argopy.set_options(src = "erddap") + +# With remote, online data sources, +# it may happens that the data server is experiencing down time. +print(argopy.status()) + +# Dataset # +# Argo data are distributed as a single dataset. +# It is referenced at https://doi.org/10.17882/42182. +# But they are several Argo missions with specific files and parameters +# that need special handling by argopy, namely: +# - the core Argo Mission: from floats that measure temperature, +# salinity, pressure down to 2000m, +# - the Deep Argo Mission: from floats that measure temperature, +# salinity, pressure down to 6000m, +# - and the BGC-Argo Mission: from floats that measure temperature, +# salinity, pressure and oxygen, pH, nitrate, chlorophyll, +# backscatter, irradiance down to 2000m. +# You can choose between phy or bgc +if (params is None): + argopy.set_options(dataset="phy") +else: + argopy.set_options(dataset="bgc") + if (measured != ['None'] and measured is not None): + argo_data = argopy.DataFetcher(params=params, measured=measured) + if (mode == "region"): + argo_data = argo_data.region([cardinal_1, cardinal_2, + cardinal_3, cardinal_4, + pressure_1, pressure_2, + date_1, date_2]) + elif (mode == "profile"): + argo_data = argo_data.profile(wmo, profile) + else: + argo_data = argo_data.float(wmo) + else: + argo_data = argopy.DataFetcher(params=params, measured=None) + if (mode == "region"): + argo_data = argo_data.region([cardinal_1, cardinal_2, + cardinal_3, cardinal_4, + pressure_1, pressure_2, + date_1, date_2]) + elif (mode == "profile"): + argo_data = argo_data.profile(wmo, profile) + else: + argo_data = argo_data.float(wmo) + +# Data fetching # +# To fetch (i.e. access, download, format) Argo data, +# argopy provides the DataFetcher class. +# Several DataFetcher arguments exist to help you select the dataset, +# the data source and the user mode the most suited for your applications; +# and also to improve performances. + +# You define the selection of data you want to fetch +# with one of the DataFetcher methods: region, float or profile. +# 2 lines to download Argo data: import and fetch ! + +argo_data = argo_data.load().data +argo_data.to_netcdf("argo_data.nc") + +# argo_metadata = argo_data.to_index() + +print(argo_data) diff -r 000000000000 -r 484930fdc002 divandfull.jl --- /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); diff -r 000000000000 -r 484930fdc002 divandfull.xml --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/divandfull.xml Thu Aug 01 09:46:44 2024 +0000 @@ -0,0 +1,140 @@ + + Data-Interpolating Variational Analysis in n dimensions + + 0.1.0 + 0 + + + julia + julia-divand + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + [0-9a-zA-Z_]+ + + + + + + + + + + + + + + + + + + + + + + + + + + + + `_ but +not further developed anymore). + +The method bears some similarities and equivalences with Optimal Interpolation or Krigging in that it allows to create a smooth +and continous field from a collection of observations, observations which can be affected by errors. The analysis method is however +different in practise, allowing to take into account topological features, physical constraints etc in a natural way. +The method was initially developped with ocean data in mind, but it can be applied to any field where localized observations have +to be used to produce gridded fields which are "smooth". + +DIVAndrun is the core analysis function in n dimensions. It does not know anything about the physical parameters or units you work with. +Coordinates can also be very general. The only constraint is that the metrics (pm,pn,po,...) when multiplied by the corresponding length +scales len lead to non-dimensional parameters. Furthermore the coordinates of the output grid (xi,yi,zi,...) need to have the same units +as the observation coordinates (x,y,z,...). + +DIVAndfun is a version with a minimal set of parameters (the coordinates and values of observations, i.e. (x,f), the remaining parameters +being optional) and provides an interpolation function rather than an already gridded field. + +diva3D is a higher-level function specifically designed for climatological analysis of data on Earth, using longitude/latitude/depth/time +coordinates and correlations length in meters. It makes the necessary preparation of metrics, parameter optimizations etc you normally would +program yourself before calling the analysis function DIVAndrun. + +DIVAnd_heatmap can be used for additive data and produces Kernel Density Estimations. + +DIVAndgo is only needed for very large problems when a call to DIVAndrun leads to memory or CPU time problems. This function tries to decide +which solver (direct or iterative) to use and how to make an automatic domain decomposition. Not all options from DIVAndrun are available. + +If you want to try out multivariate approaches, you can look at DIVAnd_multivarEOF and DIVAnd_multivarJAC + +If you want more informations about the functions and parameters see also the `documentations here `_. + + ]]> + + doi:10.5194/gmd-7-225-2014 + + diff -r 000000000000 -r 484930fdc002 macro.xml --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/macro.xml Thu Aug 01 09:46:44 2024 +0000 @@ -0,0 +1,24 @@ + + 0.1.15 + 0 + + + argopy + + + + + + + + + + + + + + + doi:10.21105/joss.02425 + + + diff -r 000000000000 -r 484930fdc002 test-data/argo_data.netcdf Binary file test-data/argo_data.netcdf has changed diff -r 000000000000 -r 484930fdc002 test-data/data_from_Eutrophication_Med_profiles_2022_unrestricted.nc Binary file test-data/data_from_Eutrophication_Med_profiles_2022_unrestricted.nc has changed