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.
 
 
 

411 lines
12 KiB

#pragma once
#include "MString.h"
#include "traits.h"
#include <algorithm>
#include <netcdf.h>
#include <utility>
#include <vector>
using michlib::MString;
class NCFileWBase
{
template<class T, class Dummy> struct NCTypeD
{
static constexpr nc_type nc = NC_NAT;
};
template<class Dummy> struct NCTypeD<michlib::int1, Dummy>
{
static constexpr nc_type nc = NC_BYTE;
};
template<class Dummy> struct NCTypeD<michlib::uint1, Dummy>
{
static constexpr nc_type nc = NC_UBYTE;
static int put_var(int nc, int vid, const michlib::uint1* data) { return nc_put_var_ubyte(nc, vid, data); }
};
template<class Dummy> struct NCTypeD<michlib::int2, Dummy>
{
static constexpr nc_type nc = NC_SHORT;
static int put_var(int nc, int vid, const michlib::int2* data) { return nc_put_var_short(nc, vid, data); }
};
template<class Dummy> struct NCTypeD<michlib::uint2, Dummy>
{
static constexpr nc_type nc = NC_USHORT;
static int put_var(int nc, int vid, const michlib::uint2* data) { return nc_put_var_ushort(nc, vid, data); }
};
template<class Dummy> struct NCTypeD<michlib::int4, Dummy>
{
static constexpr nc_type nc = NC_INT;
static int put_var(int nc, int vid, const michlib::int4* data) { return nc_put_var_int(nc, vid, data); }
};
template<class Dummy> struct NCTypeD<michlib::uint4, Dummy>
{
static constexpr nc_type nc = NC_UINT;
static int put_var(int nc, int vid, const michlib::uint4* data) { return nc_put_var_uint(nc, vid, data); }
};
template<class Dummy> struct NCTypeD<michlib::int8, Dummy>
{
static constexpr nc_type nc = NC_INT64;
};
template<class Dummy> struct NCTypeD<michlib::uint8, Dummy>
{
static constexpr nc_type nc = NC_UINT64;
};
template<class Dummy> struct NCTypeD<float, Dummy>
{
static constexpr nc_type nc = NC_FLOAT;
static int put_var(int nc, int vid, const float* data) { return nc_put_var_float(nc, vid, data); }
};
template<class Dummy> struct NCTypeD<double, Dummy>
{
static constexpr nc_type nc = NC_DOUBLE;
static int put_var(int nc, int vid, const double* data) { return nc_put_var_double(nc, vid, data); }
};
template<class Dummy> struct NCTypeD<MString, Dummy>
{
static constexpr nc_type nc = NC_CHAR;
};
template<class Dummy> struct NCTypeD<char*, Dummy>
{
static constexpr nc_type nc = NC_CHAR;
};
template<class T> using NCType = NCTypeD<T, void>;
public:
template<class T> static constexpr nc_type Type2NCType = NCType<T>::nc;
// Error class
class Error
{
int err;
public:
constexpr Error(): err(NC_NOERR) {}
constexpr Error(int n_err): err(n_err) {}
explicit operator bool() const { return err == NC_NOERR; }
const char* ErrMessage() const { return nc_strerror(err); }
int Err() const { return err; }
bool IsErr() const { return err != NC_NOERR; }
void Reset() { err = NC_NOERR; }
void Reset(int n_err) { err = n_err; }
void Reset(const Error& e) { err = e.err; }
};
~NCFileWBase() { Close(); }
int Id() const { return ncid; }
int Err() const { return err.Err(); }
const char* ErrMessage() const { return err.ErrMessage(); }
explicit operator bool() const { return Ok(); }
bool Ok() const { return opened && !err.IsErr(); }
operator int() const { return ncid; }
void Close()
{
if(opened) nc_close(ncid);
err.Reset();
opened = false;
}
void Open(const MString& path)
{
if(opened) return;
err.Reset(nc_create(path.Buf(), NC_CLOBBER | NC_NETCDF4, &ncid));
opened = !err.IsErr();
}
void AddDim(const MString& name, size_t len)
{
if(err.IsErr()) return;
static int dimid;
err.Reset(nc_def_dim(ncid, name.Buf(), len, &dimid));
}
template<class T>
requires(Type2NCType<T> != NC_NAT)
void AddAtt(const MString& name, const T& val)
{
if(err.IsErr()) return;
static constexpr auto nct = Type2NCType<T>;
static constexpr bool ischar = std::is_same_v<T, char*>;
if constexpr(nct != NC_CHAR)
err.Reset(nc_put_att(ncid, NC_GLOBAL, name.Buf(), nct, 1, &val));
else if constexpr(ischar)
err.Reset(nc_put_att_text(ncid, NC_GLOBAL, name.Buf(), strlen(val) + 1, val));
else
err.Reset(nc_put_att_text(ncid, NC_GLOBAL, name.Buf(), val.Len() + 1, val.Buf()));
}
template<class T>
requires(Type2NCType<T> != NC_NAT)
void AddAtt(const MString& vname, const MString& name, const T& val)
{
if(err.IsErr()) return;
static constexpr auto nct = Type2NCType<T>;
static constexpr bool ischar = std::is_same_v<T, char*>;
int varid;
err.Reset(nc_inq_varid(ncid, vname.Buf(), &varid));
if(err.IsErr()) return;
if constexpr(nct != NC_CHAR)
err.Reset(nc_put_att(ncid, varid, name.Buf(), nct, 1, &val));
else if constexpr(ischar)
err.Reset(nc_put_att_text(ncid, varid, name.Buf(), strlen(val) + 1, val));
else
err.Reset(nc_put_att_text(ncid, varid, name.Buf(), val.Len() + 1, val.Buf()));
}
void AddAtt(const MString& vname, const MString& name, const char* val)
{
if(err.IsErr()) return;
int varid;
err.Reset(nc_inq_varid(ncid, vname.Buf(), &varid));
if(err.IsErr()) return;
err.Reset(nc_put_att_text(ncid, varid, name.Buf(), strlen(val) + 1, val));
}
template<class... T> void AddVar(const MString& vname, nc_type type, const T&... dnames) { AddVar(vname, type, std::vector<MString>({dnames...})); }
void EndDef() const { nc_enddef(ncid); }
void SetComp(const MString& vname, int comp)
{
if(err.IsErr()) return;
int varid;
err.Reset(nc_inq_varid(ncid, vname.Buf(), &varid));
if(err.IsErr()) return;
err.Reset(nc_def_var_deflate(ncid, varid, 1, 1, comp));
}
template<class T>
requires requires(int nc, int vid, const T* d) {
{
NCTypeD<T, void>::put_var(nc, vid, d)
} -> std::same_as<int>;
}
void WriteVar(const MString& vname, const T* data)
{
if(err.IsErr()) return;
int varid;
err.Reset(nc_inq_varid(ncid, vname.Buf(), &varid));
if(err.IsErr()) return;
err.Reset(NCTypeD<T, void>::put_var(ncid, varid, data));
}
private:
// Members
int ncid;
Error err = NC_NOERR;
bool opened = false;
void AddVar(const MString& vname, nc_type type, std::vector<MString>&& dnames)
{
if(err.IsErr()) return;
std::vector<int> dimids(dnames.size());
int varid;
for(size_t i = 0; i < dnames.size(); i++)
{
err.Reset(nc_inq_dimid(ncid, dnames[i].Buf(), &dimids[i]));
if(err.IsErr()) return;
}
err.Reset(nc_def_var(ncid, vname.Buf(), type, dimids.size(), dimids.data(), &varid));
}
};
class NCFileW: public NCFileWBase
{
enum Type
{
UNKNOWN,
PSET,
GPSET,
GRID,
GGRID,
RGRID,
GRGRID
};
static constexpr auto fill = std::numeric_limits<float>::max();
Type type = UNKNOWN;
int compress;
template<class D> static constexpr Type DetType()
{
if constexpr(ReadIs2DGeoRectArray<D>) return GRGRID;
if constexpr(ReadIs2DGeoArray<D>) return GGRID;
if constexpr(ReadIs1DGeoArray<D>) return GPSET;
if constexpr(ReadIs2DXYRectArray<D>) return RGRID;
if constexpr(ReadIs2DXYArray<D>) return GRID;
if constexpr(ReadIs1DArray<D>) return PSET;
return UNKNOWN;
};
template<class D> static constexpr bool Is1DType() { return Is1DType(DetType<D>()); }
static constexpr bool IsGeoType(Type t) { return t == GPSET || t == GGRID || t == GRGRID; }
static constexpr bool Is1DType(Type t) { return t == PSET || t == GPSET; }
MString CreateFile(Type stype, const MString& name, const MString& history, int compression, size_t nx, size_t ny);
public:
static constexpr auto Fill() { return fill; }
template<class D> MString Create(const D& data, const MString& name, const MString& history, int compression)
{
if(type != UNKNOWN) return "File already created";
if constexpr(Is1DType<D>())
return CreateFile(DetType<D>(), name, history, compression, data.N(), 0);
else
return CreateFile(DetType<D>(), name, history, compression, data.Nx(), data.Ny());
}
void Close()
{
NCFileWBase::Close();
type = UNKNOWN;
}
MString AddVariable(const MString& name, const MString& stname, const MString& lname, const MString& units, const MString& comment);
template<class D, class Op> MString WriteVariable(const D& data, const MString& name, Op op)
{
static constexpr auto dtype = DetType<D>();
if(type == UNKNOWN) return "File not open";
if(type != dtype) return "Incompatible data type";
if constexpr(dtype == UNKNOWN) return "Unknown data type";
EndDef();
if constexpr(Is1DType(dtype))
{
const size_t c = data.N();
float buf[c];
for(size_t ix = 0; ix < c; ix++) buf[ix] = data.IsFill(ix) ? fill : op(ix);
WriteVar(name, buf);
}
else
{
const size_t c[2] = {data.Ny(), data.Nx()};
float buf[c[0] * c[1]];
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);
WriteVar(name, buf);
}
if(!*this) return "Can't write variable " + name + ": " + ErrMessage();
return "";
}
template<class D> MString WriteVariable(const D& data, const MString& name)
{
if constexpr(Is1DType(DetType<D>()))
return WriteVariable(data, name, [&data = std::as_const(data)](size_t i) { return data(i); });
else
return WriteVariable(data, name, [&data = std::as_const(data)](size_t i, size_t j) { return data(i, j); });
}
template<class D> MString WriteGrid(const D& data)
{
static constexpr auto dtype = DetType<D>();
if(type == UNKNOWN) return "File not open";
if(type != dtype) return "Incompatible data type";
EndDef();
if constexpr(dtype == UNKNOWN)
return "Unknown data type";
else if constexpr(dtype == PSET)
{
const size_t c = data.N();
float bufx[c];
float bufy[c];
for(size_t ix = 0; ix < c; ix++)
{
bufx[ix] = data.X(ix);
bufy[ix] = data.Y(ix);
}
WriteVar("x", bufx);
WriteVar("y", bufy);
}
else if constexpr(dtype == GPSET)
{
const size_t c = data.N();
float bufx[c];
float bufy[c];
for(size_t ix = 0; ix < c; ix++)
{
bufx[ix] = data.Lon(ix);
bufy[ix] = data.Lat(ix);
}
WriteVar("longitude", bufx);
WriteVar("latitude", bufy);
}
else if constexpr(dtype == RGRID)
{
const size_t cx = data.Nx(), cy = data.Ny();
float bufx[cx];
float bufy[cy];
for(size_t ix = 0; ix < cx; ix++) bufx[ix] = data.Ix2X(ix);
for(size_t iy = 0; iy < cy; iy++) bufy[iy] = data.Iy2Y(iy);
WriteVar("x", bufx);
WriteVar("y", bufy);
}
else if constexpr(dtype == GRGRID)
{
const size_t cx = data.Nx(), cy = data.Ny();
float bufx[cx];
float bufy[cy];
for(size_t ix = 0; ix < cx; ix++) bufx[ix] = data.Ix2Lon(ix);
for(size_t iy = 0; iy < cy; iy++) bufy[iy] = data.Iy2Lat(iy);
WriteVar("longitude", bufx);
WriteVar("latitude", bufy);
}
else if constexpr(dtype == GRID)
{
const size_t c[2] = {data.Ny(), data.Nx()};
float bufx[c[0] * c[1]];
float bufy[c[0] * c[1]];
for(size_t iy = 0; iy < c[0]; iy++)
for(size_t ix = 0; ix < c[1]; ix++)
{
bufx[iy * c[1] + ix] = data.X(ix, iy);
bufy[iy * c[1] + ix] = data.Y(ix, iy);
}
WriteVar("x", bufx);
WriteVar("y", bufy);
}
else if constexpr(dtype == GGRID)
{
const size_t c[2] = {data.Ny(), data.Nx()};
float bufx[c[0] * c[1]];
float bufy[c[0] * c[1]];
for(size_t iy = 0; iy < c[0]; iy++)
for(size_t ix = 0; ix < c[1]; ix++)
{
bufx[iy * c[1] + ix] = data.Lon(ix, iy);
bufy[iy * c[1] + ix] = data.Lat(ix, iy);
}
WriteVar("longitude", bufx);
WriteVar("latitude", bufy);
}
else
return "Unknown data type";
if(!*this) return MString("Can't write grid: ") + ErrMessage();
return "";
}
};