File: mir-statistics.cc

package info (click to toggle)
metview 5.10.2-1
  • links: PTS, VCS
  • area: main
  • in suites: bullseye
  • size: 242,296 kB
  • sloc: cpp: 437,117; ansic: 41,433; xml: 19,944; f90: 13,059; sh: 6,562; python: 3,953; yacc: 1,774; lex: 1,121; perl: 701; makefile: 92
file content (218 lines) | stat: -rw-r--r-- 7,950 bytes parent folder | download
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();
}