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.

220 lines
7.7 KiB

#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);
}