1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218
|
/*
* (C) Copyright 1996- ECMWF.
*
* This software is licensed under the terms of the Apache Licence Version 2.0
* which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
*
* In applying this licence, ECMWF does not waive the privileges and immunities
* granted to it by virtue of its status as an intergovernmental organisation nor
* does it submit to any jurisdiction.
*/
#include <memory>
#include "eckit/exception/Exceptions.h"
#include "eckit/log/Log.h"
#include "eckit/option/CmdArgs.h"
#include "eckit/option/FactoryOption.h"
#include "eckit/option/SimpleOption.h"
#include "eckit/utils/StringTools.h"
#include "eckit/utils/Tokenizer.h"
#include "mir/action/context/Context.h"
#include "mir/data/MIRField.h"
#include "mir/input/GribFileInput.h"
#include "mir/output/GribFileOutput.h"
#include "mir/param/CombinedParametrisation.h"
#include "mir/param/ConfigurationWrapper.h"
#include "mir/param/DefaultParametrisation.h"
#include "mir/repres/Representation.h"
#include "mir/stats/Method.h"
#include "mir/stats/Statistics.h"
#include "mir/tools/MIRTool.h"
#include "mir/util/MIRStatistics.h"
#include "mir/util/Pretty.h"
using namespace mir;
using prec_t = decltype(std::cout.precision());
class MIRStatistics : public mir::tools::MIRTool {
private:
void execute(const eckit::option::CmdArgs&);
void usage(const std::string& tool) const;
int minimumPositionalArguments() const { return 1; }
struct PerPointStatistics {
static void list(std::ostream& out) { out << eckit::StringTools::join(", ", perPointStats()) << std::endl; }
static std::vector<std::string> perPointStats() { return {"mean", "variance", "stddev"}; }
};
public:
MIRStatistics(int argc, char** argv) : MIRTool(argc, argv) {
using eckit::option::FactoryOption;
using eckit::option::SimpleOption;
options_.push_back(new FactoryOption<stats::StatisticsFactory>(
"statistics", "Statistics methods for interpreting field values"));
options_.push_back(new SimpleOption<double>("counter-lower-limit", "count lower limit"));
options_.push_back(new SimpleOption<double>("counter-upper-limit", "count upper limit"));
options_.push_back(new FactoryOption<PerPointStatistics>(
"output", "/-separated list of per-point statistics (output GRIB to <statistics>"));
options_.push_back(new SimpleOption<size_t>("precision", "Output precision"));
}
};
void MIRStatistics::usage(const std::string& tool) const {
eckit::Log::info() << "\n"
"Calculate field statistics, or per-point statistics when specifying output."
"\n"
"\n"
"Usage: "
<< tool
<< " [--statistics=option] file.grib [file2.grib [...]]"
"\n"
"\n"
"Examples:"
"\n"
" % "
<< tool
<< " file.grib"
"\n"
" % "
<< tool
<< " --statistics=scalar file1.grib file2.grib file3.grib"
"\n"
" % "
<< tool
<< " --statistics=spectral file.grib"
"\n"
" % "
<< tool << " --output=mean/min/max file1.grib file2.grib file3.grib" << std::endl;
}
void MIRStatistics::execute(const eckit::option::CmdArgs& args) {
// Build CombinedParametrisation from a few parts:
// - wrap the arguments, so that they behave as a MIRParametrisation
// - input grib metadata
// - lookup configuration for input metadata-specific parametrisation
const param::ConfigurationWrapper args_wrap(args);
const param::DefaultParametrisation defaults;
auto& log = eckit::Log::info();
prec_t precision;
auto old = args.get("precision", precision) ? log.precision(precision) : log.precision();
// on 'output' option, calculate per-point statistics
std::string output;
if (args_wrap.get("output", output)) {
// read in first field to reset statistics and reference representation
input::GribFileInput firstGribFile(args(0));
ASSERT(firstGribFile.next());
const input::MIRInput& firstInput(firstGribFile);
repres::RepresentationHandle reference(firstInput.field().representation());
size_t Nfirst = firstInput.field().values(0).size();
ASSERT(Nfirst > 0);
log << "Using " << Pretty(Nfirst, {"grid point"}) << std::endl;
// set paramId/metadata-specific method per-point statistics
std::string statistics = "scalar";
args_wrap.get("statistics", statistics);
// per-point statistics
std::unique_ptr<param::MIRParametrisation> param(
new param::CombinedParametrisation(args_wrap, firstGribFile, defaults));
std::unique_ptr<stats::Method> pps(stats::MethodFactory::build(statistics, *param));
pps->resize(Nfirst);
for (auto& arg : args) {
input::GribFileInput grib(arg);
const input::MIRInput& input = grib;
size_t count = 0;
while (grib.next()) {
log << "\n'" << arg << "' #" << ++count << std::endl;
repres::RepresentationHandle repres(input.field().representation());
if (!repres->sameAs(*reference)) {
eckit::Log::error() << "Input not expected,"
"\n"
"expected "
<< *reference
<< "\n"
"but got "
<< *repres;
throw eckit::UserError("input not of the expected format");
}
pps->execute(input.field());
}
}
// Write statistics
eckit::Tokenizer parse("/");
std::vector<std::string> outputs;
parse(output, outputs);
for (auto& j : outputs) {
log << "Writing '" << j << "'" << std::endl;
util::MIRStatistics stats;
context::Context ctx(firstGribFile, stats);
auto& f = ctx.field();
j == "mean" ? pps->mean(f)
: j == "variance" ? pps->variance(f)
: j == "stddev" ? pps->stddev(f)
: throw eckit::UserError("Output " + j + "' not supported");
std::unique_ptr<output::MIROutput> out(new output::GribFileOutput(j));
out->save(*param, ctx);
}
log.precision(old);
return;
}
// per-field statistics
for (auto& arg : args) {
input::GribFileInput grib(arg);
const input::MIRInput& input = grib;
size_t count = 0;
while (grib.next()) {
log << "\n'" << arg << "' #" << ++count << std::endl;
// paramId/metadata-specific method
std::string statistics = "scalar";
args_wrap.get("statistics", statistics);
// Calculate and show statistics
std::unique_ptr<param::MIRParametrisation> param(
new param::CombinedParametrisation(args_wrap, grib, defaults));
std::unique_ptr<stats::Statistics> stats(stats::StatisticsFactory::build(statistics, *param));
stats->execute(input.field());
log << *stats << std::endl;
}
}
log.precision(old);
}
int main(int argc, char** argv) {
MIRStatistics tool(argc, argv);
return tool.start();
}
|