michael_ul
/
odm
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity
Browse Source

Added support for time-splitted data (HYCOM Reanalysis)

interpolate
Michael Uleysky 2 years ago
parent
3649579610
commit
5359665083
2 changed files with 136 additions and 81 deletions
Whitespace
Show all changes Ignore whitespace when comparing lines Ignore changes in amount of whitespace Ignore changes in whitespace at EOL
Split View
Diff Options
Show Stats Download Patch File Download Diff File
  1. 93
      include/layereddata.h
  2. 124
      src/layereddata.cpp

93
include/layereddata.h
Unescape View File

@ -3,6 +3,7 @@
#include "ParseArgs.h"
#include "mdatetime.h"
#include "simple2ddata.h"
#include <algorithm>
#include <memory>
#include <set>
@ -22,8 +23,67 @@ class LayeredData
using Data = Simple2DData;
private:
std::vector<NCFileA> nc;
std::vector<MString> urls;
class NC
{
MString url;
NCFileA nc;
std::vector<MDateTime> times;
public:
NC(MString&& newurl): url(std::move(newurl)) { nc.Reset(url + "#cache&noprefetch"); }
MString ReadTimes()
{
if(!nc) return "File not open";
auto time = nc.VR("time");
if(!time) return "Can't read times";
MDateTime refdate;
{
auto units = nc.Attribute<MString>("time", "units");
if(!units) 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)) 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;
return "";
}
const NCFileA* operator->() const { return &nc; }
explicit operator bool() const { return nc; }
MDateTime Begin() const { return times.front(); }
MDateTime End() const { return times.back(); }
const std::vector<MDateTime>& Times() const { return times; }
size_t Index(MDateTime tm) const
{
if(tm < Begin() || tm > End()) return 0;
size_t b = 0, e = times.size() - 1;
if(tm == times[b]) return b + 1;
if(tm == times[e]) return e + 1;
while(e - b > 1)
{
size_t c = (e + b) / 2;
if(tm == times[c]) return c + 1;
if(tm > times[c])
b = c;
else
e = c;
}
return 0;
}
};
std::vector<NC> nc;
std::vector<real> depths;
std::vector<MDateTime> times;
MString lonname, latname;
@ -124,21 +184,38 @@ class LayeredData
}
private:
template<class DataType> Data ReadVarRaw(const NCFileA& f, const MString& name, size_t i, bool nodepth, const struct Parameters* p) const;
template<class DataType> Data ReadVarRaw(const NC& f, const MString& name, size_t i, bool nodepth, const struct Parameters* p) const;
bool HaveVar(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 true;
}
}
return false;
}
std::pair<MString, size_t> VarNameLoc(const MString vname) const
std::tuple<MString, size_t, size_t> VarNameLoc(const MString vname, MDateTime tm) const
{
for(size_t i = 0; i < nc.size(); i++)
{
auto head = nc[i].Header();
auto head = nc[i]->Header();
auto tind = nc[i].Index(tm);
if(tind == 0) continue;
for(const auto& v: head.Variables())
{
auto stname = nc[i].A<MString>(v.Name(), "standard_name");
auto stname = nc[i]->A<MString>(v.Name(), "standard_name");
if(!stname) continue;
if(StName2Name(stname) == vname) return {v.Name(), i};
if(StName2Name(stname) == vname) return {v.Name(), i, tind - 1};
}
}
return {"", 0};
return {"", 0, 0};
}
static MString StName2Name(const MString& stname)

124
src/layereddata.cpp
Unescape View File

@ -10,12 +10,12 @@ MString LayeredData::Info() const
std::set<MString> vars;
for(const auto& f: nc)
{
auto head = f.Header();
auto head = f->Header();
for(const auto& v: head.Variables())
{
auto ret = f.A<MString>(v.Name(), "standard_name");
auto ret = f->A<MString>(v.Name(), "standard_name");
if(!ret) continue;
vars.emplace(StName2Name(ret));
if(StName2Name(ret).Exist()) vars.emplace(StName2Name(ret));
}
}
MString svars;
@ -48,31 +48,34 @@ MString LayeredData::Open(const MString& dataset)
MString proxyurl = GPL.ParameterSValue("USEPROXY", "");
if(proxyurl.Exist()) proxy.Activate("all_proxy", proxyurl);
urls.clear();
nc.clear();
size_t i = 1;
while(true)
{
MString url = GPL.ParameterSValue(dataset + "_URL" + i, "");
if(url.Exist())
urls.emplace_back(std::move(url));
{
//michlib::message("Open "+url);
nc.emplace_back(std::move(url));
if(!nc.back())
{
nc.clear();
return "Can't connect to url " + url;
}
MString ret = nc.back().ReadTimes();
if(ret.Exist())
{
nc.clear();
return "Can't connect to url " + url + ": " + ret;
}
}
else
break;
i++;
}
if(urls.size() == 0) return "No urls for dataset " + dataset + " specified in config";
if(nc.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();
auto head = nc[0]->Header();
lonname = latname = "";
for(const auto& dim: head.Dimensions())
{
@ -93,7 +96,15 @@ MString LayeredData::Open(const MString& dataset)
return "Can't find longitude/latitude";
}
auto rdepths = nc[0].VR("depth");
// Read times
{
for(const auto& f: nc) times.insert(times.end(), f.Times().begin(), f.Times().end());
std::sort(times.begin(), times.end());
auto last = std::unique(times.begin(), times.end());
times.erase(last, times.end());
}
auto rdepths = nc[0]->VR("depth");
if(!rdepths)
{
nc.clear();
@ -102,41 +113,8 @@ MString LayeredData::Open(const MString& dataset)
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);
auto lons = nc[0]->VR(lonname);
auto lats = nc[0]->VR(latname);
if(!(lons && lats))
{
nc.clear();
@ -158,7 +136,7 @@ std::pair<const BaseParameters*, MString> LayeredData::Parameters(michlib_intern
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(!HaveVar(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);
@ -212,32 +190,32 @@ 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();
auto p = dynamic_cast<const struct Parameters*>(ip);
auto [name, id, tid] = VarNameLoc(p->varname, times[i]);
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);
if(v.Type().Id() == NC_SHORT) return ReadVarRaw<int2>(nc[id], name, tid, nodepth, p);
if(v.Type().Id() == NC_FLOAT) return ReadVarRaw<float>(nc[id], name, tid, 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
template<class DataType> LayeredData::Data LayeredData::ReadVarRaw(const NC& 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");
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;
@ -246,15 +224,15 @@ template<class DataType> LayeredData::Data LayeredData::ReadVarRaw(const NCFileA
if(a_scale_f) scale = a_scale_f;
}
auto unit = f.A<MString>(name, "units");
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});
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();
@ -267,10 +245,10 @@ template<class DataType> LayeredData::Data LayeredData::ReadVarRaw(const NCFileA
}
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});
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++)

Write Preview
Loading…
Cancel
Save