Browse Source

Merged HYCOM and NEMO classes.

Configuration file name changed.
URLs for data are now read from the configuration.
Variables are searched for by the standard_name attribute.
interpolate
Michael Uleysky 2 years ago
parent
commit
4dc5eda1df
  1. 390
      include/HYCOM.h
  2. 441
      include/NEMO.h
  3. 155
      include/layereddata.h
  4. 2
      src/ParseArgs.cpp
  5. 290
      src/layereddata.cpp

390
include/HYCOM.h

@ -1,26 +1,7 @@
#pragma once #pragma once
#include "DataAdapters/ncfilealt.h" #include "layereddata.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; class HYCOMData: public LayeredData
using michlib::DetGeoDomain;
using michlib::errmessage;
using michlib::Floor;
using michlib::GPL;
using michlib::int2;
using michlib::MDateTime;
using michlib::NCFileA;
using michlib::ToGeoDomain;
class HYCOMData
{ {
enum Type enum Type
{ {
@ -30,132 +11,17 @@ class HYCOMData
TYPE_FORECAST TYPE_FORECAST
}; };
struct Parameters: public BaseParameters
{
size_t xb, yb, xe, ye, layer;
MString varname;
virtual ~Parameters() override = default;
};
std::vector<NCFileA> nc;
std::vector<MString> urls;
std::vector<real> depths;
std::vector<MDateTime> times;
MString lonname, latname;
size_t nx, ny;
real lonb, latb, lone, late;
real lonstep, latstep;
Type type = TYPE_UNKNOWN; Type type = TYPE_UNKNOWN;
public: MString DataTitle() const
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_d = f.A<double>(name, "add_offset");
auto a_scale_d = f.A<double>(name, "scale_factor");
auto a_offset_f = f.A<float>(name, "add_offset");
auto a_scale_f = f.A<float>(name, "scale_factor");
if(!a_fill) return Data();
fill = a_fill;
if(a_offset_d) offset = a_offset_d;
if(a_scale_d) scale = a_scale_d;
if(a_offset_f) offset = a_offset_f;
if(a_scale_f) scale = a_scale_f;
}
auto unit = f.A<MString>(name, "units");
if(unit && (unit.Get() == "m s-1" || unit.Get() == "m/s")) 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, {lonname, p->xb, p->xe - p->xb + 1}, {latname, p->yb, p->ye - p->yb + 1}, {"time", i, 1})
: f.V<DataType>(name, {lonname, p->xb, p->xe - p->xb + 1}, {latname, 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, {lonname, p->xb}, {latname, p->yb, p->ye - p->yb + 1}, {"time", i, 1})
: f.V<DataType>(name, {lonname, p->xb}, {latname, p->yb, p->ye - p->yb + 1}, {"time", i, 1}, {"depth", p->layer, 1});
auto var2 = nodepth ? f.V<DataType>(name, {lonname, 0, p->xe + 1}, {latname, p->yb, p->ye - p->yb + 1}, {"time", i, 1})
: f.V<DataType>(name, {lonname, 0, p->xe + 1}, {latname, 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 "water_temp";
if(vname == "sal") return "salinity";
if(vname == "ssh") return "surf_el";
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; switch(type)
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(); case(TYPE_REANALYSIS): return "GOFS 3.1: 41-layer HYCOM + NCODA Global 1/12 Reanalysis";
for(const auto& v: head.Variables()) case(TYPE_HINDCAST): return "GOFS 3.1: 41-layer HYCOM + NCODA Global 1/12 Analysis Hindcast";
if(v.Name() == name) return true; case(TYPE_FORECAST): return "GOFS 3.1: 41-layer HYCOM + NCODA Global 1/12 Analysis Forecast";
default: return "No title";
} }
return false;
} }
public: public:
@ -166,251 +32,17 @@ class HYCOMData
{ {
MString dataset = args.contains("dataset") ? args.at("dataset") : "Forecast"; MString dataset = args.contains("dataset") ? args.at("dataset") : "Forecast";
nc.clear(); GPL.UsePrefix("HYCOM");
if(dataset == "Forecast") if(dataset == "Forecast")
{
urls = {"https://tds.hycom.org/thredds/dodsC/GLBy0.08/expt_93.0/FMRC/GLBy0.08_930_FMRC_best.ncd"};
type = TYPE_FORECAST; type = TYPE_FORECAST;
}
else if(dataset == "Hindcast") else if(dataset == "Hindcast")
{
urls = {"https://tds.hycom.org/thredds/dodsC/GLBy0.08/expt_93.0"};
type = TYPE_HINDCAST; type = TYPE_HINDCAST;
}
else if(dataset == "Reanalysis") else if(dataset == "Reanalysis")
{
urls.clear();
for(size_t i = 1994; i <= 2015; i++) urls.push_back("https://tds.hycom.org/thredds/dodsC/GLBv0.08/expt_53.X/data/" + MString(i));
type = TYPE_REANALYSIS; type = TYPE_REANALYSIS;
}
else else
return "Unknown dataset: " + dataset; return "Unknown dataset: " + dataset;
nc.resize(urls.size()); SetTitle(DataTitle());
for(size_t i = 0; i < urls.size(); i++) return LayeredData::Open(args, dataset);
{
nc[i].Reset(urls[i] + "#cache&noprefetch");
if(!nc[i])
{
nc.clear();
return "Can't connect to database";
}
}
auto head = nc[0].Header();
lonname = latname = "";
for(const auto& dim: head.Dimensions())
{
if(dim.Name() == "lon" || dim.Name() == "longitude")
{
lonname = dim.Name();
nx = dim.Len();
}
if(dim.Name() == "lat" || dim.Name() == "latitude")
{
latname = dim.Name();
ny = dim.Len();
}
}
if(!(lonname.Exist() && latname.Exist()))
{
nc.clear();
return "Can't find longitude/latitude";
}
auto rdepths = nc[0].VR("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 time = nc[0].VR("time");
if(!time)
{
nc.clear();
return "Can't read times";
}
MDateTime refdate;
{
auto units = nc[0].Attribute<MString>("time", "units");
if(!units)
{
nc.clear();
return "Can't read refdate";
}
MString rstr;
auto words = michlib::Split_on_words(units);
auto ci = words.begin();
if(ci != words.end()) ci++; // skip "hours"
if(ci != words.end()) ci++; // skip "since"
if(ci != words.end()) rstr = *ci; // Day
if(ci != words.end()) ci++;
if(ci != words.end()) rstr += " " + *ci; // Hours
if(!refdate.FromString(rstr))
{
nc.clear();
return "Can't parse " + rstr + " to refdate";
}
}
times.resize(time.DimLen(0));
for(size_t i = 0; i < times.size(); i++) times[i] = refdate + static_cast<time_t>(time(i)) * 3600;
auto lons = nc[0].VR(lonname);
auto lats = nc[0].VR(latname);
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() == "water_temp") vars.emplace("temp");
if(v.Name() == "salinity") vars.emplace("sal");
if(v.Name() == "surf_el") vars.emplace("ssh");
}
}
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)"};
{
auto dom = DetGeoDomain(lonb, lone);
real lon1 = ToGeoDomain(args.at("lonb").ToReal(), dom);
real lon2 = ToGeoDomain(args.at("lone").ToReal(), dom);
real lat1 = args.at("latb").ToReal();
real lat2 = 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 > ny - 1) ppar->ye = ny - 1;
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() ? (lonb + ix * lonstep) : -1000.0; }
real Lat(size_t iy) const { return isOk() ? (latb + iy * latstep) : -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_REANALYSIS): return "GOFS 3.1: 41-layer HYCOM + NCODA Global 1/12 Reanalysis";
case(TYPE_HINDCAST): return "GOFS 3.1: 41-layer HYCOM + NCODA Global 1/12 Analysis Hindcast";
case(TYPE_FORECAST): return "GOFS 3.1: 41-layer HYCOM + NCODA Global 1/12 Analysis Forecast";
default: return "No title";
}
}
static real Fillval() { return Data::Fillval(); }
}; };

441
include/NEMO.h

@ -1,26 +1,7 @@
#pragma once #pragma once
#include "DataAdapters/ncfilealt.h" #include "layereddata.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; class NEMOData: public LayeredData
using michlib::DetGeoDomain;
using michlib::errmessage;
using michlib::Floor;
using michlib::GPL;
using michlib::int2;
using michlib::MDateTime;
using michlib::NCFileA;
using michlib::ToGeoDomain;
class NEMOData
{ {
enum Type enum Type
{ {
@ -30,170 +11,17 @@ class NEMOData
TYPE_NRT6 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<MString> urls;
std::vector<real> depths;
std::vector<MDateTime> times;
MString lonname, latname;
size_t nx, ny;
real lonb, latb, lone, late;
real lonstep, latstep;
Type type = TYPE_UNKNOWN; Type type = TYPE_UNKNOWN;
class EnvVar MString DataTitle() const
{
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_d = f.A<double>(name, "add_offset");
auto a_scale_d = f.A<double>(name, "scale_factor");
auto a_offset_f = f.A<float>(name, "add_offset");
auto a_scale_f = f.A<float>(name, "scale_factor");
if(!a_fill) return Data();
fill = a_fill;
if(a_offset_d) offset = a_offset_d;
if(a_scale_d) scale = a_scale_d;
if(a_offset_f) offset = a_offset_f;
if(a_scale_f) scale = a_scale_f;
}
auto unit = f.A<MString>(name, "units");
if(unit && (unit.Get() == "m s-1" || unit.Get() == "m/s")) 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, {lonname, p->xb, p->xe - p->xb + 1}, {latname, p->yb, p->ye - p->yb + 1}, {"time", i, 1}) switch(type)
: f.V<DataType>(name, {lonname, p->xb, p->xe - p->xb + 1}, {latname, 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, {lonname, p->xb}, {latname, p->yb, p->ye - p->yb + 1}, {"time", i, 1})
: f.V<DataType>(name, {lonname, p->xb}, {latname, p->yb, p->ye - p->yb + 1}, {"time", i, 1}, {"depth", p->layer, 1});
auto var2 = nodepth ? f.V<DataType>(name, {lonname, 0, p->xe + 1}, {latname, p->yb, p->ye - p->yb + 1}, {"time", i, 1})
: f.V<DataType>(name, {lonname, 0, p->xe + 1}, {latname, 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(); case(TYPE_DT): return "NEMO Delayed time, daily mean (DT)";
for(const auto& v: head.Variables()) case(TYPE_NRT): return "NEMO Near-real time, daily mean (NRT)";
if(v.Name() == name) return true; case(TYPE_NRT6): return "NEMO Near-real time, 6h resolution (NRT6)";
default: return "No title";
} }
return false;
} }
public: public:
@ -202,266 +30,19 @@ class NEMOData
// TODO: RetVal // TODO: RetVal
MString Open(const CLArgs& args) MString Open(const CLArgs& args)
{ {
MString oldprefix = GPL.UsePrefix("AVISO");
MString dataset = args.contains("dataset") ? args.at("dataset") : "DT"; 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);
GPL.UsePrefix("NEMO");
if(dataset == "DT") if(dataset == "DT")
{
urls = {"https://" + cred + "@my.cmems-du.eu/thredds/dodsC/cmems_mod_glo_phy_my_0.083_P1D-m"};
type = TYPE_DT; type = TYPE_DT;
}
else if(dataset == "NRT") 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; type = TYPE_NRT;
}
else if(dataset == "NRT6") 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; type = TYPE_NRT6;
}
else else
return "Unknown dataset: " + dataset; return "Unknown dataset: " + dataset;
nc.resize(urls.size()); SetTitle(DataTitle());
for(size_t i = 0; i < urls.size(); i++) return LayeredData::Open(args, dataset);
{
nc[i].Reset(urls[i] + "#cache&noprefetch");
if(!nc[i])
{
nc.clear();
return "Can't connect to database";
}
}
auto head = nc[0].Header();
lonname = latname = "";
for(const auto& dim: head.Dimensions())
{
if(dim.Name() == "lon" || dim.Name() == "longitude")
{
lonname = dim.Name();
nx = dim.Len();
}
if(dim.Name() == "lat" || dim.Name() == "latitude")
{
latname = dim.Name();
ny = dim.Len();
}
}
if(!(lonname.Exist() && latname.Exist()))
{
nc.clear();
return "Can't find longitude/latitude";
}
auto rdepths = nc[0].VR("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 time = nc[0].VR("time");
if(!time)
{
nc.clear();
return "Can't read times";
}
MDateTime refdate;
{
auto units = nc[0].Attribute<MString>("time", "units");
if(!units)
{
nc.clear();
return "Can't read refdate";
}
MString rstr;
auto words = michlib::Split_on_words(units);
auto ci = words.begin();
if(ci != words.end()) ci++; // skip "hours"
if(ci != words.end()) ci++; // skip "since"
if(ci != words.end()) rstr = *ci; // Day
if(ci != words.end()) ci++;
if(ci != words.end()) rstr += " " + *ci; // Hours
if(!refdate.FromString(rstr))
{
nc.clear();
return "Can't parse " + rstr + " to refdate";
}
}
times.resize(time.DimLen(0));
for(size_t i = 0; i < times.size(); i++) times[i] = refdate + static_cast<time_t>(time(i)) * 3600;
auto lons = nc[0].VR(lonname);
auto lats = nc[0].VR(latname);
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)"};
{
auto dom = DetGeoDomain(lonb, lone);
real lon1 = ToGeoDomain(args.at("lonb").ToReal(), dom);
real lon2 = ToGeoDomain(args.at("lone").ToReal(), dom);
real lat1 = args.at("latb").ToReal();
real lat2 = 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 > ny - 1) ppar->ye = ny - 1;
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() ? (lonb + ix * lonstep) : -1000.0; }
real Lat(size_t iy) const { return isOk() ? (latb + iy * latstep) : -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(); }
}; };

155
include/layereddata.h

@ -0,0 +1,155 @@
#pragma once
#include "DataAdapters/ncfilealt.h"
#include "ParseArgs.h"
#include "mdatetime.h"
#include "simple2ddata.h"
#include <memory>
#include <set>
using michlib::Ceil;
using michlib::DetGeoDomain;
using michlib::Floor;
using michlib::GPL;
using michlib::int2;
using michlib::MDateTime;
using michlib::MString;
using michlib::NCFileA;
using michlib::ToGeoDomain;
class LayeredData
{
public:
using Data = Simple2DData;
private:
std::vector<NCFileA> nc;
std::vector<MString> urls;
std::vector<real> depths;
std::vector<MDateTime> times;
MString lonname, latname;
size_t nx, ny;
real lonb, latb, lone, late;
real lonstep, latstep;
MString title;
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;
protected:
struct Parameters: public BaseParameters
{
size_t xb, yb, xe, ye, layer;
MString varname;
virtual ~Parameters() override = default;
};
// TODO: RetVal
MString Open(const CLArgs& args, const MString& dataset);
void SetTitle(const MString& newtitle) { title = newtitle; }
public:
MString Info() const;
std::pair<const BaseParameters*, MString> Parameters(michlib_internal::ParameterListEx& pars, const CLArgs& args) const;
Data Read(const BaseParameters* ip, size_t i) const;
bool isOk() const { return nc.size() > 0; }
explicit operator bool() const { return nc.size() > 0; }
real Depth(size_t l) const { return isOk() ? depths[l] : -1000.0; }
real Lon(size_t ix) const { return isOk() ? (lonb + ix * lonstep) : -1000.0; }
real Lat(size_t iy) const { return isOk() ? (latb + iy * latstep) : -1000.0; }
size_t NDepths() const { return depths.size(); }
size_t NTimes() const { return times.size(); }
MDateTime Time(size_t i) const
{
if(!isOk() || i >= times.size()) return MDateTime();
return times[i];
}
time_t Timestep() const { return isOk() ? (times[1] - times[0]) : 0; }
MString Title() const { return title; }
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
}
private:
template<class DataType> Data ReadVarRaw(const NCFileA& f, const MString& name, size_t i, bool nodepth, const struct Parameters* p) const;
std::pair<MString, size_t> VarNameLoc(const MString vname) const
{
for(size_t i = 0; i < nc.size(); i++)
{
auto head = nc[i].Header();
for(const auto& v: head.Variables())
{
auto stname = nc[i].A<MString>(v.Name(), "standard_name");
if(!stname) continue;
if(StName2Name(stname) == vname) return {v.Name(), i};
}
}
return {"", 0};
}
static MString StName2Name(const MString& stname)
{
if(stname == "sea_water_potential_temperature") return "temp";
if(stname == "sea_water_temperature") return "temp";
if(stname == "sea_water_salinity") return "sal";
if(stname == "ocean_mixed_layer_thickness_defined_by_sigma_theta") return "mld";
if(stname == "sea_surface_height_above_geoid") return "ssh";
if(stname == "sea_surface_elevation") return "ssh";
if(stname == "upward_sea_water_velocity") return "w";
return "";
}
};

2
src/ParseArgs.cpp

@ -28,7 +28,7 @@ CLArgs ParseArgs(int argc, char** argv)
} }
SList sl; SList sl;
michlib_internal::ParseParameterFile("/etc/michlibconf", sl, 0); michlib_internal::ParseParameterFile("/etc/odm.conf", sl, 0);
return out; return out;
} }

290
src/layereddata.cpp

@ -0,0 +1,290 @@
#define MICHLIB_NOSOURCE
#include "layereddata.h"
MString LayeredData::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())
{
auto ret = f.A<MString>(v.Name(), "standard_name");
if(!ret) continue;
vars.emplace(StName2Name(ret));
}
}
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
}
MString LayeredData::Open(const CLArgs& args, const MString& dataset)
{
nc.clear();
MString proxyurl = GPL.ParameterSValue("USEPROXY", "");
if(proxyurl.Exist()) proxy.Activate("all_proxy", proxyurl);
urls.clear();
size_t i = 1;
while(true)
{
MString url = GPL.ParameterSValue(dataset + "_URL" + i, "");
if(url.Exist())
urls.emplace_back(std::move(url));
else
break;
i++;
}
if(urls.size() == 0) return "No urls for dataset " + dataset + " specified in config";
nc.resize(urls.size());
for(size_t i = 0; i < urls.size(); i++)
{
nc[i].Reset(urls[i] + "#cache&noprefetch");
if(!nc[i])
{
nc.clear();
return "Can't connect to database";
}
}
auto head = nc[0].Header();
lonname = latname = "";
for(const auto& dim: head.Dimensions())
{
if(dim.Name() == "lon" || dim.Name() == "longitude")
{
lonname = dim.Name();
nx = dim.Len();
}
if(dim.Name() == "lat" || dim.Name() == "latitude")
{
latname = dim.Name();
ny = dim.Len();
}
}
if(!(lonname.Exist() && latname.Exist()))
{
nc.clear();
return "Can't find longitude/latitude";
}
auto rdepths = nc[0].VR("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 time = nc[0].VR("time");
if(!time)
{
nc.clear();
return "Can't read times";
}
MDateTime refdate;
{
auto units = nc[0].Attribute<MString>("time", "units");
if(!units)
{
nc.clear();
return "Can't read refdate";
}
MString rstr;
auto words = michlib::Split_on_words(units);
auto ci = words.begin();
if(ci != words.end()) ci++; // skip "hours"
if(ci != words.end()) ci++; // skip "since"
if(ci != words.end()) rstr = *ci; // Day
if(ci != words.end()) ci++;
if(ci != words.end()) rstr += " " + *ci; // Hours
if(!refdate.FromString(rstr))
{
nc.clear();
return "Can't parse " + rstr + " to refdate";
}
}
times.resize(time.DimLen(0));
for(size_t i = 0; i < times.size(); i++) times[i] = refdate + static_cast<time_t>(time(i)) * 3600;
auto lons = nc[0].VR(lonname);
auto lats = nc[0].VR(latname);
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 "";
}
std::pair<const BaseParameters*, MString> LayeredData::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(!VarNameLoc(ppar->varname).first.Exist()) 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)"};
{
auto dom = DetGeoDomain(lonb, lone);
real lon1 = ToGeoDomain(args.at("lonb").ToReal(), dom);
real lon2 = ToGeoDomain(args.at("lone").ToReal(), dom);
real lat1 = args.at("latb").ToReal();
real lat2 = 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 > ny - 1) ppar->ye = ny - 1;
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(), ""};
}
LayeredData::Data LayeredData::Read(const BaseParameters* ip, size_t i) const
{
if(!isOk()) return Data();
bool nodepth = false;
auto p = dynamic_cast<const struct Parameters*>(ip);
auto [name, id] = VarNameLoc(p->varname);
auto head = nc[id].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>(nc[id], name, i, nodepth, p);
if(v.Type().Id() == NC_FLOAT) return ReadVarRaw<float>(nc[id], name, i, nodepth, p);
}
return Data();
}
template<class DataType> LayeredData::Data LayeredData::ReadVarRaw(const NCFileA& f, const MString& name, size_t i, bool nodepth, const struct LayeredData::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_d = f.A<double>(name, "add_offset");
auto a_scale_d = f.A<double>(name, "scale_factor");
auto a_offset_f = f.A<float>(name, "add_offset");
auto a_scale_f = f.A<float>(name, "scale_factor");
if(!a_fill) return Data();
fill = a_fill;
if(a_offset_d) offset = a_offset_d;
if(a_scale_d) scale = a_scale_d;
if(a_offset_f) offset = a_offset_f;
if(a_scale_f) scale = a_scale_f;
}
auto unit = f.A<MString>(name, "units");
if(unit && (unit.Get() == "m s-1" || unit.Get() == "m/s")) 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, {lonname, p->xb, p->xe - p->xb + 1}, {latname, p->yb, p->ye - p->yb + 1}, {"time", i, 1})
: f.V<DataType>(name, {lonname, p->xb, p->xe - p->xb + 1}, {latname, 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, {lonname, p->xb}, {latname, p->yb, p->ye - p->yb + 1}, {"time", i, 1})
: f.V<DataType>(name, {lonname, p->xb}, {latname, p->yb, p->ye - p->yb + 1}, {"time", i, 1}, {"depth", p->layer, 1});
auto var2 = nodepth ? f.V<DataType>(name, {lonname, 0, p->xe + 1}, {latname, p->yb, p->ye - p->yb + 1}, {"time", i, 1})
: f.V<DataType>(name, {lonname, 0, p->xe + 1}, {latname, 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;
}
Loading…
Cancel
Save