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.
125 Commits
3 Branches
0 Tags
1.7 MiB
 
 
 
Tree: 85212ae018
Branches Tags
${ item.name }
Create tag ${ searchTerm }
Create branch ${ searchTerm }
from '85212ae018'
${ noResults }
odm/include/ncfuncs.h

232 lines
8.9 KiB
Raw Blame History

#pragma once
#include "DataAdapters/ncfilealt.h"
#include "basedata.h"
#include "mdatetime.h"
#include "nczarr.h"
#include <map>
#include <set>
#include <tuple>
using michlib::MDateTime;
using michlib::MString;
using michlib::NCFileA;
class NCFuncs
{
public:
struct CoordNames
{
MString lonname, latname, depthname, timename;
};
struct DimNames
{
MString lonname, latname, depthname, timename;
size_t nx, ny, nz, nt;
};
static MString StName2Name(const MString& stname);
static MString Name2StName(const MString& name);
static MString Name2LongName(const MString& name);
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 DimNames GetDNames(const NCFileA& nc);
static bool HaveVar(const NCFileA& nc, const MString& vname);
static void GetVars(const NCZarr& nc, std::set<MString>& vars);
static CoordNames GetCNames(const NCZarr& nc);
static DimNames GetDNames(const NCZarr& nc);
static bool HaveVar(const NCZarr& nc, const MString& vname);
template<class HV> static VarPresence CheckVar(const MString& vname, HV hv)
{
if(!hv(vname))
{
auto varexist = VarPresence::NONE;
if(vname == "temp" && hv("ptemp") && hv("sal")) varexist = VarPresence::DERIVED;
if(vname == "ptemp" && hv("temp") && hv("sal")) varexist = VarPresence::DERIVED;
if(vname == "pdens" && (hv("ptemp") || hv("temp")) && hv("sal")) varexist = VarPresence::DERIVED;
if((vname == "U" || vname == "U2") && hv("u") && hv("v")) varexist = VarPresence::DERIVED;
if((vname == "ugeo" || vname == "vgeo") && hv("ssh")) varexist = VarPresence::DERIVED;
return varexist;
}
return VarPresence::INTERNAL;
}
template<class D> static bool TransformationRead(const D* data, const MString& vname, std::map<MString, typename D::Data>& cache, const BaseParameters* ip, size_t i)
{
if constexpr(requires(const D* d, const BaseParameters* p) {
{ d->Depth(p) } -> std::convertible_to<real>;
})
{
real depth = data->Depth(ip);
// ptemp from temp and sal
if(vname == "ptemp")
{
if(!(data->Read("temp", cache, ip, i) && data->Read("sal", cache, ip, i))) return false;
cache[vname] = cache.at("temp").CopyGrid();
auto& ptemp = cache.at(vname);
const auto& temp = cache.at("temp");
const auto& sal = cache.at("sal");
ptemp.SetUnit("degrees_C");
for(size_t ind = 0; ind < ptemp.N(); ind++)
{
if(temp.IsFill(ind) || sal.IsFill(ind))
ptemp.V(ind) = ptemp.Fillval();
else
ptemp.V(ind) = Temp2PTemp(temp.V(ind), sal.V(ind), depth, ptemp.Lon(ind), ptemp.Lat(ind));
}
return true;
}
// temp from ptemp and sal
if(vname == "temp")
{
if(!(data->Read("ptemp", cache, ip, i) && data->Read("sal", cache, ip, i))) return false;
cache[vname] = cache.at("ptemp").CopyGrid();
auto& temp = cache.at(vname);
const auto& ptemp = cache.at("ptemp");
const auto& sal = cache.at("sal");
temp.SetUnit("degrees_C");
for(size_t ind = 0; ind < temp.N(); ind++)
{
if(ptemp.IsFill(ind) || sal.IsFill(ind))
temp.V(ind) = temp.Fillval();
else
temp.V(ind) = PTemp2Temp(ptemp.V(ind), sal.V(ind), depth, temp.Lon(ind), temp.Lat(ind));
}
return true;
}
// pdens from temp and sal
if(vname == "pdens")
{
bool tempispot = data->CheckVar("ptemp") == VarPresence::INTERNAL;
if(!(data->Read(tempispot ? "ptemp" : "temp", cache, ip, i) && data->Read("sal", cache, ip, i))) return false;
cache[vname] = cache.at("sal").CopyGrid();
auto& pdens = cache.at(vname);
const auto& temp = cache.at(tempispot ? "ptemp" : "temp");
const auto& sal = cache.at("sal");
pdens.SetUnit("kg m-3");
for(size_t ind = 0; ind < pdens.N(); ind++)
{
if(temp.IsFill(ind) || sal.IsFill(ind))
pdens.V(ind) = pdens.Fillval();
else
pdens.V(ind) =
tempispot ? PTemp2PDens(temp.V(ind), sal.V(ind), depth, pdens.Lon(ind), pdens.Lat(ind)) : Temp2PDens(temp.V(ind), sal.V(ind), depth, pdens.Lon(ind), pdens.Lat(ind));
}
return true;
}
}
// U and U2 from u and v
if(vname == "U" || vname == "U2")
{
bool square = vname == "U2";
if(!(data->Read("u", cache, ip, i) && data->Read("v", cache, ip, i))) return false;
cache[vname] = cache.at("u").CopyGrid();
auto& U = cache.at(vname);
const auto& u = cache.at("u");
const auto& v = cache.at("v");
if(u.Unit().Exist()) U.SetUnit(square ? ("(" + u.Unit() + ")2") : u.Unit());
for(size_t ind = 0; ind < U.N(); ind++)
{
if(u.IsFill(ind) || v.IsFill(ind))
U.V(ind) = U.Fillval();
else
U.V(ind) = square ? (u(ind) * u(ind) + v(ind) * v(ind)) : michlib::Hypot(u(ind), v(ind));
}
return true;
}
// ugeo
if(vname == "ugeo")
{
if(!data->Read("ssh", cache, ip, i)) return false;
const auto& ssh = cache.at("ssh");
const real mul = (ssh.Unit() == "cm") ? 1.0 : ((ssh.Unit() == "m") ? 100.0 : -1.0);
if(mul < 0.0) return false; // Unknown units
cache[vname] = cache.at("ssh").CopyGrid();
auto& ugeo = cache.at(vname);
ugeo.SetUnit("cm/s");
for(size_t iy = 0; iy < ssh.Ny(); iy++)
{
const real phi = michlib::M_PI * ssh.Iy2Lat(iy) / 180.0;
const real sphi = michlib::Sin(phi);
const real s2phi = michlib::Sin(2.0 * phi);
const real g = 978.0318 * (1.0 + 0.005302 * sphi * sphi - 0.000006 * s2phi * s2phi); // g in cm/s
const real f = 2.0 * 7.29211 * sphi; // Coriolis parameter multiplied by 100'000
const real db = (iy == 0) ? 0.0 : (michlib::M_PI * 6371.0 * (ssh.Iy2Lat(iy) - ssh.Iy2Lat(iy - 1)) / 180.0); // Distance in km, compensated by multiplier in Coriolis parameter
const real du =
(iy == ssh.Ny() - 1) ? 0.0 : (michlib::M_PI * 6371.0 * (ssh.Iy2Lat(iy + 1) - ssh.Iy2Lat(iy)) / 180.0); // Distance in km, compensated by multiplier in Coriolis parameter
for(size_t ix = 0; ix < ssh.Nx(); ix++)
{
const bool bok = iy != 0 && !ssh.IsFill(ix, iy - 1) && !ssh.IsFill(ix, iy);
const bool uok = iy != ssh.Ny() - 1 && !ssh.IsFill(ix, iy + 1) && !ssh.IsFill(ix, iy);
if(!(bok || uok))
ugeo(ix, iy) = ugeo.Fillval();
else if(!bok)
ugeo(ix, iy) = -g * mul * (ssh(ix, iy + 1) - ssh(ix, iy)) / (f * du);
else if(!uok)
ugeo(ix, iy) = -g * mul * (ssh(ix, iy) - ssh(ix, iy - 1)) / (f * db);
else
{
const real bder = -g * mul * (ssh(ix, iy) - ssh(ix, iy - 1)) / (f * db);
const real uder = -g * mul * (ssh(ix, iy + 1) - ssh(ix, iy)) / (f * du);
ugeo(ix, iy) = 0.5 * (bder + uder);
}
}
}
return true;
}
// vgeo
if(vname == "vgeo")
{
if(!data->Read("ssh", cache, ip, i)) return false;
const auto& ssh = cache.at("ssh");
const real mul = (ssh.Unit() == "cm") ? 1.0 : ((ssh.Unit() == "m") ? 100.0 : -1.0);
if(mul < 0.0) return false; // Unknown units
cache[vname] = cache.at("ssh").CopyGrid();
auto& vgeo = cache.at(vname);
vgeo.SetUnit("cm/s");
for(size_t iy = 0; iy < ssh.Ny(); iy++)
{
const real phi = michlib::M_PI * ssh.Iy2Lat(iy) / 180.0;
const real sphi = michlib::Sin(phi);
const real cphi = michlib::Cos(phi);
const real s2phi = michlib::Sin(2.0 * phi);
const real g = 978.0318 * (1.0 + 0.005302 * sphi * sphi - 0.000006 * s2phi * s2phi); // g in cm/s
const real f = 2.0 * 7.29211 * sphi; // Coriolis parameter multiplied by 100'000
for(size_t ix = 0; ix < ssh.Nx(); ix++)
{
const real dl =
(ix == 0) ? 0.0 : (michlib::M_PI * 6371.0 * (ssh.Ix2Lon(ix) - ssh.Ix2Lon(ix - 1)) / 180.0) * cphi; // Distance in km, compensated by multiplier in Coriolis parameter
const real dr = (ix == ssh.Nx() - 1)
? 0.0
: (michlib::M_PI * 6371.0 * (ssh.Ix2Lon(ix + 1) - ssh.Ix2Lon(ix)) / 180.0) * cphi; // Distance in km, compensated by multiplier in Coriolis parameter
const bool lok = ix != 0 && !ssh.IsFill(ix - 1, iy) && !ssh.IsFill(ix, iy);
const bool rok = ix != ssh.Nx() - 1 && !ssh.IsFill(ix + 1, iy) && !ssh.IsFill(ix, iy);
if(!(lok || rok))
vgeo(ix, iy) = vgeo.Fillval();
else if(!lok)
vgeo(ix, iy) = g * mul * (ssh(ix + 1, iy) - ssh(ix, iy)) / (f * dr);
else if(!rok)
vgeo(ix, iy) = g * mul * (ssh(ix, iy) - ssh(ix - 1, iy)) / (f * dl);
else
{
const real lder = g * mul * (ssh(ix, iy) - ssh(ix - 1, iy)) / (f * dl);
const real rder = g * mul * (ssh(ix + 1, iy) - ssh(ix, iy)) / (f * dr);
vgeo(ix, iy) = 0.5 * (lder + rder);
}
}
}
return true;
}
return false;
}
};