Compare commits

...

4 Commits

  1. 219
      actions/actiongrad.cpp
  2. 146
      actions/actiongrad.h
  3. 7
      actions/actionuv.cpp
  4. 2
      include/actiondep.h
  5. 18
      include/ncfilew.h
  6. 129
      sources/TSCDATA.cpp
  7. 35
      sources/TSCDATA.h
  8. 12
      src/ncfilew.cpp

219
actions/actiongrad.cpp

@ -0,0 +1,219 @@
#define MICHLIB_NOSOURCE
#include "actiongrad.h"
MString GradMethods::NCFileW::Create(const MString& name, const MString& history, const std::vector<MString>& vnames, const std::vector<MString>& lnames,
const std::vector<GradMethods::DataType>& lons, const std::vector<GradMethods::DataType>& lats, int compress)
{
int ret;
const float node_offset = 0.0;
auto nx = lons.size();
auto ny = lats.size();
int xid, yid;
MString text;
struct
{
union
{
struct
{
int ydimid, xdimid;
};
int dimid[2];
};
} dim;
struct Keeper
{
bool active;
int ncid;
~Keeper()
{
if(active) nc_close(ncid);
}
operator int() { return ncid; }
} newnc;
const auto fill = static_cast<GradMethods::DataType>(std::numeric_limits<GradMethods::Matrix::MDataType>::max());
// Creating
ret = nc_create(name.Buf(), NC_CLOBBER | NC_NETCDF4, &newnc.ncid);
if(ret != NC_NOERR) return "Can't create netcdf file: " + name;
newnc.active = true;
// Global attributes
ret = nc_put_att_text(newnc, NC_GLOBAL, "history", history.Len() + 1, history.Buf());
if(ret != NC_NOERR) return "Can't write history attribute in the netcdf file";
ret = nc_put_att(newnc, NC_GLOBAL, "node_offset", NC_FLOAT, 1, &node_offset);
if(ret != NC_NOERR) return "Can't write history attribute in the netcdf file";
// Dimensions
ret = nc_def_dim(newnc, "longitude", nx, &dim.xdimid);
if(ret != NC_NOERR) return "Can't create x-dimension in the netcdf file";
ret = nc_def_dim(newnc, "latitude", ny, &dim.ydimid);
if(ret != NC_NOERR) return "Can't create y-dimension in the netcdf file";
// lon, lat variables
ret = nc_def_var(newnc, "longitude", NC_FLOAT, 1, &dim.xdimid, &xid);
if(ret != NC_NOERR) return "Can't create longitude variable in the netcdf file";
ret = nc_def_var(newnc, "latitude", NC_FLOAT, 1, &dim.ydimid, &yid);
if(ret != NC_NOERR) return "Can't create latitude variable in the netcdf file";
ret = nc_def_var_deflate(newnc, xid, 1, 1, compress);
if(ret != NC_NOERR) return "Can't set deflate parameters for longitude variable in the netcdf file";
ret = nc_def_var_deflate(newnc, yid, 1, 1, compress);
if(ret != NC_NOERR) return "Can't set deflate parameters for latitude variable in the netcdf file";
text = "longitude";
ret = nc_put_att_text(newnc, xid, "standard_name", text.Len() + 1, text.Buf());
if(ret != NC_NOERR) return "Can't write standard_name attribute of longitude variable in the netcdf file";
text = "latitude";
ret = nc_put_att_text(newnc, yid, "standard_name", text.Len() + 1, text.Buf());
if(ret != NC_NOERR) return "Can't write standard_name attribute of latitude variable in the netcdf file";
text = "Longitude";
ret = nc_put_att_text(newnc, xid, "long_name", text.Len() + 1, text.Buf());
if(ret != NC_NOERR) return "Can't write long_name attribute of longitude variable in the netcdf file";
text = "Latitude";
ret = nc_put_att_text(newnc, yid, "long_name", text.Len() + 1, text.Buf());
if(ret != NC_NOERR) return "Can't write long_name attribute of latitude variable in the netcdf file";
// Variables
for(size_t i = 0; i < vnames.size(); i++)
{
int vid;
ret = nc_def_var(newnc, vnames[i].Buf(), NC_FLOAT, 2, dim.dimid, &vid);
if(ret != NC_NOERR) return "Can't create " + vnames[i] + " variable in the netcdf file";
ret = nc_def_var_deflate(newnc, vid, 1, 1, compress);
if(ret != NC_NOERR) return "Can't set deflate parameters for " + vnames[i] + " variable in the netcdf file";
if(lnames[i].Exist()) ret = nc_put_att_text(newnc, vid, "long_name", lnames[i].Len() + 1, lnames[i].Buf());
if(ret != NC_NOERR) return "Can't write long_name attribute of " + vnames[i] + " variable in the netcdf file";
ret = nc_put_att_float(newnc, vid, "_FillValue", NC_FLOAT, 1, &fill);
if(ret != NC_NOERR) return "Can't write _FillValue attribute of " + vnames[i] + " variable in the netcdf file";
}
// End definitions
nc_enddef(newnc);
// Writing lon, lat
const size_t i = 0;
ret = nc_put_vara_float(newnc, xid, &i, &nx, lons.data());
if(ret != NC_NOERR) return "Can't write longitude variable in the netcdf file";
ret = nc_put_vara_float(newnc, yid, &i, &ny, lats.data());
if(ret != NC_NOERR) return "Can't write latitude variable in the netcdf file";
ncid = newnc;
newnc.active = false;
return "";
}
MString GradMethods::NCFileW::WriteVar(const MString& name, const std::vector<GradMethods::Matrix::MDataType>& data)
{
int ret;
int vid;
int xid, yid;
size_t nx, ny;
ret = nc_inq_varid(ncid, name.Buf(), &vid);
if(ret != NC_NOERR) return "Can't find variable " + name + " in the netcdf file";
ret = nc_inq_dimid(ncid, "longitude", &xid);
if(ret != NC_NOERR) return "Can't find longitude dimension in the netcdf file";
ret = nc_inq_dimid(ncid, "latitude", &yid);
if(ret != NC_NOERR) return "Can't find latitude dimension in the netcdf file";
ret = nc_inq_dimlen(ncid, xid, &nx);
if(ret != NC_NOERR) return "Can't find longitude size in the netcdf file";
ret = nc_inq_dimlen(ncid, yid, &ny);
if(ret != NC_NOERR) return "Can't find latitude size in the netcdf file";
const size_t i[2] = {0, 0};
const size_t c[2] = {ny, nx};
GradMethods::DataType buf[nx * ny];
for(size_t i = 0; i < nx * ny; i++) buf[i] = static_cast<GradMethods::DataType>(data[i]);
ret = nc_put_vara_float(ncid, vid, i, c, buf);
if(ret != NC_NOERR) return "Can't write " + name + " variable in the netcdf file";
return "";
}
GradMethods::Matrix::Matrix(const std::vector<GradMethods::DataType>& in, size_t nx_, size_t ny_, struct GradMethods::MinMax minmax): nx(nx_), ny(ny_), data(nx_ * ny_)
{
if(minmax.automin || minmax.automax)
{
DataType min = in[0];
DataType max = in[0];
for(size_t i = 1; i < in.size(); i++)
if(in[i] != minmax.fill)
{
min = std::min(min, in[i]);
max = std::max(max, in[i]);
}
if(minmax.automin) minmax.min = min;
if(minmax.automax) minmax.max = max;
}
if(minmax.log)
{
minmax.min = michlib_internal::RealType<sizeof(DataType)>::Log(minmax.min);
minmax.max = michlib_internal::RealType<sizeof(DataType)>::Log(minmax.max);
}
DataType a = (std::numeric_limits<MDataType>::max() - 1) / (minmax.max - minmax.min);
for(size_t i = 1; i < in.size(); i++)
{
DataType v = minmax.log ? michlib_internal::RealType<sizeof(DataType)>::Log(in[i]) : in[i];
if(in[i] == minmax.fill)
data[i] = std::numeric_limits<MDataType>::max();
else if(v <= minmax.min)
data[i] = 0;
else if(v >= minmax.max)
data[i] = std::numeric_limits<MDataType>::max() - 1;
else
data[i] = static_cast<MDataType>(michlib_internal::RealType<sizeof(DataType)>::Round(a * (v - minmax.min)));
}
}
void GradMethods::Matrix::Grad()
{
std::vector<MDataType> out(data.size());
const auto bad = std::numeric_limits<MDataType>::max();
for(size_t iy = 0; iy < ny; iy++)
for(size_t ix = 0; ix < nx; ix++)
{
if(iy < 1 || ix < 1 || iy > ny - 2 || ix > nx - 2)
out[iy * nx + ix] = bad;
else if(V(ix - 1, iy - 1) == bad || V(ix, iy - 1) == bad || V(ix + 1, iy - 1) == bad || V(ix - 1, iy) == bad || V(ix, iy) == bad || V(ix + 1, iy) == bad ||
V(ix - 1, iy + 1) == bad || V(ix, iy + 1) == bad || V(ix + 1, iy + 1) == bad)
out[iy * nx + ix] = bad;
else
{
using IT = michlib::int4;
// Possible but unlikely overflow
const IT m1 = -1;
const IT m2 = -2;
const IT p1 = 1;
const IT p2 = 2;
IT gx = m1 * V(ix - 1, iy + 1) + p1 * V(ix + 1, iy + 1) + m2 * V(ix - 1, iy) + p2 * V(ix + 1, iy) + m1 * V(ix - 1, iy - 1) + p1 * V(ix + 1, iy - 1);
IT gy = m1 * V(ix - 1, iy - 1) + p1 * V(ix - 1, iy + 1) + m2 * V(ix, iy - 1) + p2 * V(ix, iy + 1) + m1 * V(ix + 1, iy - 1) + p1 * V(ix + 1, iy + 1);
auto sq = static_cast<IT>(michlib::Round(michlib::Hypot(gx, gy)));
if(sq >= bad) sq = bad - 1;
out[iy * nx + ix] = static_cast<MDataType>(sq);
}
}
data = std::move(out);
}

146
actions/actiongrad.h

@ -0,0 +1,146 @@
#pragma once
#include "actiondep.h"
#include "ncfilew.h"
#include "ncfuncs.h"
class GradMethods
{
public:
using DataType = float;
static DataType ToNum(const MString& str) { return michlib_internal::RealType<sizeof(DataType)>::String2Real(str.Buf()); }
struct MinMax
{
bool automin, automax, log;
DataType min, max;
DataType fill;
};
class Matrix;
class NCFileW;
};
class GradMethods::Matrix
{
public:
using MDataType = michlib::uint2;
private:
size_t nx, ny;
std::vector<MDataType> data;
public:
Matrix(const std::vector<DataType>& in, size_t nx_, size_t ny_, struct MinMax minmax);
void Grad();
auto Nx() const { return nx; }
auto Ny() const { return ny; }
const auto& V(size_t ix, size_t iy) const { return data[iy * nx + ix]; }
auto& V(size_t ix, size_t iy) { return data[iy * nx + ix]; }
const auto& Data() const { return data; }
};
class GradMethods::NCFileW
{
bool opened;
int ncid;
public:
NCFileW(): opened(false) {}
MString Create(const MString& name, const MString& history, const std::vector<MString>& vnames, const std::vector<MString>& lnames, const std::vector<GradMethods::DataType>& lons,
const std::vector<GradMethods::DataType>& lats, int compress);
MString WriteVar(const MString& name, const std::vector<GradMethods::Matrix::MDataType>& data);
~NCFileW() { Close(); }
void Close()
{
if(opened) nc_close(ncid);
opened = false;
}
};
template<class T>
concept GradSupported = requires(T t, const MString& vname) {
{
t.ReadVar(vname)
} -> std::same_as<std::vector<GradMethods::DataType>>;
};
ADD_ACTION(GRAD, grad, GradSupported<Source>, GradMethods);
template<class D> MString ActionGRAD::DoAction(const CLArgs& args, D& ds)
{
auto resop = ds.Open(args);
if(resop.Exist()) return "Can't open source: " + resop;
MString name = args.contains("out") ? args.at("out") : "out.nc";
MString min = args.contains("min") ? args.at("min") : "auto";
MString max = args.contains("max") ? args.at("max") : "auto";
int compress = args.contains("compress") ? args.at("compress").ToInt() : 3;
std::vector<Matrix> data;
std::vector<MString> lnames;
// Read data
for(size_t i = 0; i < ds.NVar(); i++)
{
const MString& name = ds.VarNames()[i];
const MString& lname = ds.LongNames()[i];
bool hmin = args.contains(name + "_min");
bool hmax = args.contains(name + "_max");
struct MinMax minmax;
minmax.log = args.contains(name + "_log");
if(hmin)
{
MString vmin = args.at(name + "_min");
minmax.automin = (vmin == "auto");
minmax.min = ToNum(vmin);
}
else
{
minmax.automin = (min == "auto");
minmax.min = ToNum(min);
}
if(hmax)
{
MString vmax = args.at(name + "_max");
minmax.automax = (vmax == "auto");
minmax.max = ToNum(vmax);
}
else
{
minmax.automax = (max == "auto");
minmax.max = ToNum(max);
}
minmax.fill = ds.FillVal(name);
data.emplace_back(ds.ReadVar(name), ds.Nx(), ds.Ny(), minmax);
lnames.emplace_back(lname + ", gradient");
}
NCFileW fw;
fw.Create(name, (ds.History().Exist() ? (ds.History() + "; ") : "") + args.at("_cmdline"), ds.VarNames(), lnames, ds.ReadLons(), ds.ReadLats(), compress);
for(size_t i = 0; i < ds.NVar(); i++)
{
const MString& name = ds.VarNames()[i];
data[i].Grad();
fw.WriteVar(name, data[i].Data());
}
return "";
};

7
actions/actionuv.cpp

@ -49,6 +49,13 @@ MString UVMethods::StPoints::WriteNcFile(const MString& name, const MString& his
text = "latitude";
ret = nc_put_att_text(ncid, yid, "standard_name", text.Len() + 1, text.Buf());
if(ret != NC_NOERR) return "Can't write standard_name attribute of latitude variable in the netcdf file";
text = "Longitude";
ret = nc_put_att_text(ncid, xid, "long_name", text.Len() + 1, text.Buf());
if(ret != NC_NOERR) return "Can't write long_name attribute of longitude variable in the netcdf file";
text = "Latitude";
ret = nc_put_att_text(ncid, yid, "long_name", text.Len() + 1, text.Buf());
if(ret != NC_NOERR) return "Can't write long_name attribute of latitude variable in the netcdf file";
}
ret = nc_def_var_deflate(ncid, xid, 1, 1, comp);

2
include/actiondep.h

@ -12,7 +12,7 @@ using michlib::MDateTime;
#define ADD_ACTION(actclass, actname, suptest, ...) ADD ACTION CLASS: actclass
#else
#define ADD_ACTION(actclass, actname, suptest, ...) \
class Action##actclass __VA_OPT__( : public) __VA_ARGS__ \
class Action##actclass __VA_OPT__( : protected) __VA_ARGS__ \
{ \
public: \
static constexpr const char* name = #actname; \

18
include/ncfilew.h

@ -31,7 +31,7 @@ class NCFileW
{
struct
{
int xdimid, ydimid;
int ydimid, xdimid;
};
int dimid[2];
};
@ -114,11 +114,11 @@ class NCFileW
else
{
const size_t i[2] = {0, 0};
const size_t c[2] = {data.Nx(), data.Ny()};
const size_t c[2] = {data.Ny(), data.Nx()};
float buf[c[0] * c[1]];
for(size_t ix = 0; ix < c[0]; ix++)
for(size_t iy = 0; iy < c[1]; iy++) buf[ix * c[1] + iy] = data.IsFill(ix, iy) ? fill : op(ix, iy);
for(size_t iy = 0; iy < c[0]; iy++)
for(size_t ix = 0; ix < c[1]; ix++) buf[iy * c[1] + ix] = data.IsFill(ix, iy) ? fill : op(ix, iy);
ret = nc_put_vara_float(ncid, vars[v].id, i, c, buf);
if(ret != NC_NOERR) return "Can't write " + vars[v].name + " variable in the netcdf file";
@ -179,14 +179,14 @@ class NCFileW
else if constexpr(dtype == G2V2)
{
const size_t i[2] = {0, 0};
const size_t c[2] = {data.Nx(), data.Ny()};
const size_t c[2] = {data.Ny(), data.Nx()};
float bufx[c[0] * c[1]];
float bufy[c[0] * c[1]];
for(size_t ix = 0; ix < c[0]; ix++)
for(size_t iy = 0; iy < c[1]; iy++)
for(size_t iy = 0; iy < c[0]; iy++)
for(size_t ix = 0; ix < c[1]; ix++)
{
bufx[ix * c[1] + iy] = data.Lon(ix, iy);
bufy[ix * c[1] + iy] = data.Lat(ix, iy);
bufx[iy * c[1] + ix] = data.Lon(ix, iy);
bufy[iy * c[1] + ix] = data.Lat(ix, iy);
}
ret = nc_put_vara_float(ncid, xid, i, c, bufx);
if(ret != NC_NOERR) return "Can't write longitude variable in the netcdf file";

129
sources/TSCDATA.cpp

@ -0,0 +1,129 @@
#define MICHLIB_NOSOURCE
#include "TSCDATA.h"
MString TSCDATAData::Open(const CLArgs& args)
{
if(!args.contains("dataset")) return "path to data not specified";
MString dataset = args.at("dataset");
michlib::NCFileA newnc;
std::vector<MString> newvnames, newlnames;
MString newhistory;
newnc.Reset(dataset);
if(!newnc) return "Can't open file " + dataset;
auto head = newnc.Header();
if(head.Dimensions().size() != 2) return "Unsupported number of dimensions";
if((head.Dimensions()[0].Name() != "longitude" || head.Dimensions()[1].Name() != "latitude") &&
(head.Dimensions()[1].Name() != "longitude" || head.Dimensions()[0].Name() != "latitude"))
return "Unsupported dimensions names";
if(head.Dimensions()[0].Name() == "longitude")
{
nx = head.Dimensions()[0].Len();
ny = head.Dimensions()[1].Len();
}
else
{
ny = head.Dimensions()[0].Len();
nx = head.Dimensions()[1].Len();
}
{
bool lonfound = false, latfound = false;
for(const auto& v: head.Variables())
if(v.Dimensions().size() == 1 && v.Type().Id() == NC_FLOAT)
{
lonfound = lonfound || (v.Dimensions()[0].Name() == "longitude" && v.Name() == "longitude");
latfound = latfound || (v.Dimensions()[0].Name() == "latitude" && v.Name() == "latitude");
}
if(!lonfound) return "Longitude not found";
if(!latfound) return "Latitude not found";
}
for(const auto& v: head.Variables())
{
if(v.Dimensions().size() != 2) continue;
if(v.Type().Id() != NC_FLOAT) continue;
if((v.Dimensions()[0].Name() != "longitude" || v.Dimensions()[1].Name() != "latitude") && (v.Dimensions()[1].Name() != "longitude" || v.Dimensions()[0].Name() != "latitude"))
continue;
newvnames.push_back(v.Name());
auto lname = newnc.A<MString>(v.Name(), "long_name");
newlnames.push_back(lname ? lname.Get() : "");
}
if(newvnames.size() == 0) return "No variables found";
{
auto his = newnc.A<MString>("history");
if(his) newhistory = his;
}
history = std::move(newhistory);
vnames = std::move(newvnames);
lnames = std::move(newlnames);
nc = std::move(newnc);
return "";
}
MString TSCDATAData::Info() const
{
if(!nc) return "";
MString out;
out += MString("Dimensions: ") + nx + " X " + ny + "\n";
out += MString("Variables: ");
for(size_t i = 0; i < vnames.size(); i++) out += ((i == 0) ? "" : ", ") + vnames[i];
out += "\n";
auto his = nc.A<MString>("history");
if(his) out += "Creator: " + his.Get() + "\n";
return out;
}
std::vector<TSCDATAData::DataType> TSCDATAData::ReadLons() const
{
std::vector<DataType> out;
auto lons = nc.V<DataType>("longitude");
if(lons)
{
out.resize(lons.DimLen(0));
for(size_t i = 0; i < out.size(); i++) out[i] = lons(i);
}
return out;
}
std::vector<TSCDATAData::DataType> TSCDATAData::ReadLats() const
{
std::vector<DataType> out;
auto lats = nc.V<DataType>("latitude");
if(lats)
{
out.resize(lats.DimLen(0));
for(size_t i = 0; i < out.size(); i++) out[i] = lats(i);
}
return out;
}
std::vector<TSCDATAData::DataType> TSCDATAData::ReadVar(const MString& name) const
{
std::vector<DataType> out;
bool havevar = false;
for(size_t i = 0; i < vnames.size(); i++) havevar = havevar || (vnames[i] == name);
if(!havevar) return out;
auto var = nc.V<DataType>(name, "latitude", "longitude");
if(var)
{
out.resize(nx * ny);
for(size_t iy = 0; iy < ny; iy++)
for(size_t ix = 0; ix < nx; ix++) out[iy * nx + ix] = var(iy, ix);
}
return out;
}
TSCDATAData::DataType TSCDATAData::FillVal(const MString& name) const
{
auto fill = nc.A<DataType>(name, "_FillValue");
return fill ? fill.Get() : 0.0;
}

35
sources/TSCDATA.h

@ -0,0 +1,35 @@
#pragma once
#include "DataAdapters/ncfilealt.h"
#include "ParseArgs.h"
#include "simple2ddata.h"
class TSCDATAData
{
michlib::NCFileA nc;
std::vector<MString> vnames, lnames;
size_t nx, ny;
MString history;
using DataType = float;
public:
static constexpr const char* name = "TSCDATA";
MString Info() const;
MString Open(const CLArgs& args);
std::vector<DataType> ReadLons() const;
std::vector<DataType> ReadLats() const;
std::vector<DataType> ReadVar(const MString& name) const;
const auto& VarNames() const { return vnames; }
const auto& LongNames() const { return lnames; }
size_t Nx() const { return nx; }
size_t Ny() const { return ny; }
size_t NVar() const { return vnames.size(); }
const auto& History() const { return history; }
DataType FillVal(const MString& name) const;
};

12
src/ncfilew.cpp

@ -10,12 +10,17 @@ MString NCFileW::CreateFile(NCFileW::Type stype, const MString& name, const MStr
int ret;
MString text;
const float node_offset = 0.0;
ret = nc_create(name.Buf(), NC_CLOBBER | NC_NETCDF4, &ncid);
if(ret != NC_NOERR) return "Can't create netcdf file: " + name;
ret = nc_put_att_text(ncid, NC_GLOBAL, "history", history.Len() + 1, history.Buf());
if(ret != NC_NOERR) return "Can't write history attribute in the netcdf file";
ret = nc_put_att(ncid, NC_GLOBAL, "node_offset", NC_FLOAT, 1, &node_offset);
if(ret != NC_NOERR) return "Can't write node_offset attribute in the netcdf file";
switch(stype)
{
case(G1V1):
@ -67,6 +72,13 @@ MString NCFileW::CreateFile(NCFileW::Type stype, const MString& name, const MStr
ret = nc_put_att_text(ncid, yid, "standard_name", text.Len() + 1, text.Buf());
if(ret != NC_NOERR) return "Can't write standard_name attribute of latitude variable in the netcdf file";
text = "Longitude";
ret = nc_put_att_text(ncid, xid, "long_name", text.Len() + 1, text.Buf());
if(ret != NC_NOERR) return "Can't write long_name attribute of longitude variable in the netcdf file";
text = "Latitude";
ret = nc_put_att_text(ncid, yid, "long_name", text.Len() + 1, text.Buf());
if(ret != NC_NOERR) return "Can't write long_name attribute of latitude variable in the netcdf file";
type = stype;
return "";
}

Loading…
Cancel
Save