michael_ul
/
odm
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity
You can not select more than 25 topics Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
8 Commits
3 Branches
0 Tags
1.7 MiB
 
 
 
Tree: 4f6782ea09
Branches Tags
${ item.name }
Create tag ${ searchTerm }
Create branch ${ searchTerm }
from '4f6782ea09'
${ noResults }
odm/include/NEMO.h

372 lines
12 KiB
Raw Blame History

#include "DataAdapters/ncfilealt.h"
#include "mdatetime.h"
#include "simple2ddata.h"
#include "specfunc.h"
#include "vartype.h"
#include <regex.h>
#include <stdlib.h>
#include <sys/types.h>
#if !defined(M__NEMO)
#define M__NEMO
using michlib::Ceil;
using michlib::Floor;
using michlib::int2;
using michlib::MDateTime;
using michlib::NCFileA;
using michlib::ToGeoDomainNeg;
class NEMOData
{
enum Type
{
TYPE_UNKNOWN,
TYPE_DT,
TYPE_NRT,
TYPE_NRT6
};
std::vector<NCFileA> nc;
size_t xb, yb, xe, ye, layer;
std::vector<real> depths;
std::vector<MDateTime> times;
size_t nx;
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
{
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((xb < xe) ? (xe - xb + 1) : (nx + xe - xb + 1), ye - yb + 1, Lonb(), Latb());
if(xb < xe)
{
auto var = nodepth ? f.V<DataType>(name, {"longitude", xb, xe - xb + 1}, {"latitude", yb, ye - yb + 1}, {"time", i, 1})
: f.V<DataType>(name, {"longitude", xb, xe - xb + 1}, {"latitude", yb, ye - yb + 1}, {"time", i, 1}, {"depth", 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", xb}, {"latitude", yb, ye - yb + 1}, {"time", i, 1})
: f.V<DataType>(name, {"longitude", xb}, {"latitude", yb, ye - yb + 1}, {"time", i, 1}, {"depth", layer, 1});
auto var2 = nodepth ? f.V<DataType>(name, {"longitude", 0, xe + 1}, {"latitude", yb, ye - yb + 1}, {"time", i, 1})
: f.V<DataType>(name, {"longitude", 0, xe + 1}, {"latitude", yb, ye - yb + 1}, {"time", i, 1}, {"depth", 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;
}
public:
Data ReadVar(const MString& name, size_t i) const
{
if(!isOk()) return Data();
bool nodepth = false;
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);
if(v.Type().Id() == NC_FLOAT) return ReadVarRaw<float>(f, name, i, nodepth);
}
}
return Data();
}
Data ReadVar(const MString& name, const std::vector<size_t>& tindex) const
{
Data out;
if(tindex.size() == 0 || !isOk()) return out;
std::vector<size_t> count;
for(size_t i = 0; i < tindex.size(); i++)
{
Data dat = ReadVar(name, tindex[i]);
if(!dat) return Data();
out.Add(dat, count);
}
out.Div(count);
return out;
}
public:
NEMOData() = default;
// TODO: RetVal
bool Open(const MString& stype, const MString& cred, const MString& proxyurl = "")
{
nc.clear();
if(proxyurl.Exist()) proxy.Activate("all_proxy", proxyurl);
if(stype == "DT")
{
NCFileA newnc;
newnc.Reset("https://" + cred + "@my.cmems-du.eu/thredds/dodsC/cmems_mod_glo_phy_my_0.083_P1D-m");
if(!newnc) return false;
nc.push_back(std::move(newnc));
type = TYPE_DT;
}
if(stype == "NRT")
{
std::vector<MString> 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"};
for(const auto& url: urls)
{
NCFileA newnc;
newnc.Reset(url);
if(!newnc)
{
nc.clear();
return false;
}
nc.push_back(std::move(newnc));
}
type = TYPE_NRT;
}
if(stype == "NRT6")
{
std::vector<MString> 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"};
for(const auto& url: urls)
{
NCFileA newnc;
newnc.Reset(url);
if(!newnc)
{
nc.clear();
return false;
}
nc.push_back(std::move(newnc));
}
type = TYPE_NRT6;
}
auto rdepths = nc[0].V<float>("depth");
if(!rdepths)
{
nc.clear();
return false;
}
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 false;
}
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;
return true;
}
MString Info() const
{
MString d;
for(size_t i = 0; i < NDepths(); i++) d += MString(" ") + "(" + i + " " + Depth(i) + ")";
// 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" +
" Depths:" + d;
// clang-format on
}
MString Dump() const
{
// clang-format off
return
"Current settings:\n" + MString() +
" Longitudes: from " + Lonb() + " (" + xb + ") to "+ Lone() + " (" + xe + ")\n" +
" Latitudes: from " + Latb() + " (" + yb + ") to "+ Late() + " (" + ye + ")\n" +
" Depth: layer " + layer + ", depth " + Depth() + " m\n";
// clang-format on
}
bool SetRegion(real lonbin, real lonein, real latb, real late, real depth)
{
if(!isOk()) return false;
real lonb = ToGeoDomainNeg(lonbin), lone = ToGeoDomainNeg(lonein);
size_t xb_, xe_, yb_, ye_, layer_ = 0;
yb_ = static_cast<size_t>(Floor((latb + 80.0) * 12.0));
ye_ = static_cast<size_t>(Ceil((late + 80.0) * 12.0));
if(ye_ > 2040) ye_ = 2040;
if(yb_ >= ye_) return false;
xb_ = static_cast<size_t>(Floor((lonb + 180.0) * 12.0));
xe_ = static_cast<size_t>(Ceil((lone + 180.0) * 12.0));
if(xb_ == xe_) return false;
if(depth < 0.0 || depth > depths.back())
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])
{
layer_ = (depth - depths[i] <= depths[i + 1] - depth) ? i : (i + 1);
break;
}
}
xb = xb_;
xe = xe_;
yb = yb_;
ye = ye_;
layer = layer_;
return true;
}
real Lonb() const { return isOk() ? (-180.0 + xb / 12.0) : -1000.0; }
real Lone() const { return isOk() ? (-180.0 + xe / 12.0) : -1000.0; }
real Latb() const { return isOk() ? (-80.0 + yb / 12.0) : -1000.0; }
real Late() const { return isOk() ? (-80.0 + ye / 12.0) : -1000.0; }
size_t Xb() const { return isOk() ? xb : 0; }
size_t Xe() const { return isOk() ? xe : 0; }
size_t Yb() const { return isOk() ? yb : 0; }
size_t Ye() const { return isOk() ? ye : 0; }
size_t Nx() const { return isOk() ? ((xb < xe) ? (xe - xb + 1) : (nx + xe - xb + 1)) : 0; }
size_t Ny() const { return isOk() ? (ye - yb + 1) : 0; }
real Depth(size_t l) const { return isOk() ? depths[l] : -1000.0; }
real Depth() const { return Depth(layer); }
size_t Layer() const { return layer; }
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(); }
template<vartype::Vartype vt> Data Read(size_t it) const = delete;
template<vartype::Vartype vt> Data Read(const std::vector<size_t>& tindex) const = delete;
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(); }
};
template<> inline NEMOData::Data NEMOData::Read<vartype::Vartype::U>(size_t it) const { return ReadVar("uo", it); }
template<> inline NEMOData::Data NEMOData::Read<vartype::Vartype::V>(size_t it) const { return ReadVar("vo", it); }
template<> inline NEMOData::Data NEMOData::Read<vartype::Vartype::TEMP>(size_t it) const { return ReadVar("thetao", it); }
template<> inline NEMOData::Data NEMOData::Read<vartype::Vartype::SAL>(size_t it) const { return ReadVar("so", it); }
template<> inline NEMOData::Data NEMOData::Read<vartype::Vartype::MLD>(size_t it) const { return ReadVar("mlotst", it); }
template<> inline NEMOData::Data NEMOData::Read<vartype::Vartype::SSH>(size_t it) const { return ReadVar("zos", it); }
template<> inline NEMOData::Data NEMOData::Read<vartype::Vartype::W>(size_t it) const { return ReadVar("wo", it); }
template<> inline NEMOData::Data NEMOData::Read<vartype::Vartype::U>(const std::vector<size_t>& tindex) const { return ReadVar("uo", tindex); }
template<> inline NEMOData::Data NEMOData::Read<vartype::Vartype::V>(const std::vector<size_t>& tindex) const { return ReadVar("vo", tindex); }
template<> inline NEMOData::Data NEMOData::Read<vartype::Vartype::TEMP>(const std::vector<size_t>& tindex) const { return ReadVar("thetao", tindex); }
template<> inline NEMOData::Data NEMOData::Read<vartype::Vartype::SAL>(const std::vector<size_t>& tindex) const { return ReadVar("so", tindex); }
template<> inline NEMOData::Data NEMOData::Read<vartype::Vartype::MLD>(const std::vector<size_t>& tindex) const { return ReadVar("mlotst", tindex); }
template<> inline NEMOData::Data NEMOData::Read<vartype::Vartype::SSH>(const std::vector<size_t>& tindex) const { return ReadVar("zos", tindex); }
template<> inline NEMOData::Data NEMOData::Read<vartype::Vartype::W>(const std::vector<size_t>& tindex) const { return ReadVar("wo", tindex); }
#endif