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

Functions for filling Adapter from CF-compatible datasets

master
Michael Uleysky 4 weeks ago
parent
14c1244820
commit
b71dfc8340
2 changed files with 535 additions and 0 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. 66
      include/CF.h
  2. 469
      src/CF.cpp

66
include/CF.h
Unescape View File

@ -0,0 +1,66 @@
#pragma once
#include "Adapter.h"
class CF
{
public:
struct AxisInfo
{
MString tvarname, zvarname, yvarname, xvarname;
MString tdimname, zdimname, ydimname, xdimname;
uint tdims = 0, zdims = 0, ydims = 0, xdims = 0;
};
// Get time variable name
static MString GetTVarName(const NCZarr& nc);
// Get info about dimensions
static struct AxisInfo GetAxisInfo(const NCZarr& nc);
// Parse reference date and time step from time description
static std::tuple<MDateTime, time_t, bool> ParseRefDate(const MString& refdate);
// Read time moments from file
static RetVal<std::vector<MDateTime>> ReadTimes(const NCZarr& nc, const MString& tname);
// Check, if vector data lays on uniform grid
static bool CheckUniform(std::ranges::random_access_range auto&& data, real tolerance = 1e-8)
requires std::convertible_to<std::ranges::range_value_t<decltype(data)>, real>
{
if(data.size() < 2) return false;
if(data.size() == 2) return true;
real step = data[1] - data[0];
for(size_t i = 1; i < data.size() - 1; i++)
if(michlib::Abs(1.0 - (data[i + 1] - data[i]) / step) > tolerance) { return false; }
return true;
}
static Error FillAdapterFromCF(const CLArgs& args, michlib_internal::ParameterListEx& pars, Adapter& ad, std::unique_ptr<NCZarr>&& pnc);
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 == "upward_sea_water_velocity") return "w";
if(stname == "surface_geostrophic_eastward_sea_water_velocity") return "ugs";
if(stname == "surface_geostrophic_northward_sea_water_velocity") return "vgs";
if(stname == "mass_concentration_of_chlorophyll_a_in_sea_water") return "chl";
if(stname == "mole_concentration_of_nitrate_in_sea_water") return "NO3";
if(stname == "net_primary_production_of_biomass_expressed_as_carbon_per_unit_volume_in_sea_water") return "prprod";
if(stname == "mole_concentration_of_phytoplankton_expressed_as_carbon_in_sea_water") return "Cchl";
if(stname == "mole_concentration_of_phosphate_in_sea_water") return "PO4";
if(stname == "mole_concentration_of_silicate_in_sea_water") return "Si";
if(stname == "mole_concentration_of_dissolved_molecular_oxygen_in_sea_water") return "O2";
return "";
}
//static RetVal<std::shared_ptr<Data2D>> PTemp3DReader(const Adapter& ad, const void* vinfo, std::shared_ptr<Projection> proj, std::shared_ptr<Vertical> vert, size_t it);
};

469
src/CF.cpp
Unescape View File

@ -0,0 +1,469 @@
#define MICHLIB_NOSOURCE
#include "CF.h"
using michlib::DetGeoDomain;
MString CF::GetTVarName(const NCZarr& nc)
{
for(const auto& v: nc.VarNames()) // Try to define coordinates by attribute standard_name or attribute axis
{
auto havestname = nc.HasAtt(v, "standard_name");
auto haveaxis = nc.HasAtt(v, "axis");
if(!(havestname || haveaxis)) continue;
auto stname = nc.AttString(v, "standard_name");
auto axis = nc.AttString(v, "axis");
if(stname == "time" || axis == "T") return v;
}
// If time not found just check variable "time"
if(nc.HasVar("time")) return "time";
return "";
}
struct CF::AxisInfo CF::GetAxisInfo(const NCZarr& nc)
{
struct AxisInfo out;
for(const auto& v: nc.VarNames()) // Try to define coordinates by attribute standard_name or attribute axis
{
auto havestname = nc.HasAtt(v, "standard_name");
auto haveaxis = nc.HasAtt(v, "axis");
if(!(havestname || haveaxis)) continue;
auto stname = nc.AttString(v, "standard_name");
auto axis = nc.AttString(v, "axis");
bool islon = false, islat = false, isdepth = false, istime = false;
if(stname == "longitude") islon = true;
if(stname == "latitude") islat = true;
if(stname == "depth") isdepth = true;
if(stname == "time") istime = true;
if(!out.xvarname.Exist() && axis == "X") islon = true;
if(!out.yvarname.Exist() && axis == "Y") islat = true;
if(!out.zvarname.Exist() && axis == "Z") isdepth = true;
if(!out.tvarname.Exist() && axis == "T") istime = true;
if(islon) out.xvarname = v;
if(islat) out.yvarname = v;
if(isdepth) out.zvarname = v;
if(istime) out.tvarname = v;
}
// If time not found just check variable "time"
if(!out.tvarname.Exist() && nc.HasVar("time")) out.tvarname = "time";
// Dimensions
out.tdims = out.tvarname.Exist() ? nc.NDim(out.tvarname) : 0;
out.zdims = out.zvarname.Exist() ? nc.NDim(out.zvarname) : 0;
out.ydims = out.yvarname.Exist() ? nc.NDim(out.yvarname) : 0;
out.xdims = out.xvarname.Exist() ? nc.NDim(out.xvarname) : 0;
if(out.tdims == 1) out.tdimname = nc.DimNames(out.tvarname)[0];
if(out.zdims == 1) out.zdimname = nc.DimNames(out.zvarname)[0];
if(out.ydims == 1) out.ydimname = nc.DimNames(out.yvarname)[0];
if(out.xdims == 1) out.xdimname = nc.DimNames(out.xvarname)[0];
return out;
}
std::tuple<MDateTime, time_t, bool> CF::ParseRefDate(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 == "seconds") step = 1;
if(*ci == "minutes") step = 60;
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()) rstr += " " + *ci; // Hours
bool success = out.FromString(rstr);
return {out, step, success};
}
RetVal<std::vector<MDateTime>> CF::ReadTimes(const NCZarr& nc, const MString& tname)
{
static const MString pref = "CF::ReadTimes";
std::vector<double> time;
std::vector<MDateTime> out;
{
auto ret = nc.Read(tname, time);
if(!ret) return ret.Add(pref, "Can't read time");
}
MDateTime refdate;
time_t step = 0;
{
auto units = nc.AttString(tname, "units");
if(!units.Exist()) return {pref, "Can't read refdate"};
auto [rd, st, suc] = ParseRefDate(units);
if(!suc) return {pref, "Can't parse " + units + " to refdate"};
if(st == 0) return {pref, "Can't get timestep from string " + units};
refdate = rd;
step = st;
}
out.resize(time.size());
for(size_t i = 0; i < time.size(); i++) out[i] = refdate + static_cast<time_t>(time[i] * step);
return out;
}
Error CF::FillAdapterFromCF(const CLArgs& args, michlib_internal::ParameterListEx& pars, Adapter& ad, std::unique_ptr<NCZarr>&& pnc)
{
static const MString pref = "CF::FillAdapterFromCF";
bool debug = args.contains("debug");
auto axisinfo = GetAxisInfo(*pnc);
if(!axisinfo.tvarname.Exist()) return {pref, "Can't find time variable in the dataset"};
if(!axisinfo.yvarname.Exist()) return {pref, "Can't find latitude variable in the dataset"};
if(!axisinfo.xvarname.Exist()) return {pref, "Can't find longitude variable in the dataset"};
if(axisinfo.tdims != 1) return {pref, "Unsupported number of time dimensions: " + MString(axisinfo.tdims)};
if(axisinfo.ydims != 1) return {pref, "Unsupported number of latitude dimensions: " + MString(axisinfo.ydims)};
if(axisinfo.xdims != 1) return {pref, "Unsupported number of longitude dimensions: " + MString(axisinfo.xdims)};
if(axisinfo.zvarname.Exist() && axisinfo.zdims != 1) return {pref, "Unsupported number of depth dimensions: " + MString(axisinfo.zdims)};
if(debug)
{
pars.AddParameter("T variable", axisinfo.tvarname + "/" + axisinfo.tdims + "/" + axisinfo.tdimname);
pars.AddParameter("Z variable", axisinfo.zvarname + "/" + axisinfo.zdims + "/" + axisinfo.zdimname);
pars.AddParameter("Y variable", axisinfo.yvarname + "/" + axisinfo.ydims + "/" + axisinfo.ydimname);
pars.AddParameter("X variable", axisinfo.xvarname + "/" + axisinfo.xdims + "/" + axisinfo.xdimname);
}
std::shared_ptr<Projection> proj;
std::shared_ptr<Vertical> vert;
// Forming ReadInfo
std::shared_ptr<Adapter::ReadInfo> rinfo(new Adapter::ReadInfo);
rinfo->xdname = axisinfo.xdimname;
rinfo->ydname = axisinfo.ydimname;
rinfo->zdname = axisinfo.zdimname;
rinfo->tdname = axisinfo.tdimname;
std::vector<double> depths;
if(axisinfo.zvarname.Exist())
{
auto ret = pnc->Read(axisinfo.zvarname, depths);
if(!ret) return ret.Add(pref, "Can't read depths");
}
// Filling xb, yb, xe, ye
// We assume that latitude and longitude are sorted in ascending order
{
std::vector<double> lons, lats;
{
auto ret = pnc->Read(axisinfo.xvarname, lons);
if(!ret) return ret.Add(pref, "Can't read longitudes");
}
{
auto ret = pnc->Read(axisinfo.yvarname, lats);
if(!ret) return ret.Add(pref, "Can't read latitudes");
}
real lonb, latb, lone, late;
bool global = lons.back() - lons.front() + 1.5 * (lons.back() - lons.front()) / (lons.size() - 1) > 360.0;
auto dom = DetGeoDomain(lons.front(), lons.back());
lonb = ToGeoDomain(args.contains("lonb") ? args.at("lonb").ToReal() : lons.front(), dom);
lone = ToGeoDomain(args.contains("lone") ? args.at("lone").ToReal() : lons.back(), dom);
latb = args.contains("latb") ? args.at("latb").ToReal() : lats.front();
late = args.contains("late") ? args.at("late").ToReal() : lats.back();
if(global && lonb > lons.back()) return {pref, "Argument lonb must be lesser then " + MString(lons.back())};
if(global && lone < lons.front()) return {pref, "Argument lone must be greater then " + MString(lons.front())};
if(latb > lats.back()) return {pref, "Argument latb must be lesser then " + MString(lats.back())};
if(late < lats.front()) return {pref, "Argument late must be greater then " + MString(lats.front())};
if(global)
{
// Special case when the longitude lies in a small sector between the end and the start
if(lonb < lons.front()) lonb = lons.back();
if(lone > lons.back()) lone = lons.front();
}
else
{
if(lonb < lons.front()) lonb = lons.front();
if(lone > lons.back()) lone = lons.back();
}
if(latb < lats.front()) latb = lats.front();
if(late > lats.back()) late = lats.back();
for(rinfo->xb = 0; rinfo->xb + 1 < lons.size() && lonb >= lons[rinfo->xb + 1]; rinfo->xb++);
for(rinfo->yb = 0; rinfo->yb + 1 < lats.size() && latb >= lats[rinfo->yb + 1]; rinfo->yb++);
for(rinfo->xe = lons.size() - 1; rinfo->xe >= 1 && lone <= lons[rinfo->xe - 1]; rinfo->xe--);
for(rinfo->ye = lats.size() - 1; rinfo->ye >= 1 && late <= lats[rinfo->ye - 1]; rinfo->ye--);
if(debug)
{
pars.AddParameter("Requested lonb", lonb);
pars.AddParameter("Requested lone", lone);
pars.AddParameter("Requested latb", latb);
pars.AddParameter("Requested late", late);
pars.AddParameter("xb", rinfo->xb);
pars.AddParameter("xe", rinfo->xe);
pars.AddParameter("yb", rinfo->yb);
pars.AddParameter("ye", rinfo->ye);
}
if(rinfo->xb == rinfo->xe) return {pref, "Longitude interval too small"};
if(rinfo->yb == rinfo->ye) return {pref, "Latitude interval too small"};
if(rinfo->xe == 0) return {pref, "Sorry, special case then xe=0 does'nt handle. Increase lone."};
if(rinfo->xb == lons.size() - 1) return {pref, "Sorry, special case then xb=lons.size()-1 does'nt handle. Decrease lonb."};
lonb = lons[rinfo->xb];
lone = lons[rinfo->xe];
latb = lats[rinfo->yb];
late = lats[rinfo->ye];
if(lonb > lone)
{
if(lonb > 0.0)
lonb -= 360.0;
else
lone += 360.0;
}
pars.AddParameter("lonb", michlib::ToGeoDomainPos(lonb));
pars.AddParameter("lone", michlib::ToGeoDomainPos(lone));
pars.AddParameter("latb", latb);
pars.AddParameter("late", late);
decltype(lats) rlats = lats | std::views::drop(rinfo->yb) | std::views::take(rinfo->ye - rinfo->yb + 1) | std::ranges::to<std::vector>();
decltype(lons) rlons;
if(rinfo->xe >= rinfo->xb)
rlons = lons | std::views::drop(rinfo->xb) | std::views::take(rinfo->xe - rinfo->xb + 1) | std::ranges::to<std::vector>();
else
{
rlons.insert(rlons.begin(), lons.begin() + rinfo->xb, lons.end());
std::ranges::transform(lons.begin(), lons.begin() + rinfo->xe + 1, std::back_inserter(rlons), [](auto lon) { return lon + 360.0; });
}
// Check uniform and create projection
if(CF::CheckUniform(rlons, 0.0005) && CF::CheckUniform(rlats, 0.0005))
{
proj.reset(new Projection(Projection::Create<Projection::Type::EQC>(rlons, rlats)));
pars.AddParameter("lonstep", (rlons.back() - rlons.front()) / (rlons.size() - 1));
pars.AddParameter("latstep", (rlats.back() - rlats.front()) / (rlats.size() - 1));
}
else
{
michlib::errmessage("Warning! Nonuniform grid");
proj.reset(new Projection(Projection::Create<Projection::Type::IRREGULARXY>(std::move(rlons), std::move(rlats))));
}
}
// ReadInfo for the plain data
std::shared_ptr<Adapter::ReadInfo> rinfoplane(new Adapter::ReadInfo(*rinfo));
rinfoplane->zdname = "";
rinfoplane->zb = rinfoplane->ze = rinfo->zb = rinfo->ze = 0;
// Filling zb, ze
bool depthinv = false;
if(depths.size() > 1)
{
if(depths.back() <= 0 && depths.front() <= 0) std::ranges::transform(depths, depths.begin(), std::negate{});
if(depths.back() < depths.front() && depths.size() > 1)
{
depthinv = true;
for(size_t i = 0; i < depths.size() - i - 1; i++) std::swap(depths[i], depths[depths.size() - i - 1]);
rinfoplane->zb = rinfoplane->ze = rinfo->zb = rinfo->ze = depths.size() - 1;
}
if(args.contains("layer") || args.contains("depth"))
{
size_t layer = 0;
if(args.contains("layer"))
layer = args.at("layer").ToInteger<size_t>();
else
{
auto d = args.at("depth").ToReal();
if(d < depths.front()) d = depths.front();
if(d > depths.back()) d = depths.back();
for(size_t i = 0; i < depths.size() - 1; i++)
if(d >= depths[i] && d <= depths[i + 1])
{
layer = (d - depths[i] <= depths[i + 1] - d) ? i : (i + 1);
break;
}
if(debug) pars.AddParameter("Requested depth", d);
}
if(layer > depths.size() - 1) return {pref, "Layer must be lesser then " + MString(depths.size())};
rinfoplane->zb = rinfoplane->ze = rinfo->zb = rinfo->ze = depthinv ? depths.size() - layer - 1 : layer;
vert.reset(new Vertical(Vertical::Create<Vertical::Type::FIXED>(depths[layer])));
pars.AddParameter("layer", layer);
pars.AddParameter("depth", depths[layer]);
if(debug) pars.AddParameter("zb", rinfo->zb);
}
else if(args.contains("layers") || args.contains("depths"))
{
size_t layer1 = 0, layer2 = 0;
if(args.contains("layers"))
{
if(args.at("layers") == "ALL" || args.at("layers") == "all")
{
layer1 = 0;
layer2 = depths.size() - 1;
}
else
{
auto layerstr = args.at("layers").Split(":");
if(layerstr.size() > 2 || layerstr.size() == 0) return {pref, "Layers must have format layer1:layer2"};
if(layerstr.size() == 1)
layer1 = layer2 = layerstr[0].ToInteger<size_t>();
else
{
layer1 = layerstr[0].ToInteger<size_t>();
layer2 = layerstr[1].ToInteger<size_t>();
}
}
}
else
{
real d1 = 0.0, d2 = 0.0;
auto depthstr = args.at("depths").Split(":");
if(depthstr.size() > 2 || depthstr.size() == 0) return {pref, "Depths must have format depth1:depth2"};
if(depthstr.size() == 1)
d1 = d2 = depthstr[0].ToReal();
else
{
d1 = depthstr[0].ToReal();
d2 = depthstr[1].ToReal();
}
if(d1 < depths.front()) d1 = depths.front();
if(d1 > depths.back()) d1 = depths.back();
if(d2 < depths.front()) d2 = depths.front();
if(d2 > depths.back()) d2 = depths.back();
for(size_t i = 0; i < depths.size() - 1; i++)
if(d1 >= depths[i] && d1 <= depths[i + 1])
{
layer1 = (d1 - depths[i] <= depths[i + 1] - d1) ? i : (i + 1);
break;
}
for(size_t i = 0; i < depths.size() - 1; i++)
if(d2 >= depths[i] && d2 <= depths[i + 1])
{
layer2 = (d2 - depths[i] <= depths[i + 1] - d2) ? i : (i + 1);
break;
}
if(debug)
{
pars.AddParameter("Requested depth1", d1);
pars.AddParameter("Requested depth2", d2);
}
}
if(layer1 > layer2) return {pref, "Layer2 must be greater then layer1"};
if(layer1 > depths.size() - 1) return {pref, "Layer1 must be lesser then " + MString(depths.size())};
if(layer2 > depths.size() - 1) return {pref, "Layer2 must be lesser then " + MString(depths.size())};
rinfoplane->zb = rinfoplane->ze = depthinv ? depths.size() - layer1 - 1 : layer1;
rinfo->zb = depthinv ? depths.size() - layer1 - 1 : layer1;
rinfo->ze = depthinv ? depths.size() - layer2 - 1 : layer2;
pars.AddParameter("layer1", layer1);
pars.AddParameter("depth1", depths[layer1]);
pars.AddParameter("layer2", layer2);
pars.AddParameter("depth2", depths[layer2]);
if(debug)
{
pars.AddParameter("zb", rinfo->zb);
pars.AddParameter("ze", rinfo->ze);
}
// Create vertical projection
{
decltype(depths) rdepths(depths.begin() + layer1, depths.begin() + layer2 + 1);
vert.reset(new Vertical(Vertical::Create<Vertical::Type::HORIZONTS>(rdepths)));
}
}
else
{
pars.AddParameter("layer", 0);
pars.AddParameter("depth", depths[0]);
if(debug) pars.AddParameter("zb", rinfo->zb);
vert.reset(new Vertical(Vertical::Create<Vertical::Type::FIXED>(depths[0])));
}
}
if(debug) pars.AddParameter("depthinv", depthinv);
//michlib::message("depth1=", depths[depthinv ? depths.size() - 1 - rinfo->zb : rinfo->zb], ", zb=", rinfo->zb);
//michlib::message("depth2=", depths[depthinv ? depths.size() - 1 - rinfo->ze : rinfo->ze], ", ze=", rinfo->ze);
//michlib::message(depthinv ? "Depths inverted" : "Depths not inverted");
// Read variables information
for(const auto& v: pnc->VarNames())
{
Adapter::VInfo vinfo;
auto dnames = pnc->DimNames(v);
if(dnames.size() != 3 && dnames.size() != 4) continue;
if(dnames.size() == 3 && (dnames[0] != axisinfo.tdimname || dnames[1] != axisinfo.ydimname || dnames[2] != axisinfo.xdimname)) continue;
if(dnames.size() == 4 && (dnames[0] != axisinfo.tdimname || dnames[1] != axisinfo.zdimname || dnames[2] != axisinfo.ydimname || dnames[3] != axisinfo.xdimname)) continue;
vinfo.name = v;
vinfo.unit = "1";
if(dnames.size() == 4)
vinfo.rinfo = rinfo;
else
vinfo.rinfo = rinfoplane;
if(pnc->HasAtt(v, "units")) vinfo.unit = pnc->AttString(v, "units");
if(pnc->HasAtt(v, "standard_name")) vinfo.stname = pnc->AttString(v, "standard_name");
if(pnc->HasAtt(v, "long_name")) vinfo.lname = pnc->AttString(v, "long_name");
vinfo.offset = pnc->HasAtt(v, "add_offset") ? pnc->AttReal(v, "add_offset") : 0.0;
vinfo.scale = pnc->HasAtt(v, "scale_factor") ? pnc->AttReal(v, "scale_factor") : 1.0;
vinfo.fillval = std::visit(
[](auto r) -> real
{
if constexpr(std::is_convertible_v<decltype(r), real>)
return static_cast<real>(r);
else
return NAN;
},
pnc->VarFill(v));
vinfo.targetunit = vinfo.unit;
if(vinfo.stname.SubStr(vinfo.stname.Len() - 7, 8) == "velocity") vinfo.targetunit = "cm/s";
MString name = StName2Name(vinfo.stname).Exist() ? StName2Name(vinfo.stname) : v;
if(dnames.size() == 3) { ad.Add2DVariable(name, std::move(vinfo), proj); }
else { ad.Add3DVariable(name, std::move(vinfo), proj, vert); }
}
return Error();
}
/*
RetVal<std::shared_ptr<Data2D>> CF::PTemp3DReader(const Adapter& ad, const void* vinfo, std::shared_ptr<Projection> proj, std::shared_ptr<Vertical> vert, size_t it)
{
static const MString pref = "CF::PTemp3DReader";
auto v = michlib::pointer_cast<const struct Adapter::VInfo*>(vinfo);
const struct Adapter::ReadInfo* ri = v->rinfo.get();
std::shared_ptr<Data3D> out(new Data3D(proj, vert, {.unit = v->targetunit, .stname = v->stname, .lname = v->lname, .comment = v->comment}));
}
*/
Write Preview
Loading…
Cancel
Save