Browse Source

Added source AVISOLOCAL and made a some code rearrangement.

interpolate
Michael Uleysky 2 years ago
parent
commit
712be94466
  1. 71
      include/AVISOLOCAL.h
  2. 72
      include/layereddata.h
  3. 40
      include/ncfuncs.h
  4. 10
      include/odm.h
  5. 3
      include/simple2ddata.h
  6. 258
      src/AVISOLOCAL.cpp
  7. 125
      src/layereddata.cpp
  8. 148
      src/ncfuncs.cpp
  9. 11
      src/odm.cpp

71
include/AVISOLOCAL.h

@ -0,0 +1,71 @@
#pragma once
#include "GPL.h"
#include "ParseArgs.h"
#include "mregex.h"
#include "ncfuncs.h"
#include "simple2ddata.h"
#include "uvdata.h"
#include <dirent.h>
#include <utility>
using michlib::Ceil;
using michlib::DetGeoDomain;
using michlib::Floor;
using michlib::GPL;
using michlib::int2;
using michlib::int4;
using michlib::MString;
using michlib::real;
using michlib::RegExp;
class AVISOLOCALData: public NCFuncs
{
std::vector<MDateTime> times;
MString datapath;
struct Parameters: public BaseParameters
{
real lonb, latb, lone, late;
virtual ~Parameters() override = default;
};
static MString Title() { return "AVISO local mirror"; }
public:
using Data = Simple2DData;
AVISOLOCALData() = default;
MString Info() const;
// TODO: RetVal
MString Open(const CLArgs& args);
std::pair<const BaseParameters*, MString> Parameters(michlib_internal::ParameterListEx& pars, const CLArgs& args) const;
bool isOk() const { return times.size() > 0; }
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; }
explicit operator bool() const { return times.size() > 0; }
bool CheckVar(const MString& vname) const
{
NCFileA nc;
nc.Reset(datapath + "/uv-" + times[0].ToString() + ".nc");
return NCFuncs::CheckVar(vname, [&nc = std::as_const(nc)](const MString& vn) { return HaveVar(nc, vn); });
}
Data Read(const MString& vname, const BaseParameters* ip, size_t i) const;
template<class DataType> Data ReadVarRaw(const NCFileA& nc, const MString& name, const struct Parameters* p) const;
UVData ReadUV(const BaseParameters* ip, size_t i) const;
};

72
include/layereddata.h

@ -3,6 +3,7 @@
#include "ParseArgs.h"
#include "gsw.h"
#include "mdatetime.h"
#include "ncfuncs.h"
#include "simple2ddata.h"
#include "uvdata.h"
#include <algorithm>
@ -14,12 +15,13 @@ using michlib::DetGeoDomain;
using michlib::Floor;
using michlib::GPL;
using michlib::int2;
using michlib::int4;
using michlib::MDateTime;
using michlib::MString;
using michlib::NCFileA;
using michlib::ToGeoDomain;
class LayeredData
class LayeredData: public NCFuncs
{
public:
using Data = Simple2DData;
@ -34,32 +36,22 @@ class LayeredData
public:
NC(MString&& newurl): url(std::move(newurl)) { nc.Reset(url + "#cache&noprefetch"); }
MString ReadTimes()
MString ReadTimes(const MString& tname)
{
if(!nc) return "File not open";
auto time = nc.VR("time");
auto time = nc.VR(tname);
if(!time) return "Can't read times";
MDateTime refdate;
time_t step = 0;
{
auto units = nc.Attribute<MString>("time", "units");
auto units = nc.Attribute<MString>(tname, "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())
{
if(*ci == "hours") step = 3600;
if(*ci == "days") step = 3600 * 24;
ci++;
}
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";
auto [rd, st, suc] = Refdate(units);
if(!suc) return "Can't parse " + units.Get() + " to refdate";
if(st == 0) return "Can't get timestep from string " + units.Get();
refdate = rd;
step = st;
}
times.resize(time.DimLen(0));
@ -72,6 +64,8 @@ class LayeredData
MDateTime Begin() const { return times.front(); }
MDateTime End() const { return times.back(); }
const NCFileA& Get() const { return nc; }
const std::vector<MDateTime>& Times() const { return times; }
size_t Index(MDateTime tm) const
{
@ -94,8 +88,7 @@ class LayeredData
std::vector<NC> nc;
std::vector<real> depths;
std::vector<MDateTime> times;
MString lonname, latname;
size_t nx, ny;
struct CoordNames dname;
real lonb, latb, lone, late;
real lonstep, latstep;
MString title;
@ -194,16 +187,7 @@ class LayeredData
bool CheckVar(const MString& vname) const
{
if(!HaveVar(vname))
{
bool varexist = false;
if(vname == "temp" && HaveVar("ptemp") && HaveVar("sal")) varexist = true;
if(vname == "ptemp" && HaveVar("temp") && HaveVar("sal")) varexist = true;
if(vname == "pdens" && (HaveVar("ptemp") || HaveVar("temp")) && HaveVar("sal")) varexist = true;
if((vname == "U" || vname == "U2") && HaveVar("u") && HaveVar("v")) varexist = true;
if(!varexist) return false;
}
return true;
return NCFuncs::CheckVar(vname, [this](const MString& vn) { return HaveVar(vn); });
}
private:
@ -212,15 +196,7 @@ class LayeredData
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;
}
}
if(NCFuncs::HaveVar(nc[i].Get(), vname)) return true;
return false;
}
@ -240,20 +216,4 @@ class LayeredData
}
return {"", 0, 0};
}
static MString StName2Name(const MString& stname)
{
if(stname == "sea_water_potential_temperature") return "ptemp";
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 == "eastward_sea_water_velocity") return "u";
if(stname == "northward_sea_water_velocity") return "v";
if(stname == "surface_geostrophic_eastward_sea_water_velocity") return "u";
if(stname == "surface_geostrophic_northward_sea_water_velocity") return "v";
if(stname == "upward_sea_water_velocity") return "w";
return "";
}
};

40
include/ncfuncs.h

@ -0,0 +1,40 @@
#pragma once
#include "DataAdapters/ncfilealt.h"
#include "StringFunctions.h"
#include "mdatetime.h"
#include <set>
#include <tuple>
using michlib::MDateTime;
using michlib::MString;
using michlib::NCFileA;
class NCFuncs
{
public:
struct CoordNames
{
MString lonname, latname, depthname, timename;
size_t nx, ny, nz, nt;
};
static MString StName2Name(const MString& stname);
static std::tuple<MDateTime, time_t, bool> Refdate(const MString& refdate);
static void GetVars(const NCFileA& nc, std::set<MString>& vars);
static CoordNames GetCNames(const NCFileA& nc);
static CoordNames GetDNames(const NCFileA& nc);
static bool HaveVar(const NCFileA& nc, const MString& vname);
template<class HV> static bool CheckVar(const MString& vname, HV hv)
{
if(!hv(vname))
{
bool varexist = false;
if(vname == "temp" && hv("ptemp") && hv("sal")) varexist = true;
if(vname == "ptemp" && hv("temp") && hv("sal")) varexist = true;
if(vname == "pdens" && (hv("ptemp") || hv("temp")) && hv("sal")) varexist = true;
if((vname == "U" || vname == "U2") && hv("u") && hv("v")) varexist = true;
if(!varexist) return false;
}
return true;
}
};

10
include/odm.h

@ -1,5 +1,6 @@
#pragma once
#include "AVISO.h"
#include "AVISOLOCAL.h"
#include "BFileW.h"
#include "HYCOM.h"
#include "NEMO.h"
@ -11,7 +12,7 @@ using michlib::errmessage;
using michlib::GPL;
using michlib::message;
using DataVariants = std::variant<NEMOData, HYCOMData, AVISOData>;
using DataVariants = std::variant<NEMOData, HYCOMData, AVISOData, AVISOLOCALData>;
template<class T>
concept InfoSupported = requires(T t) {
@ -154,6 +155,13 @@ class Data: public DataVariants
if(res.Exist()) return "Can't open source " + src + ":\n" + res;
*this = Data(std::move(data));
}
else if(src == "AVISOLOCAL")
{
AVISOLOCALData data;
auto res = data.Open(args);
if(res.Exist()) return "Can't open source " + src + ":\n" + res;
*this = Data(std::move(data));
}
else if(src == "HYCOM")
{
HYCOMData data;

3
include/simple2ddata.h

@ -43,4 +43,7 @@ class Simple2DData: public BaseData
real Lon(size_t i) const { return Lon(i % nx, i / nx); }
real Lat(size_t i) const { return Lat(i % nx, i / nx); }
real XStep() const { return xstep; }
real YStep() const { return ystep; }
};

258
src/AVISOLOCAL.cpp

@ -0,0 +1,258 @@
#define MICHLIB_NOSOURCE
#include "AVISOLOCAL.h"
MString AVISOLOCALData::Info() const
{
if(!isOk()) return "";
NCFileA nc;
struct CoordNames cn, dn;
std::set<MString> vars;
nc.Reset(datapath + "/uv-" + times[0].ToString() + ".nc");
if(!nc) return "Can't open file " + datapath + "/uv-" + times[0].ToString() + ".nc";
GetVars(nc, vars);
dn = GetDNames(nc);
cn = GetCNames(nc);
MString svars;
{
bool first = true;
for(const auto& v: vars)
{
svars += (first ? "" : ", ") + v;
first = false;
}
}
auto lons = nc.VR(cn.lonname);
auto lats = nc.VR(cn.latname);
if(!(lons && lats)) return "Can't get longitudes/latitudes";
real lonb = lons(0);
real latb = lats(0);
real lone = lons(dn.nx - 1);
real late = lats(dn.ny - 1);
real lonstep = (lone - lonb) / (dn.nx - 1);
real latstep = (late - latb) / (dn.ny - 1);
// 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: " + dn.nx + "x" + dn.ny + " (" + lonstep + " x " + latstep + ")\n" +
" Supported variables: " + svars;
// clang-format on
}
MString AVISOLOCALData::Open(const CLArgs& args)
{
GPL.UsePrefix("AVISOLOCAL");
datapath = GPL.ParameterSValue("Datapath", "");
RegExp regex("uv-([0-9]{4}-[0-9]{2}-[0-9]{2}).nc");
regex.Compile();
DIR* dir = opendir(datapath.Buf());
struct dirent* de;
if(nullptr == dir) return "Can't open directory " + datapath;
while((de = readdir(dir)))
{
if(!regex.Match(de->d_name)) continue;
times.emplace_back(MString(de->d_name + regex.Off(1), regex.Len(1)));
}
closedir(dir);
std::sort(times.begin(), times.end());
return "";
}
std::pair<const BaseParameters*, MString> AVISOLOCALData::Parameters(michlib_internal::ParameterListEx& pars, const CLArgs& args) const
{
std::unique_ptr<struct Parameters> ppar{new struct Parameters};
if(!(args.contains("lonb") && args.contains("lone") && args.contains("latb") && args.contains("late"))) return {nullptr, "Region not specified (lonb, lone, latb, late)"};
ppar->lonb = args.at("lonb").ToReal();
ppar->lone = args.at("lone").ToReal();
ppar->latb = args.at("latb").ToReal();
ppar->late = args.at("late").ToReal();
pars.SetParameter("lonb", ppar->lonb);
pars.SetParameter("latb", ppar->latb);
pars.SetParameter("lone", ppar->lone);
pars.SetParameter("late", ppar->late);
return {ppar.release(), ""};
}
AVISOLOCALData::Data AVISOLOCALData::Read(const MString& vname, const BaseParameters* ip, size_t i) const
{
if(!isOk()) return Data();
auto p = dynamic_cast<const struct Parameters*>(ip);
NCFileA nc;
MString name = "";
bool isfloat = false, isint2 = false, isint4 = false;
nc.Reset(datapath + "/uv-" + times[i].ToString() + ".nc");
if(!nc) return Data();
{
auto head = nc.Header();
for(const auto& v: head.Variables())
{
auto stname = nc.A<MString>(v.Name(), "standard_name");
if(!stname) continue;
if(StName2Name(stname) == vname)
{
name = v.Name();
isint2 = v.Type().Id() == NC_SHORT;
isint4 = v.Type().Id() == NC_INT;
isfloat = v.Type().Id() == NC_FLOAT;
}
}
}
if(!name.Exist()) // Conversion read
{
// U and U2 from u and v
if(vname == "U" || vname == "U2")
{
bool square = vname == "U2";
auto u = Read("u", ip, i);
auto v = Read("v", ip, i);
if(!(u && v)) return Data();
auto out = u;
for(size_t ind = 0; ind < out.N(); ind++)
{
if(u.IsFill(ind) || v.IsFill(ind))
out.V(ind) = out.Fillval();
else
out.V(ind) = square ? (u(ind) * u(ind) + v(ind) * v(ind)) : michlib::Hypot(u(ind), v(ind));
}
return out;
}
return Data();
}
// Direct read
if(isint2) return ReadVarRaw<int2>(nc, name, p);
if(isint4) return ReadVarRaw<int4>(nc, name, p);
if(isfloat) return ReadVarRaw<float>(nc, name, p);
return Data();
}
template<class DataType> AVISOLOCALData::Data AVISOLOCALData::ReadVarRaw(const NCFileA& nc, const MString& name, const struct AVISOLOCALData::Parameters* p) const
{
real unitmul = 1.0;
DataType fill;
real offset = 0.0, scale = 1.0;
{
auto a_fill = nc.A<DataType>(name, "_FillValue");
auto a_offset_d = nc.A<double>(name, "add_offset");
auto a_scale_d = nc.A<double>(name, "scale_factor");
auto a_offset_f = nc.A<float>(name, "add_offset");
auto a_scale_f = nc.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 cn = GetCNames(nc);
auto dn = GetDNames(nc);
auto lons = nc.VR(cn.lonname);
auto lats = nc.VR(cn.latname);
if(!(lons && lats)) return Data();
auto lonb = lons(0);
auto latb = lats(0);
auto lone = lons(dn.nx - 1);
auto late = lats(dn.ny - 1);
auto lonstep = (lone - lonb) / (dn.nx - 1);
auto latstep = (late - latb) / (dn.ny - 1);
auto dom = DetGeoDomain(lonb, lone);
real lon1 = ToGeoDomain(p->lonb, dom);
real lon2 = ToGeoDomain(p->lone, dom);
real lat1 = p->latb;
real lat2 = p->late;
auto yb = static_cast<size_t>(Floor((lat1 - latb) / latstep));
auto ye = static_cast<size_t>(Ceil((lat2 - latb) / latstep));
if(ye > dn.ny - 1) ye = dn.ny - 1;
if(yb >= ye) return Data();
auto xb = static_cast<size_t>(Floor((lon1 - lonb) / lonstep));
auto xe = static_cast<size_t>(Ceil((lon2 - lonb) / lonstep));
if(xb == xe) return Data();
auto unit = nc.A<MString>(name, "units");
if(unit && (unit.Get() == "m s-1" || unit.Get() == "m/s")) unitmul = 100.0;
Data data((xb < xe) ? (xe - xb + 1) : (dn.nx + xe - xb + 1), ye - yb + 1, lons(xb), lats(yb), lonstep, latstep);
if(xb < xe)
{
auto var = nc.V<DataType>(name, {dn.lonname, xb, xe - xb + 1}, {dn.latname, yb, ye - yb + 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 = nc.V<DataType>(name, {dn.lonname, xb}, {dn.latname, yb, ye - yb + 1});
auto var2 = nc.V<DataType>(name, {dn.lonname, 0, xe + 1}, {dn.latname, yb, ye - yb + 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;
}
UVData AVISOLOCALData::ReadUV(const BaseParameters* ip, size_t i) const
{
if(!isOk()) return UVData();
auto u = Read("u", ip, i);
auto v = Read("v", ip, i);
if(!(u && v)) return UVData();
UVData out{u.Nx(), u.Ny(), u.Lon(0, 0), u.Lat(0, 0), u.XStep(), u.YStep()};
for(size_t i = 0; i < out.N(); i++)
{
if(u(i) == Data::Fillval() || v(i) == Data::Fillval())
out.U(i) = out.V(i) = UVData::Fillval();
else
{
out.U(i) = u(i);
out.V(i) = v(i);
}
}
return out;
}

125
src/layereddata.cpp

@ -8,22 +8,7 @@ MString LayeredData::Info() const
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;
if(StName2Name(ret).Exist()) vars.emplace(StName2Name(ret));
}
}
if((vars.contains("ptemp") || vars.contains("temp")) && vars.contains("sal")) vars.emplace("pdens");
if(vars.contains("ptemp") && vars.contains("sal")) vars.emplace("temp");
if(vars.contains("temp") && vars.contains("sal")) vars.emplace("ptemp");
if(vars.contains("u") && vars.contains("v")) vars.emplace("U");
if(vars.contains("u") && vars.contains("v")) vars.emplace("U2");
for(const auto& f: nc) GetVars(f.Get(), vars);
MString svars;
{
@ -43,7 +28,7 @@ MString LayeredData::Info() const
" Time step: " + Timestep() + " seconds\n" +
" Time moments: " + NTimes() + "\n" +
" Region: (" + lonb + " : " + lone + ") x (" + latb + " : " + late + ")\n" +
" Grid: " + nx + "x" + ny + " (" + lonstep + " x " + latstep + ")\n" +
" Grid: " + dname.nx + "x" + dname.ny + " (" + lonstep + " x " + latstep + ")\n" +
" Depths:" + d + "\n" +
" Supported variables: " + svars;
// clang-format on
@ -69,12 +54,6 @@ MString LayeredData::Open(const MString& dataset)
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;
@ -82,60 +61,66 @@ MString LayeredData::Open(const MString& dataset)
}
if(nc.size() == 0) return "No urls for dataset " + dataset + " specified in config";
auto head = nc[0]->Header();
lonname = latname = "";
for(const auto& dim: head.Dimensions())
dname = GetDNames(nc[0].Get());
if(!(dname.lonname.Exist() && dname.latname.Exist()))
{
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();
}
nc.clear();
return "Can't find longitude/latitude";
}
if(!(lonname.Exist() && latname.Exist()))
if(!dname.timename.Exist())
{
nc.clear();
return "Can't find longitude/latitude";
return "Can't find time";
}
auto cn = GetCNames(nc[0].Get());
// Read times
for(auto& f: nc)
{
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());
MString ret = f.ReadTimes(cn.timename);
if(ret.Exist())
{
nc.clear();
return ret;
}
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)
if(cn.depthname.Exist())
{
auto rdepths = nc[0]->VR(cn.depthname);
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);
}
else // Surface only data
{
depths.resize(1);
depths[0]=0;
depths[0] = 0;
}
auto lons = nc[0]->VR(lonname);
auto lats = nc[0]->VR(latname);
auto lons = nc[0]->VR(cn.lonname);
auto lats = nc[0]->VR(cn.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);
lone = lons(dname.nx - 1);
late = lats(dname.ny - 1);
lonstep = (lone - lonb) / (dname.nx - 1);
latstep = (late - latb) / (dname.ny - 1);
return "";
}
@ -144,7 +129,7 @@ std::pair<const BaseParameters*, MString> LayeredData::Parameters(michlib_intern
{
std::unique_ptr<struct Parameters> ppar{new struct Parameters};
if(args.contains("layer")) ppar->layer = args.at("layer").ToInteger<size_t>();
ppar->layer = args.contains("layer") ? args.at("layer").ToInteger<size_t>() : 0;
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);
@ -159,7 +144,7 @@ std::pair<const BaseParameters*, MString> LayeredData::Parameters(michlib_intern
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->ye > dname.ny - 1) ppar->ye = dname.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));
@ -277,7 +262,9 @@ LayeredData::Data LayeredData::Read(const MString& vname, const BaseParameters*
{
if(v.Dimensions().size() == 3) nodepth = true;
if(v.Type().Id() == NC_SHORT) return ReadVarRaw<int2>(nc[id], name, tid, nodepth, p);
if(v.Type().Id() == NC_INT) return ReadVarRaw<int>(nc[id], name, tid, nodepth, p);
if(v.Type().Id() == NC_FLOAT) return ReadVarRaw<float>(nc[id], name, tid, nodepth, p);
if(v.Type().Id() == NC_DOUBLE) return ReadVarRaw<double>(nc[id], name, tid, nodepth, p);
}
return Data();
@ -317,8 +304,15 @@ template<class DataType> LayeredData::Data LayeredData::ReadVarRaw(const NC& f,
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_fill)
fill = a_fill;
else
{
if constexpr(std::is_floating_point_v<DataType>)
fill = NAN;
else
fill = -1;
}
if(a_offset_d) offset = a_offset_d;
if(a_scale_d) scale = a_scale_d;
if(a_offset_f) offset = a_offset_f;
@ -328,12 +322,13 @@ template<class DataType> LayeredData::Data LayeredData::ReadVarRaw(const NC& 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);
Data data((p->xb < p->xe) ? (p->xe - p->xb + 1) : (dname.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, {dname.lonname, p->xb, p->xe - p->xb + 1}, {dname.latname, p->yb, p->ye - p->yb + 1}, {dname.timename, i, 1})
: f->V<DataType>(name, {dname.lonname, p->xb, p->xe - p->xb + 1}, {dname.latname, p->yb, p->ye - p->yb + 1}, {dname.timename, i, 1}, {dname.depthname, p->layer, 1});
if(!var) return Data();
if(var.DimLen(0) != data.Nx() || var.DimLen(1) != data.Ny()) return Data();
@ -341,28 +336,28 @@ template<class DataType> LayeredData::Data LayeredData::ReadVarRaw(const NC& f,
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);
data(ix, iy) = (v == fill || isnan(v)) ? 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});
auto var1 = nodepth ? f->V<DataType>(name, {dname.lonname, p->xb}, {dname.latname, p->yb, p->ye - p->yb + 1}, {dname.timename, i, 1})
: f->V<DataType>(name, {dname.lonname, p->xb}, {dname.latname, p->yb, p->ye - p->yb + 1}, {dname.timename, i, 1}, {dname.depthname, p->layer, 1});
auto var2 = nodepth ? f->V<DataType>(name, {dname.lonname, 0, p->xe + 1}, {dname.latname, p->yb, p->ye - p->yb + 1}, {dname.timename, i, 1})
: f->V<DataType>(name, {dname.lonname, 0, p->xe + 1}, {dname.latname, p->yb, p->ye - p->yb + 1}, {dname.timename, i, 1}, {dname.depthname, 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);
data(ix, iy) = (v == fill || isnan(v)) ? 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);
data(ix + var1.DimLen(0), iy) = (v == fill || isnan(v)) ? Data::Fillval() : ((v * scale + offset) * unitmul);
}
}
return data;

148
src/ncfuncs.cpp

@ -0,0 +1,148 @@
#define MICHLIB_NOSOURCE
#include "ncfuncs.h"
NCFuncs::CoordNames NCFuncs::GetDNames(const NCFileA& nc)
{
CoordNames out;
auto head = nc.Header();
for(const auto& dim: head.Dimensions())
{
if(dim.Name() == "lon" || dim.Name() == "longitude")
{
out.lonname = dim.Name();
out.nx = dim.Len();
}
if(dim.Name() == "lat" || dim.Name() == "latitude")
{
out.latname = dim.Name();
out.ny = dim.Len();
}
if(dim.Name() == "depth")
{
out.depthname = dim.Name();
out.nz = dim.Len();
}
if(dim.Name() == "time")
{
out.timename = dim.Name();
out.nt = dim.Len();
}
}
return out;
}
NCFuncs::CoordNames NCFuncs::GetCNames(const NCFileA& nc)
{
CoordNames out;
auto head = nc.Header();
for(const auto& v: head.Variables()) // Try to define coordinates by attribute standard_name or attribute axis
{
auto stname = nc.A<MString>(v.Name(), "standard_name");
auto axis = nc.A<MString>(v.Name(), "axis");
bool islon = false, islat = false, isdepth = false, istime = false;
if(!(stname || axis)) continue;
if(stname && stname.Get() == "longitude") islon = true;
if(stname && stname.Get() == "latitude") islat = true;
if(stname && stname.Get() == "depth") isdepth = true;
if(stname && stname.Get() == "time") istime = true;
if(!out.lonname.Exist() && axis && axis.Get() == "X") islon = true;
if(!out.latname.Exist() && axis && axis.Get() == "Y") islat = true;
if(!out.depthname.Exist() && axis && axis.Get() == "Z") isdepth = true;
if(!out.timename.Exist() && axis && axis.Get() == "T") istime = true;
if(islon) out.lonname = v.Name();
if(islat) out.latname = v.Name();
if(isdepth) out.depthname = v.Name();
if(istime) out.timename = v.Name();
if(islon) out.nx = v.Dimensions().size();
if(islat) out.ny = v.Dimensions().size();
if(isdepth) out.nz = v.Dimensions().size();
if(istime) out.nt = v.Dimensions().size();
}
return out;
}
void NCFuncs::GetVars(const NCFileA& nc, std::set<MString>& vars)
{
auto head = nc.Header();
for(const auto& v: head.Variables())
{
auto ret = nc.A<MString>(v.Name(), "standard_name");
if(!ret) continue;
if(StName2Name(ret).Exist()) vars.emplace(StName2Name(ret));
}
if((vars.contains("ptemp") || vars.contains("temp")) && vars.contains("sal")) vars.emplace("pdens");
if(vars.contains("ptemp") && vars.contains("sal")) vars.emplace("temp");
if(vars.contains("temp") && vars.contains("sal")) vars.emplace("ptemp");
if(vars.contains("u") && vars.contains("v")) vars.emplace("U");
if(vars.contains("u") && vars.contains("v")) vars.emplace("U2");
}
std::tuple<MDateTime, time_t, bool> NCFuncs::Refdate(const MString& refdate)
{
MDateTime out;
time_t step = 0;
MString rstr;
auto words = michlib::Split_on_words(refdate);
auto ci = words.begin();
if(ci != words.end())
{
if(*ci == "hours") step = 3600;
if(*ci == "days") step = 3600 * 24;
ci++;
}
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
bool success = out.FromString(rstr);
return {out, step, success};
}
MString NCFuncs::StName2Name(const MString& stname)
{
if(stname == "sea_water_potential_temperature") return "ptemp";
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 == "eastward_sea_water_velocity") return "u";
if(stname == "northward_sea_water_velocity") return "v";
if(stname == "surface_geostrophic_eastward_sea_water_velocity") return "u";
if(stname == "surface_geostrophic_northward_sea_water_velocity") return "v";
if(stname == "upward_sea_water_velocity") return "w";
return "";
}
bool NCFuncs::HaveVar(const NCFileA& nc, const MString& vname)
{
auto head = nc.Header();
for(const auto& v: head.Variables())
{
auto stname = nc.A<MString>(v.Name(), "standard_name");
if(!stname) continue;
if(StName2Name(stname) == vname) return true;
}
return false;
}
/*
template<class HV> bool NCFuncs::CheckVar(const MString& vname, HV hv)
{
if(!hv(vname))
{
bool varexist = false;
if(vname == "temp" && hv("ptemp") && hv("sal")) varexist = true;
if(vname == "ptemp" && hv("temp") && hv("sal")) varexist = true;
if(vname == "pdens" && (hv("ptemp") || hv("temp")) && hv("sal")) varexist = true;
if((vname == "U" || vname == "U2") && hv("u") && hv("v")) varexist = true;
if(!varexist) return false;
}
return true;
}
*/

11
src/odm.cpp

@ -6,7 +6,7 @@ inline void Usage(const MString& arg0)
message("Keys are:");
message(" action. What the program should do. May be: info, tsc, uv. Default: tsc.");
message(" Keys for action=info. Print some information about dataset.");
message(" source. Required. May be: NEMO, HYCOM, AVISO");
message(" source. Required. May be: NEMO, HYCOM, AVISO, AVISOLOCAL");
message(" Keys for source=NEMO");
message(" dataset. Can be DT, NRT or NRT6. Default: DT");
message(" Keys for source=HYCOM");
@ -14,7 +14,7 @@ inline void Usage(const MString& arg0)
message(" Keys for source=AVISO");
message(" dataset. Can be DT, NRT, EckmanDT or EckmanNRT. Default: DT");
message(" Keys for action=tsc. Get temperature, salinity, chlorofill from dataset.");
message(" source. Required. May be: NEMO, HYCOM, AVISO");
message(" source. Required. May be: NEMO, HYCOM, AVISO, AVISOLOCAL");
message(" var. Required. May be: U, U2, u, v, temp, ptemp, pdens, sal, chl, mld, ssh or w.");
message(" time. Time moment or regular expression. If present, timeb and timee must be absent");
message(" timeb, timee. Time interval. If present, time must be absent");
@ -30,7 +30,7 @@ inline void Usage(const MString& arg0)
message(" dataset. Can be DT, NRT, EckmanDT or EckmanNRT. Default: DT");
message(" layer and/or depth. Layer or depth of AVISO dataset. If depth is specified, layer is ignored. Both ignored for datasets DT and NRT. Default: layer=0");
message(" Keys for action=uv. Get velocity field and its derivatives.");
message(" source. Required. May be: NEMO, HYCOM, AVISO");
message(" source. Required. May be: NEMO, HYCOM, AVISO, AVISOLOCAL");
message(" time. Time moment or regular expression. If present, timeb and timee must be absent");
message(" timeb, timee. Time interval. If present, time must be absent");
message(" out. Output file for components of velocity field, divergency, rotor and Okubo-Weiss parameter. If absent, this data not calculated.");
@ -68,10 +68,7 @@ int main(int argc, char** argv)
else if(action == "tsc")
ret = data.ActionTsc(args);
else if(action == "uv")
{
args["var"] = "U";
ret = data.ActionUV(args);
}
ret = data.ActionUV(args);
else
{
errmessage("Unknown action " + action);

Loading…
Cancel
Save