|
|
|
#pragma once
|
|
|
|
#include "DataAdapters/ncfilealt.h"
|
|
|
|
#include "ParseArgs.h"
|
|
|
|
#include "mdatetime.h"
|
|
|
|
#include "simple2ddata.h"
|
|
|
|
#include "specfunc.h"
|
|
|
|
#include <memory>
|
|
|
|
#include <regex.h>
|
|
|
|
#include <set>
|
|
|
|
#include <stdlib.h>
|
|
|
|
#include <sys/types.h>
|
|
|
|
|
|
|
|
using michlib::Ceil;
|
|
|
|
using michlib::errmessage;
|
|
|
|
using michlib::Floor;
|
|
|
|
using michlib::GPL;
|
|
|
|
using michlib::int2;
|
|
|
|
using michlib::MDateTime;
|
|
|
|
using michlib::NCFileA;
|
|
|
|
using michlib::ToGeoDomainNeg;
|
|
|
|
|
|
|
|
class NEMOData
|
|
|
|
{
|
|
|
|
enum Type
|
|
|
|
{
|
|
|
|
TYPE_UNKNOWN,
|
|
|
|
TYPE_DT,
|
|
|
|
TYPE_NRT,
|
|
|
|
TYPE_NRT6
|
|
|
|
};
|
|
|
|
|
|
|
|
struct Parameters: public BaseParameters
|
|
|
|
{
|
|
|
|
size_t xb, yb, xe, ye, layer;
|
|
|
|
MString varname;
|
|
|
|
virtual ~Parameters() override = default;
|
|
|
|
};
|
|
|
|
|
|
|
|
std::vector<NCFileA> nc;
|
|
|
|
std::vector<real> depths;
|
|
|
|
std::vector<MDateTime> times;
|
|
|
|
size_t nx, ny;
|
|
|
|
real lonb, latb, lone, late;
|
|
|
|
real lonstep, latstep;
|
|
|
|
Type type = TYPE_UNKNOWN;
|
|
|
|
|
|
|
|
class EnvVar
|
|
|
|
{
|
|
|
|
MString name, oldvalue;
|
|
|
|
bool activated, saved;
|
|
|
|
|
|
|
|
public:
|
|
|
|
EnvVar(): activated(false) {}
|
|
|
|
~EnvVar() { Deactivate(); }
|
|
|
|
|
|
|
|
void Activate(const MString& var, const MString& val)
|
|
|
|
{
|
|
|
|
if(activated) Deactivate();
|
|
|
|
name = var;
|
|
|
|
char* curval = getenv(name.Buf());
|
|
|
|
if(nullptr == curval)
|
|
|
|
saved = false;
|
|
|
|
else
|
|
|
|
{
|
|
|
|
oldvalue = curval;
|
|
|
|
saved = true;
|
|
|
|
}
|
|
|
|
setenv(name.Buf(), val.Buf(), 1);
|
|
|
|
}
|
|
|
|
void Deactivate()
|
|
|
|
{
|
|
|
|
if(!activated) return;
|
|
|
|
if(saved)
|
|
|
|
setenv(name.Buf(), oldvalue.Buf(), 1);
|
|
|
|
else
|
|
|
|
unsetenv(name.Buf());
|
|
|
|
activated = false;
|
|
|
|
}
|
|
|
|
};
|
|
|
|
|
|
|
|
EnvVar proxy;
|
|
|
|
|
|
|
|
public:
|
|
|
|
using Data = Simple2DData;
|
|
|
|
|
|
|
|
private:
|
|
|
|
template<class DataType> Data ReadVarRaw(const NCFileA& f, const MString& name, size_t i, bool nodepth, const Parameters* p) const
|
|
|
|
{
|
|
|
|
real unitmul = 1.0;
|
|
|
|
DataType fill;
|
|
|
|
real offset = 0.0, scale = 1.0;
|
|
|
|
|
|
|
|
{
|
|
|
|
auto a_fill = f.A<DataType>(name, "_FillValue");
|
|
|
|
auto a_offset = f.A<double>(name, "add_offset");
|
|
|
|
auto a_scale = f.A<double>(name, "scale_factor");
|
|
|
|
if(!a_fill) return Data();
|
|
|
|
fill = a_fill;
|
|
|
|
if(a_offset) offset = a_offset;
|
|
|
|
if(a_scale) scale = a_scale;
|
|
|
|
}
|
|
|
|
|
|
|
|
auto unit = f.A<MString>(name, "units");
|
|
|
|
if(unit && unit.Get() == "m s-1") unitmul = 100.0;
|
|
|
|
|
|
|
|
Data data((p->xb < p->xe) ? (p->xe - p->xb + 1) : (nx + p->xe - p->xb + 1), p->ye - p->yb + 1, Lon(p->xb), Lat(p->yb), lonstep, latstep);
|
|
|
|
|
|
|
|
if(p->xb < p->xe)
|
|
|
|
{
|
|
|
|
auto var = nodepth ? f.V<DataType>(name, {"longitude", p->xb, p->xe - p->xb + 1}, {"latitude", p->yb, p->ye - p->yb + 1}, {"time", i, 1})
|
|
|
|
: f.V<DataType>(name, {"longitude", p->xb, p->xe - p->xb + 1}, {"latitude", p->yb, p->ye - p->yb + 1}, {"time", i, 1}, {"depth", p->layer, 1});
|
|
|
|
if(!var) return Data();
|
|
|
|
if(var.DimLen(0) != data.Nx() || var.DimLen(1) != data.Ny()) return Data();
|
|
|
|
|
|
|
|
for(size_t ix = 0; ix < var.DimLen(0); ix++)
|
|
|
|
for(size_t iy = 0; iy < var.DimLen(1); iy++)
|
|
|
|
{
|
|
|
|
DataType v = var(ix, iy);
|
|
|
|
data(ix, iy) = (v == fill) ? Data::Fillval() : ((v * scale + offset) * unitmul);
|
|
|
|
}
|
|
|
|
}
|
|
|
|
else
|
|
|
|
{
|
|
|
|
auto var1 = nodepth ? f.V<DataType>(name, {"longitude", p->xb}, {"latitude", p->yb, p->ye - p->yb + 1}, {"time", i, 1})
|
|
|
|
: f.V<DataType>(name, {"longitude", p->xb}, {"latitude", p->yb, p->ye - p->yb + 1}, {"time", i, 1}, {"depth", p->layer, 1});
|
|
|
|
auto var2 = nodepth ? f.V<DataType>(name, {"longitude", 0, p->xe + 1}, {"latitude", p->yb, p->ye - p->yb + 1}, {"time", i, 1})
|
|
|
|
: f.V<DataType>(name, {"longitude", 0, p->xe + 1}, {"latitude", p->yb, p->ye - p->yb + 1}, {"time", i, 1}, {"depth", p->layer, 1});
|
|
|
|
if(!(var1 && var2)) return Data();
|
|
|
|
if((var1.DimLen(0) + var2.DimLen(0)) != data.Nx() || var1.DimLen(1) != data.Ny() || var2.DimLen(1) != data.Ny()) return Data();
|
|
|
|
for(size_t ix = 0; ix < var1.DimLen(0); ix++)
|
|
|
|
for(size_t iy = 0; iy < var1.DimLen(1); iy++)
|
|
|
|
{
|
|
|
|
DataType v = var1(ix, iy);
|
|
|
|
data(ix, iy) = (v == fill) ? Data::Fillval() : ((v * scale + offset) * unitmul);
|
|
|
|
}
|
|
|
|
for(size_t ix = 0; ix < var2.DimLen(0); ix++)
|
|
|
|
for(size_t iy = 0; iy < var2.DimLen(1); iy++)
|
|
|
|
{
|
|
|
|
DataType v = var2(ix, iy);
|
|
|
|
data(ix + var1.DimLen(0), iy) = (v == fill) ? Data::Fillval() : ((v * scale + offset) * unitmul);
|
|
|
|
}
|
|
|
|
}
|
|
|
|
return data;
|
|
|
|
}
|
|
|
|
|
|
|
|
static MString VName2IName(const MString& vname)
|
|
|
|
{
|
|
|
|
if(vname == "temp") return "thetao";
|
|
|
|
if(vname == "sal") return "so";
|
|
|
|
if(vname == "mld") return "mlotst";
|
|
|
|
if(vname == "ssh") return "zos";
|
|
|
|
if(vname == "w") return "wo";
|
|
|
|
return "";
|
|
|
|
}
|
|
|
|
|
|
|
|
public:
|
|
|
|
Data Read(const BaseParameters* ip, size_t i) const
|
|
|
|
{
|
|
|
|
if(!isOk()) return Data();
|
|
|
|
bool nodepth = false;
|
|
|
|
|
|
|
|
auto p = dynamic_cast<const struct Parameters*>(ip);
|
|
|
|
MString name = VName2IName(p->varname);
|
|
|
|
for(const auto& f: nc)
|
|
|
|
{
|
|
|
|
auto head = f.Header();
|
|
|
|
for(const auto& v: head.Variables())
|
|
|
|
if(v.Name() == name)
|
|
|
|
{
|
|
|
|
if(v.Dimensions().size() == 3) nodepth = true;
|
|
|
|
if(v.Type().Id() == NC_SHORT) return ReadVarRaw<int2>(f, name, i, nodepth, p);
|
|
|
|
if(v.Type().Id() == NC_FLOAT) return ReadVarRaw<float>(f, name, i, nodepth, p);
|
|
|
|
}
|
|
|
|
}
|
|
|
|
return Data();
|
|
|
|
}
|
|
|
|
|
|
|
|
bool HasVar(const MString vname) const
|
|
|
|
{
|
|
|
|
MString name = VName2IName(vname);
|
|
|
|
if(name == "") return false;
|
|
|
|
|
|
|
|
for(const auto& f: nc)
|
|
|
|
{
|
|
|
|
auto head = f.Header();
|
|
|
|
for(const auto& v: head.Variables())
|
|
|
|
if(v.Name() == name) return true;
|
|
|
|
}
|
|
|
|
return false;
|
|
|
|
}
|
|
|
|
|
|
|
|
public:
|
|
|
|
NEMOData() = default;
|
|
|
|
|
|
|
|
// TODO: RetVal
|
|
|
|
MString Open(const CLArgs& args)
|
|
|
|
{
|
|
|
|
MString oldprefix = GPL.UsePrefix("AVISO");
|
|
|
|
MString dataset = args.contains("dataset") ? args.at("dataset") : "DT";
|
|
|
|
MString cred = GPL.ParameterSValue("COPERNICUS_USER", "");
|
|
|
|
MString proxyurl = GPL.ParameterSValue("COPERNICUS_PROXY", "");
|
|
|
|
|
|
|
|
GPL.UsePrefix(oldprefix);
|
|
|
|
|
|
|
|
nc.clear();
|
|
|
|
if(proxyurl.Exist()) proxy.Activate("all_proxy", proxyurl);
|
|
|
|
|
|
|
|
std::vector<MString> urls;
|
|
|
|
if(dataset == "DT")
|
|
|
|
{
|
|
|
|
urls = {"https://" + cred + "@my.cmems-du.eu/thredds/dodsC/cmems_mod_glo_phy_my_0.083_P1D-m"};
|
|
|
|
type = TYPE_DT;
|
|
|
|
}
|
|
|
|
else if(dataset == "NRT")
|
|
|
|
{
|
|
|
|
urls = {"https://" + cred + "@nrt.cmems-du.eu/thredds/dodsC/cmems_mod_glo_phy-cur_anfc_0.083deg_P1D-m",
|
|
|
|
"https://" + cred + "@nrt.cmems-du.eu/thredds/dodsC/cmems_mod_glo_phy-thetao_anfc_0.083deg_P1D-m",
|
|
|
|
"https://" + cred + "@nrt.cmems-du.eu/thredds/dodsC/cmems_mod_glo_phy-so_anfc_0.083deg_P1D-m",
|
|
|
|
"https://" + cred + "@nrt.cmems-du.eu/thredds/dodsC/cmems_mod_glo_phy_anfc_0.083deg_P1D-m",
|
|
|
|
"https://" + cred + "@nrt.cmems-du.eu/thredds/dodsC/cmems_mod_glo_phy-wcur_anfc_0.083deg_P1D-m"};
|
|
|
|
type = TYPE_NRT;
|
|
|
|
}
|
|
|
|
else if(dataset == "NRT6")
|
|
|
|
{
|
|
|
|
urls = {"https://" + cred + "@nrt.cmems-du.eu/thredds/dodsC/cmems_mod_glo_phy-cur_anfc_0.083deg_PT6H-i",
|
|
|
|
"https://" + cred + "@nrt.cmems-du.eu/thredds/dodsC/cmems_mod_glo_phy-thetao_anfc_0.083deg_PT6H-i",
|
|
|
|
"https://" + cred + "@nrt.cmems-du.eu/thredds/dodsC/cmems_mod_glo_phy-so_anfc_0.083deg_PT6H-i"};
|
|
|
|
type = TYPE_NRT6;
|
|
|
|
}
|
|
|
|
else
|
|
|
|
return "Unknown dataset: " + dataset;
|
|
|
|
|
|
|
|
nc.resize(urls.size());
|
|
|
|
for(size_t i = 0; i < urls.size(); i++)
|
|
|
|
{
|
|
|
|
nc[i].Reset(urls[i]);
|
|
|
|
if(!nc[i])
|
|
|
|
{
|
|
|
|
nc.clear();
|
|
|
|
return "Can't connect to database";
|
|
|
|
}
|
|
|
|
}
|
|
|
|
|
|
|
|
auto rdepths = nc[0].V<float>("depth");
|
|
|
|
if(!rdepths)
|
|
|
|
{
|
|
|
|
nc.clear();
|
|
|
|
return "Can't read depths";
|
|
|
|
}
|
|
|
|
depths.resize(rdepths.DimLen(0));
|
|
|
|
for(size_t i = 0; i < depths.size(); i++) depths[i] = rdepths(i);
|
|
|
|
|
|
|
|
auto timeD = nc[0].V<double>("time");
|
|
|
|
auto timeF = nc[0].V<float>("time");
|
|
|
|
if(!(timeD || timeF))
|
|
|
|
{
|
|
|
|
nc.clear();
|
|
|
|
return "Can't read times";
|
|
|
|
}
|
|
|
|
MDateTime refdate("1950-01-01");
|
|
|
|
timeD ? times.resize(timeD.DimLen(0)) : times.resize(timeF.DimLen(0));
|
|
|
|
for(size_t i = 0; i < times.size(); i++) times[i] = refdate + static_cast<time_t>(timeD ? timeD(i) : timeF(i)) * 3600;
|
|
|
|
|
|
|
|
auto nlon = nc[0].D("longitude");
|
|
|
|
auto nlat = nc[0].D("latitude");
|
|
|
|
if(!(nlon && nlat))
|
|
|
|
{
|
|
|
|
nc.clear();
|
|
|
|
return "Can't get number of longitudes/latitudes";
|
|
|
|
}
|
|
|
|
nx = nlon;
|
|
|
|
ny = nlat;
|
|
|
|
|
|
|
|
auto lons = nc[0].V<float>("longitude");
|
|
|
|
auto lats = nc[0].V<float>("latitude");
|
|
|
|
if(!(lons && lats))
|
|
|
|
{
|
|
|
|
nc.clear();
|
|
|
|
return "Can't get longitudes/latitudes";
|
|
|
|
}
|
|
|
|
lonb = lons(0);
|
|
|
|
latb = lats(0);
|
|
|
|
lone = lons(nx - 1);
|
|
|
|
late = lats(ny - 1);
|
|
|
|
lonstep = (lone - lonb) / (nx - 1);
|
|
|
|
latstep = (late - latb) / (ny - 1);
|
|
|
|
|
|
|
|
return "";
|
|
|
|
}
|
|
|
|
|
|
|
|
MString Info() const
|
|
|
|
{
|
|
|
|
if(!isOk()) return "";
|
|
|
|
MString d;
|
|
|
|
for(size_t i = 0; i < NDepths(); i++) d += MString(" ") + "(" + i + " " + Depth(i) + ")";
|
|
|
|
|
|
|
|
std::set<MString> vars;
|
|
|
|
for(const auto& f: nc)
|
|
|
|
{
|
|
|
|
auto head = f.Header();
|
|
|
|
for(const auto& v: head.Variables())
|
|
|
|
{
|
|
|
|
if(v.Name() == "thetao") vars.emplace("temp");
|
|
|
|
if(v.Name() == "so") vars.emplace("sal");
|
|
|
|
if(v.Name() == "mlotst") vars.emplace("mld");
|
|
|
|
if(v.Name() == "zos") vars.emplace("ssh");
|
|
|
|
if(v.Name() == "wo") vars.emplace("w");
|
|
|
|
}
|
|
|
|
}
|
|
|
|
MString svars;
|
|
|
|
{
|
|
|
|
bool first = true;
|
|
|
|
for(const auto& v: vars)
|
|
|
|
{
|
|
|
|
svars += (first ? "" : ", ") + v;
|
|
|
|
first = false;
|
|
|
|
}
|
|
|
|
}
|
|
|
|
|
|
|
|
// clang-format off
|
|
|
|
return
|
|
|
|
"Dataset: " + Title() + "\n" +
|
|
|
|
" Begin date: " + Time(0).ToString() + "\n" +
|
|
|
|
" End date: " + Time(NTimes()-1).ToString() + "\n" +
|
|
|
|
" Time step: " + Timestep() + " seconds\n" +
|
|
|
|
" Time moments: " + NTimes() + "\n" +
|
|
|
|
" Region: (" + lonb + " : " + lone + ") x (" + latb + " : " + late + ")\n" +
|
|
|
|
" Grid: " + nx + "x" + ny + " (" + lonstep + " x " + latstep + ")\n" +
|
|
|
|
" Depths:" + d + "\n" +
|
|
|
|
" Supported variables: " + svars;
|
|
|
|
// clang-format on
|
|
|
|
}
|
|
|
|
|
|
|
|
std::pair<const BaseParameters*, MString> Parameters(michlib_internal::ParameterListEx& pars, const CLArgs& args) const
|
|
|
|
{
|
|
|
|
std::unique_ptr<struct Parameters> ppar{new struct Parameters};
|
|
|
|
|
|
|
|
if(!args.contains("var")) return {nullptr, "Variable not specified"};
|
|
|
|
ppar->varname = args.at("var");
|
|
|
|
if(!HasVar(ppar->varname)) return {nullptr, "Variable " + ppar->varname + " not exists in this dataset"};
|
|
|
|
if(args.contains("layer")) ppar->layer = args.at("layer").ToInteger<size_t>();
|
|
|
|
if(!args.contains("depth") && ppar->layer >= NDepths()) return {nullptr, MString("Layer ") + ppar->layer + " is too deep!"};
|
|
|
|
real depth = args.contains("depth") ? args.at("depth").ToReal() : Depth(ppar->layer);
|
|
|
|
|
|
|
|
if(!(args.contains("lonb") && args.contains("lone") && args.contains("latb") && args.contains("late"))) return {nullptr, "Region not specified (lonb, lone, latb, late)"};
|
|
|
|
|
|
|
|
{
|
|
|
|
real lon1 = ToGeoDomainNeg(args.at("lonb").ToReal());
|
|
|
|
real lon2 = ToGeoDomainNeg(args.at("lone").ToReal());
|
|
|
|
real lat1 = ToGeoDomainNeg(args.at("latb").ToReal());
|
|
|
|
real lat2 = ToGeoDomainNeg(args.at("late").ToReal());
|
|
|
|
|
|
|
|
ppar->yb = static_cast<size_t>(Floor((lat1 - latb) / latstep));
|
|
|
|
ppar->ye = static_cast<size_t>(Ceil((lat2 - latb) / latstep));
|
|
|
|
if(ppar->ye > 2040) ppar->ye = 2040;
|
|
|
|
if(ppar->yb >= ppar->ye) return {nullptr, "Latb must be lesser then late"};
|
|
|
|
|
|
|
|
ppar->xb = static_cast<size_t>(Floor((lon1 - lonb) / lonstep));
|
|
|
|
ppar->xe = static_cast<size_t>(Ceil((lon2 - lonb) / lonstep));
|
|
|
|
|
|
|
|
if(ppar->xb == ppar->xe) return {nullptr, "Lonb must be not equal late"};
|
|
|
|
|
|
|
|
if(depth < 0.0 || depth > depths.back())
|
|
|
|
ppar->layer = (depth < 0.0) ? 0 : (depths.size() - 1);
|
|
|
|
else
|
|
|
|
for(size_t i = 0; i < depths.size() - 1; i++)
|
|
|
|
{
|
|
|
|
if(depth >= depths[i] && depth <= depths[i + 1])
|
|
|
|
{
|
|
|
|
ppar->layer = (depth - depths[i] <= depths[i + 1] - depth) ? i : (i + 1);
|
|
|
|
break;
|
|
|
|
}
|
|
|
|
}
|
|
|
|
}
|
|
|
|
|
|
|
|
pars.SetParameter("variable", ppar->varname);
|
|
|
|
pars.SetParameter("depth", Depth(ppar->layer));
|
|
|
|
pars.SetParameter("layer", ppar->layer);
|
|
|
|
pars.SetParameter("dataset", Title());
|
|
|
|
pars.SetParameter("lonb", Lon(ppar->xb));
|
|
|
|
pars.SetParameter("latb", Lat(ppar->yb));
|
|
|
|
pars.SetParameter("lone", Lon(ppar->xe));
|
|
|
|
pars.SetParameter("late", Lat(ppar->ye));
|
|
|
|
|
|
|
|
return {ppar.release(), ""};
|
|
|
|
}
|
|
|
|
|
|
|
|
MString Dump(const struct Parameters* ppar) const
|
|
|
|
{
|
|
|
|
// clang-format off
|
|
|
|
return
|
|
|
|
"Current settings:\n" + MString() +
|
|
|
|
" Longitudes: from " + Lon(ppar->xb) + " (" + ppar->xb + ") to "+ Lon(ppar->xe) + " (" + ppar->xe + ")\n" +
|
|
|
|
" Latitudes: from " + Lat(ppar->yb) + " (" + ppar->yb + ") to "+ Lat(ppar->ye) + " (" + ppar->ye + ")\n" +
|
|
|
|
" Depth: layer " + ppar->layer + ", depth " + Depth(ppar->layer) + " m\n";
|
|
|
|
// clang-format on
|
|
|
|
}
|
|
|
|
|
|
|
|
real Lon(size_t ix) const { return isOk() ? (-180.0 + ix / 12.0) : -1000.0; }
|
|
|
|
real Lat(size_t iy) const { return isOk() ? (-80.0 + iy / 12.0) : -1000.0; }
|
|
|
|
real Depth(size_t l) const { return isOk() ? depths[l] : -1000.0; }
|
|
|
|
time_t Timestep() const { return isOk() ? (times[1] - times[0]) : 0; }
|
|
|
|
MDateTime Time(size_t i) const
|
|
|
|
{
|
|
|
|
if(!isOk() || i >= times.size()) return MDateTime();
|
|
|
|
return times[i];
|
|
|
|
}
|
|
|
|
|
|
|
|
bool isOk() const { return nc.size() > 0; }
|
|
|
|
explicit operator bool() const { return nc.size() > 0; }
|
|
|
|
|
|
|
|
size_t NDepths() const { return depths.size(); }
|
|
|
|
size_t NTimes() const { return times.size(); }
|
|
|
|
|
|
|
|
MString Title() const
|
|
|
|
{
|
|
|
|
switch(type)
|
|
|
|
{
|
|
|
|
case(TYPE_DT): return "NEMO Delayed time, daily mean (DT)";
|
|
|
|
case(TYPE_NRT): return "NEMO Near-real time, daily mean (NRT)";
|
|
|
|
case(TYPE_NRT6): return "NEMO Near-real time, 6h resolution (NRT6)";
|
|
|
|
default: return "No title";
|
|
|
|
}
|
|
|
|
}
|
|
|
|
|
|
|
|
static real Fillval() { return Data::Fillval(); }
|
|
|
|
};
|