Tile

src/mir/CMakeLists.txt

@@ -465,6 +465,8 @@ list(APPEND mir_srcs
     output/RawOutput.h
     output/ResizableOutput.cc
     output/ResizableOutput.h
+    output/Tile.cc
+    output/Tile.h
     output/ValuesOutput.cc
     output/ValuesOutput.h
     output/VectorOutput.cc

src/mir/output/Tile.cc

@@ -0,0 +1,261 @@
+/*
+ * (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 "mir/output/Tile.h"
+
+#include <fstream>
+#include <iomanip>
+#include <iostream>
+#include <memory>
+#include <ostream>
+#include <string>
+#include <vector>
+
+#include "eckit/filesystem/PathName.h"
+#include "mir/action/area/AreaCropper.h"
+#include "mir/action/context/Context.h"
+#include "mir/data/MIRField.h"
+#include "mir/output/PNGOutput.h"
+#include "mir/param/CombinedParametrisation.h"
+#include "mir/param/DefaultParametrisation.h"
+#include "mir/param/MIRParametrisation.h"
+#include "mir/param/RuntimeParametrisation.h"
+#include "mir/repres/Representation.h"
+#include "mir/util/BoundingBox.h"
+#include "mir/util/Domain.h"
+#include "mir/util/Exceptions.h"
+#include "mir/util/Log.h"
+
+
+namespace mir::output {
+
+
+namespace {
+
+
+struct Tiling {
+    explicit Tiling(const param::MIRParametrisation& param) {
+        std::vector<long> tile;
+        ASSERT(param.get("tile", tile));
+
+        if (tile.size() == 4) {
+            single = true;
+
+            ASSERT(1 <= tile[0]);
+            ASSERT(1 <= tile[1]);
+            Ti = static_cast<size_t>(tile[0]);
+            Tj = static_cast<size_t>(tile[1]);
+
+            ASSERT(0 <= tile[2] && tile[2] < Ti);
+            ASSERT(0 <= tile[3] && tile[3] < Tj);
+            i = static_cast<size_t>(tile[2]);
+            j = static_cast<size_t>(tile[3]);
+
+            return;
+        }
+
+        if (tile.size() == 2) {
+            single = false;
+
+            ASSERT(1 <= tile[0]);
+            ASSERT(1 <= tile[1]);
+            Ti = static_cast<size_t>(tile[0]);
+            Tj = static_cast<size_t>(tile[1]);
+            i  = 0;
+            j  = 0;
+
+            return;
+        }
+
+        throw exception::BadValue("Tile: invalid tiling", Here());
+    }
+
+    friend bool operator==(const Tiling& a, const Tiling& b) {
+        return a.Ti == b.Ti && a.Tj == b.Tj && a.i == b.i && a.j == b.j && a.single == b.single;
+    }
+
+    friend std::ostream& operator<<(std::ostream& out, const Tiling& tiling) {
+        if (tiling.single) {
+            out << "Tile[Ni=" << tiling.Ti << ",Nj=" << tiling.Tj << ",i=" << tiling.i << ",j=" << tiling.j << "]";
+        }
+        else {
+            out << "Tile[Ni=" << tiling.Ti << ",Nj=" << tiling.Tj << "]";
+        }
+        return out;
+    }
+
+    struct tile_type {
+        tile_type(size_t _i, size_t _j, const util::BoundingBox& _bbox) : i(_i), j(_j), bbox(_bbox) {}
+
+        std::vector<double> make_area() const {
+            return {bbox.north().value(), bbox.west().value(), bbox.south().value(), bbox.east().value()};
+        }
+
+        size_t i;
+        size_t j;
+        util::BoundingBox bbox;
+    };
+
+    std::vector<tile_type> make_tiles(const util::BoundingBox& bbox) const {
+        std::vector<tile_type> tiles;
+
+        for (size_t j = 0; j < Tj; ++j) {
+            for (size_t i = 0; i < Ti; ++i) {
+                const auto north = bbox.north().fraction();
+                const auto dlat  = (north - bbox.south().fraction()) / Tj;
+
+                const auto west = bbox.west().fraction();
+                const auto dlon = (bbox.east().fraction() - west) / Ti;
+
+                const Latitude n  = north - j * dlat;
+                const Longitude w = west + i * dlon;
+                const Latitude s  = n - dlat;
+                const Longitude e = w + dlon;
+
+                tiles.emplace_back(i, j, util::BoundingBox{n, w, s, e});
+            }
+        }
+
+        return tiles;
+    }
+
+    size_t Ti;
+    size_t Tj;
+    size_t i;
+    size_t j;
+    bool single;
+};
+
+
+std::string get_path(const std::string& path, size_t i, size_t j) {
+    const eckit::PathName p(path);
+
+    std::ostringstream s;
+    s << p.dirName() / p.baseName(false);
+    s << "." << std::setw(4) << std::setfill('0') << i;
+    s << "." << std::setw(4) << std::setfill('0') << j;
+    s << p.extension();
+
+    return s.str();
+};
+
+
+}  // namespace
+
+
+size_t Tile::save(const param::MIRParametrisation& param, context::Context& ctx) {
+    Tiling tiling(param);
+
+    Log::info() << "Tile::save " << tiling << std::endl;
+
+    auto& field = ctx.field();
+    ASSERT(field.dimensions() == 1);
+
+    repres::RepresentationHandle repres(field.representation());
+    auto tiles = tiling.make_tiles(repres->domain());
+    ASSERT(tiles.size() == tiling.Ti * tiling.Tj);
+
+    if (tiling.single) {
+        decltype(tiles) one{tiles[tiling.j * tiling.Ti + tiling.i]};
+        ASSERT(one[0].i == tiling.i);
+        ASSERT(one[0].j == tiling.j);
+
+        one.swap(tiles);
+    }
+
+    std::string kml;
+    param.get("tile-kml", kml);
+
+    for (const auto& tile : tiles) {
+        // set context
+        param::RuntimeParametrisation runtime_user(param.userParametrisation());
+
+        const std::vector<double> area(tile.make_area());
+        runtime_user.set("area", area);
+
+        runtime_user.unset("format");
+        if (std::string format; param.get("tile-format", format)) {
+            runtime_user.set("format", format);
+        }
+
+        static const param::DefaultParametrisation defaults;
+        param::CombinedParametrisation runtime(runtime_user, param.fieldParametrisation(), defaults);
+
+        // crop
+        context::Context local(ctx);
+        std::unique_ptr<action::Action>(new action::AreaCropper(runtime))->perform(local);
+
+        // save
+        auto path = get_path(path_, tile.j, tile.i);
+        std::unique_ptr<output::MIROutput> output(output::MIROutputFactory::build(path, runtime));
+        output->save(runtime, local);
+    }
+
+    if (!kml.empty()) {
+        std::ofstream file(kml);
+
+        file << "<?xml version=\"1.0\" encoding=\"UTF-8\"?>\n"
+             << "<kml xmlns=\"http://www.opengis.net/kml/2.2\">\n"
+             << "<Document>\n"
+             << "  <name>Image Grid Overlay</name>\n";
+
+        for (const auto& tile : tiles) {
+            const auto img = get_path(path_, tile.j, tile.i);
+            file << "  <GroundOverlay>\n"
+                 << "    <name>" << img << "</name>\n"
+                 << "    <Icon>\n"
+                 << "      <href>" << img << "</href>\n"
+                 << "    </Icon>\n"
+                 << "    <LatLonBox>\n"
+                 << "      <north>" << tile.bbox.north().value() << "</north>\n"
+                 << "      <south>" << tile.bbox.south().value() << "</south>\n"
+                 << "      <east>" << tile.bbox.east().value() << "</east>\n"
+                 << "      <west>" << tile.bbox.west().value() << "</west>\n"
+                 << "    </LatLonBox>\n"
+                 << "  </GroundOverlay>\n";
+        }
+
+        // Write KML footer
+        file << "</Document>\n"
+             << "</kml>\n";
+    }
+
+    return tiles.size();
+}
+
+
+bool Tile::sameAs(const MIROutput& other) const {
+    const auto* o = dynamic_cast<const Tile*>(&other);
+    return (o != nullptr);
+}
+
+
+bool Tile::sameParametrisation(const param::MIRParametrisation& param1, const param::MIRParametrisation& param2) const {
+    return Tiling{param1} == Tiling{param2};
+}
+
+
+bool Tile::printParametrisation(std::ostream& out, const param::MIRParametrisation& param) const {
+    out << "tile=" << Tiling{param};
+    return true;
+}
+
+
+void Tile::print(std::ostream& out) const {
+    out << "Tile[]";
+}
+
+
+static const MIROutputBuilder<Tile> output("tile");
+
+
+}  // namespace mir::output

src/mir/output/Tile.h

@@ -0,0 +1,37 @@
+/*
+ * (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.
+ */
+
+
+#pragma once
+
+#include "mir/output/MIROutput.h"
+
+
+namespace mir::output {
+
+
+class Tile final : public MIROutput {
+public:
+    explicit Tile(const std::string& path) : path_{path} {}
+
+private:
+    const std::string path_;
+
+    // From MIROutput
+    size_t save(const param::MIRParametrisation&, context::Context&) override;
+    bool sameAs(const MIROutput&) const override;
+    bool sameParametrisation(const param::MIRParametrisation&, const param::MIRParametrisation&) const override;
+    bool printParametrisation(std::ostream&, const param::MIRParametrisation&) const override;
+    void print(std::ostream&) const override;
+};
+
+
+}  // namespace mir::output

src/tools/mir.cc

@@ -240,7 +240,7 @@ struct MIR : MIRTool {
             new VectorOption<double>("mode-real-values", "Statistics mode bin ranges value (default 0/1)", 0));
 
         //==============================================
-        options_.push_back(new Separator("Filtering"));
+        options_.push_back(new Separator("Filtering (gridded)"));
         options_.push_back(new FactoryOption<key::Area>("area", "area cropping/masking"));
         options_.push_back(new FactoryOption<key::Area::Mode>("area-mode", "area cropping/masking mode"));
         options_.push_back(new SimpleOption<eckit::PathName>("bitmap", "Bitmap file to apply"));
@@ -262,6 +262,8 @@ struct MIR : MIRTool {
             "Unstructured grid globalise minimum distance to insert missing values if needed (default 555975. [m])"));
         options_.push_back(new SimpleOption<bool>("unstructured", "Convert to unstructured grid"));
 
+        //==============================================
+        options_.push_back(new Separator("Filtering (spectral)"));
         options_.push_back(new SimpleOption<bool>("cesaro", "Cesàro summation filtering"));
         options_.push_back(new SimpleOption<double>("cesaro-k", "Cesàro summation k (default 2.)"));
         options_.push_back(new SimpleOption<size_t>(
@@ -355,6 +357,12 @@ struct MIR : MIRTool {
         options_.push_back(new SimpleOption<eckit::PathName>("plan-script", "File containing a plan definition"));
         options_.push_back(new SimpleOption<eckit::PathName>("dump-weights-info", "Dump weights information to file"));
 
+        //==============================================
+        options_.push_back(new Separator("Tiling"));
+        options_.push_back(new VectorOption<long>("tile", "Tiling Ni/Nj[/i/j]", 0));
+        options_.push_back(new FactoryOption<output::MIROutputFactory>("tile-format", "Output tile format"));
+        options_.push_back(new SimpleOption<std::string>("tile-kml", "Output MKL file"));
+
         //==============================================
         options_.push_back(new Separator("Caching"));
         options_.push_back(new FactoryOption<caching::matrix::MatrixLoaderFactory>(