odm/actions/actiongrad.cpp

#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)
{
 const float node_offset = 0.0;

 auto nx = lons.size();
 auto ny = lats.size();

 const auto fill = static_cast<GradMethods::DataType>(std::numeric_limits<GradMethods::Matrix::MDataType>::max());

 // Creating
 Open(name);
 if(!*this) return "Can't create netcdf file " + name + ": " + ErrMessage();

 AddAtt("history", history);
 AddAtt("node_offset", node_offset);

 AddDim("longitude", nx);
 AddDim("latitude", ny);
 AddVar("longitude", NC_FLOAT, "longitude");
 AddVar("latitude", NC_FLOAT, "latitude");
 SetComp("longitude", compress);
 SetComp("latitude", compress);
 AddAtt("longitude", "standard_name", "longitude");
 AddAtt("longitude", "long_name", "Longitude");
 AddAtt("latitude", "standard_name", "latitude");
 AddAtt("latitude", "long_name", "Latitude");

 // Variables
 for(size_t i = 0; i < vnames.size(); i++)
 {
  AddVar(vnames[i], NC_FLOAT, "latitude", "longitude");
  SetComp(vnames[i], compress);
  if(lnames[i].Exist()) AddAtt(vnames[i], "long_name", lnames[i]);
  AddAtt(vnames[i], "_FillValue", fill);
 }

 // End definitions
 EndDef();

 // Writing lon, lat
 WriteVar("longitude", lons.data());
 WriteVar("latitude", lats.data());

 if(!*this) return "Can't set grid in the netcdf file " + name + ": " + ErrMessage();

 return "";
}

MString GradMethods::NCFileW::WriteVariable(const MString& name, const GradMethods::Matrix& data)
{
 WriteVar(name, data.Data().data());
 if(!*this) return "Can't write variable " + name + ": " + ErrMessage();

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