"area.h"
#ifndef HEADER_ICON_AREA
#define HEADER_ICON_AREA
//**************************************************************************************************
// what is needed to define an area and get a forecast (i.e. a sequence of hours with a set of
// fields per hour, vith values) over all points of the area
//
// program(m)er: theodore.andreadis@yahoo.com, Sep 2023
//**************************************************************************************************
// for namespace geo2d
#include <limits>
#include <tuple>
#include "../../../openGL/la_vector.h"
// for namespace icon
#include "../../../_so_/parse-ini.h"
#include "../icon/field.h"
#include "../icon/thin.h"
namespace geo2d
//********************************************************************************** namespace geo2d
{
//*********************************************************************************** geo2d typedefs
typedef struct { float w, e, s, n; } bounds_t;
//************************************************************************** geo2d utility functions
vector<pair<float,float>> txtToCoords(const string &txtPoints);
//**************************************************************************************************
// a convex geographical area: a geometric convex shape with a human name
//************************************************************************** class geo2d::convexArea
class convexArea : public la::convex2d_t
{
protected:
string humanName;
public:
convexArea(const vector<la::vector2d_t> &vrtcs, const string &hName="")
: la::convex2d_t(vrtcs), humanName(hName) {};
convexArea(const convexArea &oth) : la::convex2d_t(oth()), humanName(oth.humanName) {};
operator string() const { return humanName+": "+la::convex2d_t::operator string(); };
bounds_t getBounds() const;
const string &name() const { return humanName; }
};
//******************************************************************************** class geo2d::Area
class Area : public vector<geo2d::convexArea>
{
protected:
string humanName;
float lon, lat; // optional "reperesentative" point of the area, e.g. for
public: // rectangular areas auto-created around a meteogram point
Area(const string &hName="") : humanName(hName), lon(0), lat(0) {};
operator string() const;
virtual void clear() { humanName=""; vector<geo2d::convexArea>::clear(); }
const string &name() const { return humanName; }
void setName(const string &newName) { if (newName.length()) humanName=newName; }
pair<float,float> getRepresentPoint() const;
void setRepresentPoint(float lo, float la) { lon=lo; lat=la; };
geo2d::bounds_t getBounds() const;
bool addConvexArea(const vector<pair<float,float>> &lonLats, const string &hName="");
int contains(const pair<float,float> &point) const;
la::vector2d_t barycenter() const;
};
//****************************************************************************** end namespace geo2d
}
namespace icon
//*********************************************************************************** namespace icon
{
//***************************************************************************** icon::area::typedefs
namespace area
{
typedef map<int,float> HourlyField; // storage for values over a region
typedef map<string, HourlyField> HourlyFields; // all fields for an hour
typedef map<int, HourlyFields> Forecast; // a collection of fcst::Hours for a set of fields
};
//**************************************************************************************************
// an area as a vector of convex subareas. The pairs of its data (i.e. map<int,int> gridIndices)
// contain the area point coords as the index in the global clon/clat arrays of the icon grid file
// and the index (in the class itself, being a vector<geo2d::convexArea>) of the subarea it
// belongs to
//******************************************************************************* class icon::Area
class Area : public geo2d::Area
{
protected:
map<int, int> gridIndices; // maps icon-grid index to convex subarea
const h5::Dataset *clon, *clat;
bool operator ()();
virtual bool loadSubAreas(const parse::IniFile &conf);
public:
static constexpr float halfSide=0.04/icon::scale; // use this to define square-areas around
// points, when auto-creating areas
Area(const string &hName="") : geo2d::Area(hName), clon(NULL), clat(NULL) {};
virtual bool load(const string &path, const string &aName="", bool autoArea=true);
virtual bool load(const float lo, const float la);
void clear() { gridIndices.clear(); geo2d::Area::clear(); clon=clat=NULL; };
bool setCoords(const h5::Dataset *clo, const h5::Dataset *cla);
const map<int,int> &index() const { return gridIndices; };
pair<float,float> barycenter() const
{ la::vector2d_t b=geo2d::Area::barycenter(); return { b[0], b[1] }; };
vector<vector<pair<float,float>>> corners() const;
};
//**************************************************************************************************
// a functor to generate a collection of the values of a parameter set over a region, for a
// sequence of hours
//***************************************************************************** class icon::AreaFcst
class AreaFcst
{
protected:
vector<int> hours;
vector<string> params;
public:
// external refs
vector<const Area*> areas;
const Analysis &analysis;
const Resolver &pathResolver;
ThinSwitcher &thinResolver;
// local functions
AreaFcst(vector<const Area*> ar, const Analysis &ana, const Resolver &rslv, ThinSwitcher &stch)
: areas(ar), analysis(ana), pathResolver(rslv), thinResolver(stch) {};
bool set(const vector<string> &strFieldIds, const vector<int> &hrs);
void clear() { hours.clear(); params.clear(); };
vector<area::Forecast> operator()() const;
};
//******************************************************************************* end namespace icon
};
//**************************************************************************************************
#endif
"area.cpp"
#include "./area.h"
//**************************************************************************************************
// receive a string of lon-lat pairs like: "( 22.85, 38.05 ) ( 23.82, 38.3 ) ( 24.07, 38.2 )" and
// parse it as as a vector of float pairs (of lon-lat)
//***************************************************************************** geo2d::txtToCoords()
vector<pair<float,float>> geo2d::txtToCoords(const string &txtPoints)
{
vector<pair<float,float>> result; string txt(txtPoints);
const char *strRegPoint="^\\s*(\\(\\s*[^)]+\\s*\\))",
*strRegLonLat="\\(\\s*([0-9]+\\.?[0-9]+)\\s*,\\s*([0-9]+\\.?[0-9]+)\\s*\\)";
regex_t regPoint, regLonLat;
regmatch_t matchesP[2], matchesL[3]; // 2 and 3, because 1st match is the whole text
if (0!=regcomp(®Point, strRegPoint, REG_EXTENDED))
{ fprintf(stderr, "regcomp %s error\n", strRegPoint); return result; }
if (0!=regcomp(®LonLat, strRegLonLat, REG_EXTENDED))
{ fprintf(stderr, "regcomp %s error\n", strRegLonLat); regfree(®Point); return result; }
while (REG_NOMATCH!=regexec(®Point, txt.c_str(), 2, matchesP, 0))
{
if (matchesP[1].rm_so==-1) break;
string point(txt.c_str()+matchesP[1].rm_so, txt.c_str()+matchesP[1].rm_eo);
txt.erase(0, matchesP[1].rm_eo);
if (REG_NOMATCH==regexec(®LonLat, point.c_str(), 3, matchesL, 0))
{ fprintf(stderr, "error oarsing %s point\n", point.c_str()); result.clear(); break; }
string sLon(point.c_str()+matchesL[1].rm_so, point.c_str()+matchesL[1].rm_eo);
string sLat(point.c_str()+matchesL[2].rm_so, point.c_str()+matchesL[2].rm_eo);
result.push_back({atof(sLon.c_str()), atof(sLat.c_str()) });
}
regfree(®Point); regfree(®LonLat);
return result;
}
//**************************************************************************************************
// get the max extends of a convex area at the four directions
//***************************************************************** geo2d::convexArea::getBounds()
geo2d::bounds_t geo2d::convexArea::getBounds() const
{
float w=numeric_limits<float>::max(), e=-w, s=numeric_limits<float>::max(), n=-s;
for (const la::vector2d_t &v : this->operator()())
{
const float &lon(v(0)), &lat(v(1));
if (lon>e) { e=lon; }; if (lon<w) { w=lon; }; if (lat<s) { s=lat; }; if (lat>n) n=lat;
}
return geo2d::bounds_t({w, e, s, n});
}
//***************************************************************** geo2d::Area::operator string()
geo2d::Area::operator string() const
{
string result="area "+humanName+": \n";
for (const geo2d::convexArea &cArea : *this) result+=" "+string(cArea)+"\n";
result.erase(result.length()-1); // remove the last \n
return result;
}
//**************************************************************************************************
// get the max extends of the union-of-convex areas at the four directions
//******************************************************* geo2d::bounds_t geo2d::Area::getBounds()
geo2d::bounds_t geo2d::Area::getBounds() const
{
float w=numeric_limits<float>::max(), e=-w, s=numeric_limits<float>::max(), n=-s;
for (const geo2d::convexArea &cArea : *this)
{
geo2d::bounds_t aBounds=cArea.getBounds();
if (aBounds.e>e) { e=aBounds.e; }; if (aBounds.w<w) w=aBounds.w;
if (aBounds.s<s) { s=aBounds.s; }; if (aBounds.n>n) n=aBounds.n;
}
return geo2d::bounds_t({w, e, s, n});
}
//******************************************************************* geo2d::Area::addConvexArea()
bool geo2d::Area::addConvexArea(const vector<pair<float,float>> &lonLats, const string &hName)
{
vector<la::vector2d_t> vectors2d;
for (const pair<float,float> &lonLat : lonLats)
vectors2d.push_back(la::vector2d_t({ lonLat.first, lonLat.second }));
try { this->push_back(geo2d::convexArea(vectors2d, hName)); }
catch (const la::EMsg *msg)
{
fprintf(stderr, "cannot add %s to %s: %s\n", hName.c_str(), humanName.c_str(), msg->why().c_str());
return false;
}
return true;
}
//**************************************************************************************************
// see which convex subarea contains a point and return its internal index; if none, return -1
//************************************************************************ geo2d::Area::contains()
int geo2d::Area::contains(const pair<float,float> &point) const
{
la::vector2d_t p({point.first, point.second });
int i=0;
for (const geo2d::convexArea &cArea : *this) { if (cArea.contains(p)) return i; i++; }
return -1;
}
//************************************************************************ geo2d::Area::barycenter()
la::vector2d_t geo2d::Area::barycenter() const
{
vector<la::vector2d_t> allVertices;
for (const geo2d::convexArea &cArea : *this)
{
vector<la::vector2d_t> subVertices=cArea();
allVertices.insert(allVertices.end(), subVertices.begin(), subVertices.end());
}
return la::barycenter(allVertices);
}
//***************************************************************** geo2d::Area::getRepresentPoint()
pair<float,float> geo2d::Area::getRepresentPoint() const
{
if (lon==0 && lat==0) // i.e. when not explicitly set, set lazily & implicitly to barycenter
{
la::vector2d_t bary=barycenter();
const_cast<geo2d::Area*>(this)->lon=bary(0);
const_cast<geo2d::Area*>(this)->lat=bary(1);
}
return make_pair(lon,lat);
}
//*********************************************************************** icon::Area::loadSubAreas()
bool icon::Area::loadSubAreas(const parse::IniFile &conf)
{
map<string, string> convexAreas=conf.sectionPairs(humanName+"_convexAreas");
if (!convexAreas.size())
{ fprintf(stderr, "%s has no subareas defined\n", humanName.c_str()); return false; }
for (const pair<string,string> &convex : convexAreas)
{
string shName(convex.first), pointsTxt(convex.second);
fprintf(stderr, "found %s subarea %s=%s\n",
humanName.c_str(), shName.c_str(), pointsTxt.c_str());
vector<pair<float,float>> points=geo2d::txtToCoords(pointsTxt);
if (points.size()<3)
{
fprintf(stderr, "failed to parse %s subarea %s points\n", humanName.c_str(), shName.c_str());
return false;
}
for (pair<float,float> &lonLat : points) // this loop is icon-specific
{ lonLat.first/=icon::scale; lonLat.second/=icon::scale; }
if (!addConvexArea(points, shName)) { clear(); return false; }
}
return true;
}
//**************************************************************************************************
// read the area from a ini-file. It should have (at least) these lines:
//
// #in the global vars, the area's name
// humanName=attica
//
// # a "convexAreas" section with a collection of convex subarea names and their lon-lat pairs. Eg:
// [ <humanName>_convexAreas ]
// north= ( 22.85, 38.05 ) ( 23.82, 38.3 ) ( 24.07, 38.2 ) ( 23.05, 37.93 )
// south= ( 24.06, 38.2 ) ( 23.50, 38.06 ) ( 23.93, 37.66 ) ( 24.06, 37.66 )
//******************************************************************************* icon::Area::load()
bool icon::Area::load(const string &path, const string &aName, bool autoArea)
{
clear();
parse::IniFile conf;
if (!conf.load(path.c_str())) return false;
map<string,string> sect=conf.sectionPairs(aName);
const map<string,string>::const_iterator it=sect.find("humanName");
humanName=(it!=sect.end() && it->second.length())?it->second:aName;
// see if there is a reperesent point
auto itLo=sect.find("lon"), itLa=sect.find("lat");
if (itLo==sect.end() || itLa==sect.end())
fprintf(stderr, "%s has no representation point defined\n", humanName.c_str());
else
{
lon=atof(itLo->second.c_str()), lat=atof(itLa->second.c_str());
fprintf(stderr, "represent point set to (%.3f, %.3f)\n", lon, lat);
lon/=icon::scale; lat/=icon::scale; // this is icon-specific
}
// subareas were explicitly
if (loadSubAreas(conf)) return true;
if (!autoArea || itLo==sect.end() || itLa==sect.end()) return false;
// auto-creation of a single subarea
load(lon, lat);
return true;
}
//***************************************************************************bool icon::Area::load()
bool icon::Area::load(const float lo, const float la)
{
lon=lo; lat=la;
float W=lon-icon::Area::halfSide, E=lon+icon::Area::halfSide,
S=lat-icon::Area::halfSide, N=lat+icon::Area::halfSide;
addConvexArea({ {W,N}, {E,N}, {E,S}, {W,S} });
return true;
}
//************************************************************************ icon::Area::setCoords()
bool icon::Area::setCoords(const h5::Dataset *clo, const h5::Dataset *cla)
{
clon=clat=NULL; gridIndices.clear();
if (!clo || !cla)
{ fprintf(stderr, "error: clon clat set to empty\n"); return false; }
if (clo->shape!=cla->shape)
{ fprintf(stderr, "error: area search with clon clat of different shapes\n"); return false; }
if (clo->shape.rank!=1)
{ fprintf(stderr, "error: area search with multidimnsional clon/clat\n"); return false; }
clon=clo; clat=cla;
return this->operator()();
}
//********************************************************************** icon::Area::operator ()()
bool icon::Area::operator()()
{
gridIndices.clear();
if (!clon || !clat) { fprintf(stderr, "error: clon clat set to empty\n"); return false; }
float *lon=(float*)clon->bytes, *lat=(float*)clat->bytes;
hsize_t count=clon->elemCount();
geo2d::bounds_t bounds=getBounds();
fprintf(stderr, "%s bound in -w %.3f -e %.3f -s %.3f -n %.3f, testing %llu points \n",
humanName.c_str(), bounds.w*icon::scale, bounds.e*icon::scale,
bounds.s*icon::scale, bounds.n*icon::scale, count);
for (hsize_t i=0; i<count; i++,lon++,lat++)
{
// avoid barocentric calcs for most cases
if (*lon<bounds.w || *lon>bounds.e) { continue; }; if (*lat<bounds.s || *lat>bounds.n) continue;
int convexAreaNo=contains( {*lon, *lat }); if (convexAreaNo<0) continue;
gridIndices[i]=convexAreaNo;
}
return true;
}
//**************************************************************************** icon::Area::corners()
vector<vector<pair<float,float>>> icon::Area::corners() const
{
vector<vector<pair<float,float>>> result; result.reserve(this->size());
for (vector<geo2d::convexArea>::const_iterator it=this->begin(); it!=this->end(); it++)
{
const vector<la::vector2d_t> &vertices=it->operator()();
result.push_back(vector<pair<float,float>>()); vector<pair<float,float>> ¤t=result.back();
for (const la::vector2d_t &v : vertices) current.push_back(make_pair(v(0), v(1)));
}
return result;
}
//******************************************************************* icon::AreaFcst::setHourRange()
bool icon::AreaFcst::set(const vector<string> &strFieldIds, const vector<int> &hrs)
{
clear(); hours=hrs;
for (const string &strId : strFieldIds)
{
string id;
try { id=string(icon::FieldId(strId)); }
catch (const ios_base::failure *err)
{ fprintf(stderr, "skip misformed param %s\n", strId.c_str()); continue; }
vector<string>::const_iterator it=find(params.begin(), params.end(), id);
if (it!=params.end())
{ fprintf(stderr, "param %s given more than one times\n", strId.c_str()); continue; }
params.push_back(id);
}
return params.size();
}
//********************************************************************* icon::AreaFcst::operator()()
vector<icon::area::Forecast> icon::AreaFcst::operator()() const
{
// prepare each slot of the result
vector<icon::area::Forecast> result;
result.reserve(areas.size());
for (long unsigned a=0; a<areas.size(); a++)
{
result.push_back(icon::area::Forecast());
for (int hour : hours) result.back()[hour]=icon::area::HourlyFields();
}
for (int hour : hours)
for (const string &strId : params)
{
// read the param's full-area values for this hour
icon::FieldId fieldId(strId); icon::FcstHour fcstHour(analysis, hour);
string path=pathResolver(fcstHour, icon::H5NameToLevType(fieldId.getLevName()));
fprintf(stderr, "reading %s, for %s, from %s\n",
strId.c_str(), string(fcstHour).c_str(), path.c_str());
icon::File dataFile(path);
if (!dataFile.isValid()) { fprintf(stderr, "unable to use %s\n", path.c_str()); continue; }
h5::Dataset param=dataFile.getVals(fieldId.getName());
if (!param.isValid())
{ fprintf(stderr, "unable to read %s from %s\n", strId.c_str(), path.c_str()); continue; }
if (!thinResolver(fieldId.getName(), dataFile))
{ fprintf(stderr, "error thinning of %s\n", strId.c_str()); continue; }
int lNo=fieldId.getLevNo();
if (fieldId.getLevName()==icon::levTypeToH5Name(icon::LevType::p)) lNo*=100;
const float *values=dataFile.floats(param, fieldId.getLevName(), lNo);
if (!values)
{ fprintf(stderr, "read error of %s\n", strId.c_str()); continue; }
icon::floatTransformFunc_t *trFunc=icon::getTransform(fieldId.getName());
// use this param's values read for all areas to fill each area's this hourly forecast
for (long unsigned a=0; a<areas.size(); a++)
{
result[a][hour][strId]=icon::area::HourlyField();
for (const pair<int,int> &indx: areas[a]->index())
{
int i=thinResolver.toThin(indx.first); if (i<0) continue; // point not mapped: for radius thining
result[a][hour][strId][indx.first]=trFunc(values[i]);
}
fprintf(stderr, "found %lu values over [%s][%d][%s]\n", result[a][hour][strId].size(),
areas[a]->name().c_str(), hour, strId.c_str());
}
}
fprintf(stderr, "read: \n");
for (long unsigned a=0; a<result.size(); a++)
for (auto h : hours)
for (auto p : params)
fprintf(stderr, " forecast[%s][%d hrs][%s] %lu values\n",
areas[a]->name().c_str(), h, p.c_str(), result[a][h][p].size());
return result;
}
"iconAreaGrid-gmt.h"
#ifndef HEADER_ICON_AREA_GRID_GMT
#define HEADER_ICON_AREA_GRID_GMT
//**************************************************************************************************
// the print functions to stderr and gmt-style textfile the iconAreaGrid-main.cpp uses
//**************************************************************************************************
#include "../icon/grid.h"
bool printRepresentPoint(const pair<float,float> &point, const string &name="", const string &outPath="");
bool printCenters(const vector<pair<float,float>> &points, const string &outPath="");
bool printTriangles(const vector<pair<float,float>> &points, const string &outPath="");
bool printCorners(const vector<vector<pair<float,float>>> &corners, const string &outPath="");
//**************************************************************************************************
#endif
"iconAreaGrid-gmt.cpp"
#include "./iconAreaGrid-gmt.h"
//**************************************************************************** printRepresentPoint()
bool printRepresentPoint(const pair<float,float> &point, const string &name, const string &outPath)
{
FILE *fout=stderr;
if (outPath.length())
{
fout=fopen(outPath.c_str(), "w");
if (fout==NULL) { fprintf(stderr, "%s: %s\n", outPath.c_str(), strerror(errno)); return false; }
}
if (fout==stderr)
fprintf(fout, "represent point: (%.4f, %.4f)\n", point.first*icon::scale, point.second*icon::scale);
else
fprintf(fout, "%.4f %.4f %s\n", point.first*icon::scale, point.second*icon::scale, name.c_str());
return true;
}
//*********************************************************************************** printCenters()
bool printCenters(const vector<pair<float,float>> &points, const string &outPath)
{
FILE *fout=stderr;
if (outPath.length())
{
fout=fopen(outPath.c_str(), "w");
if (fout==NULL) { fprintf(stderr, "%s: %s\n", outPath.c_str(), strerror(errno)); return false; }
}
if (fout==stderr) fprintf(stderr, "the area's %lu points are:\n", points.size());
for (const pair<float,float> &point: points)
fprintf(fout, "%.4f %.4f\n", point.first*icon::scale, point.second*icon::scale);
if (outPath.length()) fclose(fout);
return true;
}
//********************************************************************************* printTriangles()
bool printTriangles(const vector<pair<float,float>> &points, const string &outPath)
{
FILE *fout=NULL;
if (outPath.length())
{
fout=fopen(outPath.c_str(), "w");
if (fout==NULL) { fprintf(stderr, "%s: %s\n", outPath.c_str(), strerror(errno)); return false; }
}
int areaNo=0;
if (fout==NULL) fprintf(stderr, "the area's %lu facets are:\n", points.size()/3);
vector<pair<float,float>>::const_iterator it=points.begin();
while (it!=points.end())
{
if (fout==NULL) fprintf(stderr, "the area's facet no%d vertices are:\n", areaNo++);
else fprintf(fout, "path\n");
vector<pair<float,float>>::const_iterator it0=it;
for (int i=0; i<3; i++, it++)
{
if (fout==NULL)
{
fprintf(stderr, " vertex[%d]: ", i);
fprintf(stderr, "%.4f %.4f\n", it->first*icon::scale, it->second*icon::scale);
}
else fprintf(fout, "%.4f %.4f\n", it->first*icon::scale, it->second*icon::scale);
}
if (fout) // for gmt text file, repeat the 1st point to close the path
fprintf(fout, "%.4f %.4f\n", it0->first*icon::scale, it0->second*icon::scale);
}
if (fout) fclose(fout);
return true;
}
//*********************************************************************************** printCorners()
bool printCorners(const vector<vector<pair<float,float>>> &corners, const string &outPath)
{
FILE *fout=NULL;
if (outPath.length())
{
fout=fopen(outPath.c_str(), "w");
if (fout==NULL) { fprintf(stderr, "%s: %s\n", outPath.c_str(), strerror(errno)); return false; }
}
int convexNo=0;
for (const vector<pair<float,float>> &convex : corners)
{
if (convex.size()<3) continue;
if (!fout) fprintf(stderr, "subarea No%d corners: \n ", convexNo++);
else fprintf(fout, "path\n");
for (const pair<float,float> &point : convex)
if (!fout) fprintf(stderr, " (%.4f, %.4f) ", point.first*icon::scale, point.second*icon::scale);
else fprintf(fout, "%.4f %.4f\n", point.first*icon::scale, point.second*icon::scale);
// for gmt repeat the initial point to close path
if (!fout) fprintf(stderr, "\n");
else fprintf(fout, "%.4f %.4f\n", convex[0].first*icon::scale, convex[0].second*icon::scale);
}
if (fout) fclose(fout);
return true;
}
"reducer.h"
#ifndef HEADER_ICON_AREA_REDUCER
#define HEADER_ICON_AREA_REDUCER
//**************************************************************************************************
// what is needed to reduce area forecasts to single-point meteograms
//
// program(m)er: theodore.andreadis@yahoo.com, Sep 2023
//**************************************************************************************************
#include <limits>
#include "../meteograms/meteograms.h"
#include "./area.h"
namespace icon
//*********************************************************************************** namespace icon
{
//***************************************************************************** icon::area::typedefs
// from namespace area:
// typedef map<int,float> HourlyField; // storage for values over a region
// typedef map<string, HourlyField> HourlyFields; // all fields for an hour
// typedef map<int, HourlyFields> Forecast; // a collection of fcst::Hours for a set of fields
//***************************************************************************** namespace icon::area
// functors that reduce a set of area values to a single value
namespace area
{
//**************************************************************************************************
// typical icon::area::IReduce functors
//************************************************************************ class IReduce && children
constexpr int defaultLowBound=-273; // the value below which found param values
// found in a reducer, are considered "lowBound"
//******************************* class icon::area::IReduce, a common ancestor for all area reducers
class IReduce
{
protected:
int lowBound; // any value in a field lower than that, will be ignored
public: // while reducing. Typically also used as "error" return-value
IReduce(const int lowB=defaultLowBound) : lowBound(lowB) {};
// receive a field and do a custom reducing it to a single value
virtual float operator()(const area::HourlyField &) { return defaultLowBound; }
// an human-readable identifier of the reducer classtype
virtual string name() const { return "invalid"; };
// if the single value of the reduced field will need a new param name
virtual string newNameOf(const string ¶mName) const { return paramName; };
// for when internal state has to be cleared before nother loop
virtual void reset() {};
// if you want to know if a reducted value was found
virtual bool valid(const float value) const { return value>lowBound; };
};
//********************************************* class icon::area::ReduceAvg, return area's avg value
class ReduceAvg : public IReduce // return area's avg value
{
public:
ReduceAvg(const int lowB=defaultLowBound) : IReduce(lowB) {};
virtual float operator()(const area::HourlyField &field);
virtual string name() const { return "avg (above "+to_string(lowBound)+")"; };
};
static ReduceAvg reduceAvg, reduceAvgM1(-1);
//********************************************* class icon::area::ReduceMax, return area's max value
class ReduceMax : public IReduce
{
public:
ReduceMax(const int lowB=defaultLowBound) : IReduce(lowB) {};
virtual float operator()(const area::HourlyField &field);
virtual string name() const { return "max (above "+to_string(lowBound)+")"; };
};
static ReduceMax reduceMax, reduceMaxM1(-1);
//********************************************* class icon::area::ReduceMin, return area's min value
class ReduceMin : public IReduce
{
public:
ReduceMin(const int lowB=defaultLowBound) : IReduce(lowB) {};
virtual float operator()(const area::HourlyField &field);
virtual string name() const { return "min (above "+to_string(lowBound)+")"; };
};
static ReduceMin reduceMin, reduceMinM1(-1);
//****************************************************************** class icon::area::ReduceCumuMax
class ReduceCumuMax : public IReduce
{
protected:
const HourlyField *prevField; float prevresult;
public:
ReduceCumuMax(const int lowB=defaultLowBound) : IReduce(lowB), prevField(NULL), prevResult(0) {};
virtual string name() const { return "ReduceCumuMax"; };
virtual float operator()(const HourlyField &field);
virtual void reset() { prevField=NULL; prevResult=0; };
};
//******************************************************************* class icon::area::ReduceMasked
class ReduceMasked : public IReduce // return area's value for a field not measured everywhere
{ // lile low-cloud ceiling, where -1 means no-low-ceiling
public:
// at least that % masked, means all field as masked; e.g. >=75% skyClear means full skyClear
static constexpr float perc=0.75;
float fillValue; // what to use for a masked piunt (i.e. cloud base for skyClear)
ReduceMasked(const float fillV=0, const int lowB=defaultLowBound)
: IReduce(lowB), fillValue(fillV) {};
virtual float operator()(const area::HourlyField &field);
virtual string name() const { return "masked (eq "+to_string(lowBound)+")"; };
};
static ReduceMasked reduceLowCeiling(2000, -1);
//**************************************************************************************************
// a baseclass 2-part field reducer (e.g. wind's u,v). It will accept the two field components and
// decide and store the reduced value for each component. When the "usual" ReduceI::operator() is
// asked for each component field, it will return the component-appropriate stored value
//********************************************************************* class icon::area::IReduceDuo
class IReduceDuo : public IReduce
{
protected:
map<const void*, const void*> uv; map<const void*, pair<float,float>> uvVals;
virtual pair<float,float> uvCalc(const HourlyField *u, const HourlyField *v)
{ return { lowBound, lowBound }; };
public:
IReduceDuo(const int lowB=defaultLowBound) : IReduce(lowB) {};
virtual float operator()(const HourlyField &field);
virtual string name() const { return "invalidDuo"; };
virtual void reset() { uv.clear(); uvVals.clear(); }
// duo-reducer specifics:
bool prepare(const FieldId &uId, const FieldId &vId, const HourlyFields *fields);
};
//**************************************************************************************************
// a two-component reducer that will reduce a field to the max-force component
//****************************************************************** class icon::area::ReduceDuoMax
class ReduceDuoMax : public IReduceDuo
{
protected:
virtual pair<float,float> uvCalc(const HourlyField *u, const HourlyField *v);
public:
ReduceDuoMax(const int lowB=defaultLowBound) : IReduceDuo(lowB) {};
virtual string name() const { return "maxDuo"; };
};
//**************************************************************************************************
// a two-component reducer that will reduce a field to the min-force component
//****************************************************************** class icon::area::ReduceDuoMin
class ReduceDuoMin : public IReduceDuo
{
protected:
virtual pair<float,float> uvCalc(const HourlyField *u, const HourlyField *v);
public:
ReduceDuoMin(const int lowB=defaultLowBound) : IReduceDuo(lowB) {};
virtual string name() const { return "minDuo"; };
};
//**************************************************************************************************
// two trivial ReduceDuoMin/Max wind reducers, that just add a human meaningful name to themserves
//************************************************************** class icon::area::ReduceWindMin/Max
class ReduceWindMax : public ReduceDuoMax
{
public:
ReduceWindMax(const int lowB=defaultLowBound) : ReduceDuoMax(lowB) {};
virtual string newNameOf(const string ¶mName) const { return "max_"+paramName; };
};
static ReduceWindMax reduceWindMax;
class ReduceWindMin : public ReduceDuoMin
{
public:
ReduceWindMin(const int lowB=defaultLowBound) : ReduceDuoMin(lowB) {};
virtual string newNameOf(const string ¶mName) const { return "min_"+paramName; };
};
static ReduceWindMin reduceWindMin;
//*************************************************************************** collection of reducers
typedef map<string, IReduce*> Reducers;
Reducers getDefaultReducers();
IReduce *findReducer(const string &fieldName, Reducers &reducers);
//**************************************************************************************************
// a functor class to convert from area forecast to point meteogram using reducers
//***************************************************************** class icon::area::FcstToMeteogram
class FcstToMeteogram
{
protected:
const icon::Area &area;
float lon, lat; int elev; // the area's meteogram "representation" point
const icon::Analysis &analysis;
public:
bool copyInvalids; // copy "invalid" reducted values to meteograms? default is no
FcstToMeteogram(const icon::Area &ar, float lo, float la, int elv, const icon::Analysis &ana)
: area(ar), lon(lo), lat(la), elev(elv), analysis(ana), copyInvalids(false) {};
// create a meteogram and fill it from a forecast
Meteogram operator()(const Forecast &forecast, Reducers &reducers) const
{ Meteogram m(area.name(), lon, lat, elev, analysis); add(m, forecast, reducers); return m; }
// add another forecast to an existing meteogram
virtual int add(Meteogram &m, const Forecast &forecast, Reducers &reducers) const;
};
//***************************************************************************** namespace icon::area
};
//******************************************************************************* end namespace icon
};
//**************************************************************************************************
#endif
"reducer.cpp"
#include "./reducer.h"
//***************************************************************************** static _firstValid()
static float _firstValid(const float &lowBound, const icon::area::HourlyField *field)
{
float result=lowBound; if (!field || !field->size()) return result;
for (icon::area::HourlyField::const_iterator it=field->begin(); it!=field->end(); it++)
if (it->second>lowBound) return it->second;
return result;
}
//************************************************************************** icon::area::ReduceAvg()
float icon::area::ReduceAvg::operator()(const icon::area::HourlyField &field)
{
if (!field.size()) return icon::area::IReduce::operator()(field);
float sum=0; unsigned count=0;
for (const pair<int, float> &iVal : field) if (iVal.second>lowBound) { count++; sum+=iVal.second; }
return (count)?sum/count:lowBound;
}
//************************************************************************** icon::area::ReduceMax()
float icon::area::ReduceMax::operator()(const icon::area::HourlyField &field)
{
float result=_firstValid(lowBound, &field); if (result==lowBound) return lowBound;
for (const pair<int, float> &iVal : field)
if (iVal.second>lowBound && iVal.second>=result) result=iVal.second;
return result;
}
//************************************************************************** icon::area::ReduceMin()
float icon::area::ReduceMin::operator()(const icon::area::HourlyField &field)
{
float result=_firstValid(lowBound, &field); if (result==lowBound) return lowBound;
for (const pair<int, float> &iVal : field)
if (iVal.second>lowBound && iVal.second<=result) result=iVal.second;
return result;
}
//************************************************************ icon::area::ReduceCumuMax::operator()
float icon::area::ReduceCumuMax::operator()(const icon::area::HourlyField &field)
{
float result=0;
// find the max difference between this and previous cummulative fields
if (prevField!=NULL)
{
fprintf(stderr, "comparing %p to new %p\n", prevField, &field);
icon::area::HourlyField::const_iterator itp=prevField->begin();
for (icon::area::HourlyField::const_iterator it=field.begin(); it!=field.end(); it++,itp++)
{
if (it->second<=lowBound || itp->second<=lowBound) continue;
float diff=it->second-itp->second; if (diff>result) result=diff;
}
}
// or just the max for the 1st time (with no previous field recorded)
else
{
result=_firstValid(lowBound, &field);
if (result==lowBound) result=0;
else
for (const pair<int, float> &iVal : field)
if (iVal.second>lowBound && iVal.second>=result) result=iVal.second;
}
// increase the cummulative reduced value by the current max-diff and save for the next field
result+=prevresult; prevResult=result; prevField=&field;
return result;
}
//************************************************************* icon::area::ReduceMasked::operator()
float icon::area::ReduceMasked::operator()(const icon::area::HourlyField &field)
{
if (!field.size()) return icon::area::IReduce::operator()(field);
float sum=0; unsigned count=0;
for (const pair<int, float> &iVal : field)
{ float val=iVal.second; if (val==lowBound) { val=fillValue; count++; }; sum+=val; }
if (count>=perc*field.size()) return fillValue;
float result=sum/field.size();
if (abs(result-fillValue)<=abs((1-perc)*fillValue)) return fillValue;
return result;
};
//***************************************************************** icon::area::IReduceDuo prepare()
bool icon::area::IReduceDuo::prepare
(const icon::FieldId &uFid, const icon::FieldId &vFid, const icon::area::HourlyFields *fields)
{
string uId=string(uFid), vId=string(vFid);
icon::area::HourlyFields::const_iterator itu=fields->find(uId), itv=fields->find(vId);
if (itu==fields->end() || itv==fields->end())
{ fprintf(stderr, "not both %s & %s coords found\n", uId.c_str(), vId.c_str()); return false; }
const icon::area::HourlyField *pU=&(itu->second), *pV=&(itv->second);
// first occurence of uId, set it up
if (uv.find((void*)pU)==uv.end())
{
//fprintf(stderr, "set %s-reducer: %s (@ %p) + %s\n", name().c_str(), uId.c_str(), pU, vId.c_str());
uv[pU]=pV; uvVals[pU]=uvCalc(pU, pV); return true;
}
// else u-component already encounterd
if (pV!=uv[(void*)pU]) // pV should already be there. This assumes
{ // u-component is ALWAYS encountered first
fprintf(stderr, "error: %s resetting already set %s\n", string(uId).c_str(), name().c_str());
return false;
}
//fprintf(stderr, "%s-reducer already prepared for %s\n", name().c_str(), vId.c_str());
return true;
}
//*************************************************************** icon::area::IReduceDuo::operator()
float icon::area::IReduceDuo::operator()(const icon::area::HourlyField &field)
{
const void *p=(void*)&field, *pU=NULL;
for (map<const void*, const void*>::const_iterator it=uv.begin(); it!=uv.end(); it++)
if (p==it->first) { pU=p; break; } else if (p==it->second) { pU=it->first; break; }
if (!pU) return lowBound;
map<const void*, pair<float,float>>::const_iterator it=uvVals.find(pU);
if (it==uvVals.end()) return lowBound;
float result=(p==pU)?it->second.first:it->second.second;
if (p!=pU) // that is, we answered for the 2nd coord, assume 1st has already been answered
{
uv.erase(pU); uvVals.erase(pU);
//fprintf(stderr, "%s-reducing done with %p\n", name().c_str(), pU);
}
return result;
}
//*************************************************************** icon::area::ReduceDuoMin::uvCalc()
pair<float,float> icon::area::ReduceDuoMin::uvCalc
(const icon::area::HourlyField *u, const icon::area::HourlyField *v)
{
pair<float,float> result={ _firstValid(lowBound, u), _firstValid(lowBound, v) };
if (result.first==lowBound || result.second==lowBound) return { lowBound, lowBound };
float minLen=result.first*result.first+result.second*result.second;
icon::area::HourlyField::const_iterator itu, itv;
for (itu=u->begin(); itu!=u->end(); itu++, itv++)
{
if (itu->second<=lowBound || itv->second<=lowBound) continue;
float len=itu->second*itu->second+itv->second*itv->second;
if (len<minLen) { minLen=len; result.first=itu->second; result.second=itv->second; }
}
return result;
}
//*************************************************************** icon::area::ReduceDuoMax::uvCalc()
pair<float,float> icon::area::ReduceDuoMax::uvCalc
(const icon::area::HourlyField *u, const icon::area::HourlyField *v)
{
pair<float,float> result={ _firstValid(lowBound, u), _firstValid(lowBound, v) };
if (result.first==lowBound || result.second==lowBound) return { lowBound, lowBound };
float maxLen=result.first*result.first+result.second*result.second;
icon::area::HourlyField::const_iterator itu, itv=v->begin();
for (itu=u->begin(); itu!=u->end(); itu++, itv++)
{
if (itu->second<=lowBound || itv->second<=lowBound) continue;
float len=itu->second*itu->second+itv->second*itv->second;
if (len>maxLen) { maxLen=len; result.first=itu->second; result.second=itv->second; }
}
return result;
}
//**************************************************************************************************
// get a "usual" set of field-to-single-value reducers, to be overriden/extended as each app needs
//***************************************************************** icon::area::getDefaultReducers()
icon::area::Reducers icon::area::getDefaultReducers()
{
return
{
{ "default", &icon::area::reduceAvg },
{ clch, &icon::area::reduceMax },
{ clcl, &icon::area::reduceMax },
{ clcm, &icon::area::reduceMax },
{ gust10, &icon::area::reduceMax },
{ hbas_con, &icon::area::reduceMinM1 },
{ h_snow, &icon::area::reduceMax },
{ htop_con, &icon::area::reduceMaxM1 },
{ hzerocl, &icon::area::reduceMinM1 },
{ pres, &icon::area::reduceAvg },
{ pres_msl, &icon::area::reduceAvg },
{ pres_sfc, &icon::area::reduceAvg },
{ rain_con, &icon::area::reduceMax },
{ rh_2m, &icon::area::reduceAvg },
{ snow_con, &icon::area::reduceMax },
{ snowfrac, &icon::area::reduceMax },
{ snow_gsp, &icon::area::reduceMax },
{ sob_s, &icon::area::reduceMax },
{ t_2m, &icon::area::reduceAvg },
{ td_2m, &icon::area::reduceAvg },
{ temp, &icon::area::reduceAvg },
{ tmax_2m, &icon::area::reduceMax },
{ tmin_2m, &icon::area::reduceMin },
{ tot_prec, &icon::area::reduceMax },
{ u_10m, &icon::area::reduceAvg },
{ u, &icon::area::reduceAvg },
{ v_10m, &icon::area::reduceAvg },
{ v, &icon::area::reduceAvg },
{ clc_ceiling, &icon::area::reduceLowCeiling }
};
}
//************************************************************************************ findReducer()
icon::area::IReduce *icon::area::findReducer(const string &fieldName, icon::area::Reducers &reducers)
{
icon::area::IReduce *dfltReduce=&icon::area::reduceAvg;
icon::area::Reducers::const_iterator it=reducers.find(fieldName);
if (it==reducers.end()) { it=reducers.find("default"); if (it==reducers.end()) return dfltReduce; }
return it->second;
}
//*************************************************************** icon::area::FcstToMeteogram::add()
int icon::area::FcstToMeteogram::add(icon::Meteogram &meteogram,
const icon::area::Forecast &forecast, icon::area::Reducers &reducers) const
{
int result=0;
for (pair<string, IReduce*> r : reducers) r.second->reset();
for (auto hit=forecast.begin(); hit!=forecast.end(); hit++)
{
int hour=hit->first; const icon::area::HourlyFields *hourlyFields=&(hit->second);
for (auto it=hourlyFields->begin(); it!=hourlyFields->end(); it++)
{
icon::FieldId fieldId(it->first);
const icon::area::HourlyField *field=&(it->second);
icon::area::IReduce *reducer=icon::area::findReducer(fieldId.getName(), reducers);
if (reducer==NULL)
{ fprintf(stderr, "no reducer for %s, skipping\n", fieldId.getName().c_str()); continue; }
//fprintf(stderr, "reducing %dhr %s, at %p\n", hour, string(it->first).c_str(), field);
if (fieldId.isWind() && dynamic_cast<icon::area::IReduceDuo*>(reducer))
{
string sid=string(fieldId), sidU=sid, sidV=sid; sidU[0]='u'; sidV[0]='v';
if (!((icon::area::IReduceDuo*)reducer)->prepare(sidU, sidV, hourlyFields))
{ fprintf(stderr, "failed %s reducer, skip it\n", fieldId.getName().c_str()); continue; }
}
float value=(*reducer)(*field);
fprintf(stderr, "%s-reducer reduced %s %s-field at %d hr to %.2f\n", reducer->name().c_str(),
area.name().c_str(), fieldId.getName().c_str(), hour, value);
if (copyInvalids || reducer->valid(value))
{
icon::DataPoint point(reducer->newNameOf(fieldId.getName()), value, lon, lat,
icon::H5NameToLevType(fieldId.getLevName()), fieldId.getLevNo());
meteogram.add(hour, point); result++;
}
}
}
return result;
}
"iconAreaFcst-main.cpp"
//**************************************************************************************************
// receive a region, an analysis with a hour range and a set of parameters and either print into
// text files the parameter values over the area for the hours range, or reduce that forecast to a
// single "representative" point meteogram
//
// usage: $0 [-i inBaseDir] [-o outDir|-M(eteogram] -g gridFile -p param1[,..,paramN]
// [-y year] [-m month] [-d day] [-a analysis] [-f fromHour] [-t toHour] [-s stepHour]
// [-n areaName] areaConf | -h
//
// where:
// inBaseDir : the non-date part of the icon files directory
// outDir : where the text files will be written (one per hour and parameter). Default is $CWD
// meteogram : if -M is given, outDir is ignored and instead of output text files, a reduced
// meteogram is printed to stderr
// gridFile : the path of the icon grid file
// year,mont,day,analysis : define the analysis date
// hourFrom/To/Step : the hours range we ask data for
// areConf : the area geography definition time
// areaName : the [ areaName ] section in areaConf, if it contains more than one area defs
// param1,N : a comma separated list of params. For 3d, they are given as name:levName:LevNo
// e.g -p t_2m,temp:plev:850 where plev dataset appear in the pressure-levels file
//
// program(m)er: theodore.andreadis@yahoo.com, Sep 2023
//**************************************************************************************************
#include "reducer.h"
//************************************************************************************* getCmdLine()
bool getCmdLine(int argC, char *args[],
icon::Analysis &analysis, vector<int> &hours, string &gridPath, string &inBaseDir,
string &outDir, string &areaName, string &areaConf, vector<string> ¶ms)
{
constexpr char *usage=
(char*)"usage: %s [-i inBaseDir] [-o outDir|-M(eteogram)] -g gridFile -p param1[,..,paramN] "
"[-y year] [-m month] [-a analysis] [-f fromHour] [-t toHour] [-s stepHour] "
"[-n areaName] areaConf | -h\n";
// set defaults
gridPath=areaConf=areaName=""; inBaseDir=outDir="."; hours.clear();
unsigned y=analysis.y(), m=analysis.m(), d=analysis.d(), a=analysis.a();
string sParams=""; params.clear();
// read options
int hourFirst=0, hourLast=0, hourStep=1;
opterr=0; int opt;
while ((opt=getopt(argC, args, ":i:o:Mg:p:y:m:d:a:f:t:s:n:h"))!=-1)
{
if (opt=='i') inBaseDir=optarg;
else if (opt=='M') outDir="";
else if (opt=='o') { if (outDir.length()) outDir=optarg; }
else if (opt=='g') gridPath=optarg;
else if (opt=='p') sParams=optarg;
else if (opt=='y') y=atoi(optarg);
else if (opt=='m') m=atoi(optarg);
else if (opt=='d') d=atoi(optarg);
else if (opt=='a') a=atoi(optarg);
else if (opt=='f') hourFirst=atoi(optarg);
else if (opt=='t') hourLast=atoi(optarg);
else if (opt=='s') hourStep=atoi(optarg);
else if (opt=='n') areaName=optarg;
else if (opt=='h') { fprintf(stderr, usage, args[0]); return false; }
else { fprintf(stderr, "unknown option -%c, ignored\n", optopt); return false; }
}
if (argC-optind!=0) areaConf=args[optind];
hourFirst=abs(hourFirst); hourLast=abs(hourLast); hourStep=max(1, abs(hourStep));
if (hourFirst>hourLast) swap(hourFirst, hourLast);
for (int h=hourFirst; h<=hourLast; h+=hourStep) hours.push_back(h);
// do some checking
if (!hours.size())
{ fprintf(stderr, "no forecast hours given,\n"); fprintf(stderr, usage, args[0]); return false; }
if (!areaConf.length())
{ fprintf(stderr, "no area-conf path given,\n"); fprintf(stderr, usage, args[0]); return false; }
if (!gridPath.length())
{ fprintf(stderr, "no gridFile given,\n"); fprintf(stderr, usage, args[0]); return false; }
analysis=icon::Analysis(y,m,d,a);
if (analysis.t()-time(NULL)>10*365*24*3600)
{ fprintf(stderr, "error: +10years old analysis given\n"); return false; }
string::size_type idx;
while ((idx=sParams.find(','))!=string::npos)
{
string sParam=sParams.substr(0, idx); sParams.erase(0, idx+1);
if (sParam.length()) params.push_back(sParam);
}
if (sParams.length()) params.push_back(sParams);
if (!params.size())
{ fprintf(stderr, "no param names given, "); fprintf(stderr, usage, args[0]); return false; }
return true;
}
//****************************************************************************** _writeHourlyParam()
bool _writeHourlyParam(const string &path, const icon::area::HourlyField &areaParamIdx,
const float *clon, const float *clat)
{
if (!areaParamIdx.size())
{ fprintf(stderr, "nothing to write at %s\n", path.c_str()); return true; }
fprintf(stderr, "writing %lu points at %s\n", areaParamIdx.size(), path.c_str());
int fd=open(path.c_str(), O_CREAT|O_TRUNC|O_WRONLY, S_IRUSR|S_IWUSR|S_IRGRP);
if (fd==-1) { fprintf(stderr, "%s: %s\n", path.c_str(), strerror(errno)); return false; }
FILE *f=fdopen(fd, "w");
for (const pair<int,float> &p : areaParamIdx)
fprintf(f, "%.4f %.4f %.3f\n", clon[p.first]*icon::scale, clat[p.first]*icon::scale, p.second);
fclose(f); close(fd);
return true;
}
//*********************************************************************************** _txtFileName()
string _txtFileName(const string ¶m, int hour)
{
char buff[4]; memset(buff, 0, sizeof(buff));
snprintf(buff, sizeof(buff), "%3.3d", hour);
return param+"-"+string(buff)+"hr.txt";
}
//******************************************************************************************* main()
int main(int argC, char *args[])
{
//*** read cmdline opts
string gridPath, areaName, areaConf, inBaseDir, outDir;
vector<string> params; vector<int> hours; icon::Analysis analysis;
if (!getCmdLine(argC, args, analysis, hours, gridPath, inBaseDir, outDir,
areaName, areaConf, params)) return EXIT_FAILURE;
//*** setup the area
icon::Area area; if (!area.load(areaConf, areaName, true)) return EXIT_FAILURE;
icon::ThinSwitcher thinResolver(gridPath);
if (!area.setCoords(&(thinResolver.getClon()), &(thinResolver.getClat()))) return EXIT_FAILURE;
//*** fill the forecast
const icon::Resolver pathResolver(inBaseDir);
icon::AreaFcst areaFcst({ &area }, analysis, pathResolver, thinResolver);
areaFcst.set(params, hours);
icon::area::Forecast forecast=areaFcst().front();
//*** if we asked for text out of all points, print them in text files, one per hour and param
if (outDir.length())
{
for (int hr : hours)
for (const string ¶m : params)
{
string path=outDir+"/"+_txtFileName(param, hr);
if (!_writeHourlyParam(path, forecast[hr][param], (float*)(thinResolver.getClon().bytes),
(float*)(thinResolver.getClat().bytes))) return EXIT_FAILURE;
}
return EXIT_SUCCESS;
}
//*** else, reduce all-area-points forecast to a single-point meteogram
pair<float,float> loLa=area.getRepresentPoint();
icon::area::FcstToMeteogram converter(area, loLa.first, loLa.second, 0, analysis);
// decide how each parameter values over all-area-points get reduced to single one
icon::area::Reducers reducers=icon::area::getDefaultReducers(); // defaults to avegage-reduction
// use max-vector reducer for 10m wind
icon::area::ReduceDuoMax w10m;
reducers[u_10m]=reducers[v_10m]=(icon::area::IReduce*)&w10m;
// use max-diff reducer for all accumulated fields
icon::area::ReduceCumuMax cumuReducers[sizeof(icon::cummulativeFields)];
int i=0;
for (const char *cumuPar : icon::cummulativeFields)
{ cumuReducers[i]=icon::area::ReduceCumuMax(-1); reducers[cumuPar]=&cumuReducers[i]; i++; }
// convert to a meteogram
fprintf(stderr, "converting %s forecast to meteogram\n", area.name().c_str());
icon::Meteogram meteogram=converter(forecast, reducers);
fprintf(stdout, "%s\n", string(meteogram).c_str());
return EXIT_SUCCESS;
}