Merge branch 'feature/MIR-673' into develop

etc/mir/grib-input.yaml

@@ -19,6 +19,10 @@ md5GridSection=f5dc74ec36353f4c83f7de3bf46e1aef, centre:i=98, gridName=O640, edi
 md5GridSection=6ae0076681a2541abe053ab52fa60f30, centre:i=84/85, gridType=lambert, edition=1:
   - first_point_bottom_left: true
 
+# gridType=polar_stereographic does not support the default area method
+gridType=polar_stereographic:
+  - area-mode: mask
+
 # gridType=space_view does not support the default interpolation method
 gridType=space_view:
   - interpolation: k-nearest

src/mir/repres/regular/PolarStereographic.cc

@@ -12,7 +12,14 @@
 
 #include "mir/repres/regular/PolarStereographic.h"
 
+#include <cmath>
+#include <sstream>
+
+#include "mir/api/MIRJob.h"
+#include "mir/param/MIRParametrisation.h"
+#include "mir/util/Angles.h"
 #include "mir/util/Exceptions.h"
+#include "mir/util/Grib.h"
 
 
 namespace mir::repres::regular {
@@ -22,11 +29,67 @@ static const RepresentationBuilder<PolarStereographic> __builder("polar_stereogr
 
 
 PolarStereographic::PolarStereographic(const param::MIRParametrisation& param) :
-    RegularGrid(param, make_proj_spec(param)) {}
+    RegularGrid(param, make_proj_spec(param)),
+    LaDInDegrees_(0),
+    orientationOfTheGridInDegrees_(0),
+    southPoleOnProjectionPlane_(false) {
+    ASSERT(param.get("proj", proj_) && !proj_.empty());
+
+    param.get("LaDInDegrees", LaDInDegrees_);
+    param.get("orientationOfTheGridInDegrees", orientationOfTheGridInDegrees_);
+
+    long south = 0;
+    param.get("southPoleOnProjectionPlane", south);
+    southPoleOnProjectionPlane_ = south != 0;
+}
+
+
+void PolarStereographic::fillGrib(grib_info& info) const {
+    info.grid.grid_type = CODES_UTIL_GRID_SPEC_POLAR_STEREOGRAPHIC;
+    info.extra_set("edition", 2L);  // write LaDInDegrees
+
+    Point2 first   = {x().front(), y().front()};
+    Point2 firstLL = grid().projection().lonlat(first);
+
+    info.grid.latitudeOfFirstGridPointInDegrees  = firstLL[LLCOORDS::LAT];
+    info.grid.longitudeOfFirstGridPointInDegrees = util::normalise_longitude(firstLL[LLCOORDS::LON], 0);
+
+    info.grid.Ni = static_cast<long>(x().size());
+    info.grid.Nj = static_cast<long>(y().size());
+
+    info.grid.iScansNegatively = x().front() < x().back() ? 0L : 1L;
+    info.grid.jScansPositively = y().front() < y().back() ? 1L : 0L;
+
+    info.extra_set("DxInMetres", std::abs(x().step()));
+    info.extra_set("DyInMetres", std::abs(y().step()));
+    info.extra_set("LaDInDegrees", LaDInDegrees_);
+    info.extra_set("orientationOfTheGridInDegrees", util::normalise_longitude(orientationOfTheGridInDegrees_, 0));
+
+    // some extra keys are edition-specific, so parent call is here
+    RegularGrid::fillGrib(info);
+}
+
+
+void PolarStereographic::fillJob(api::MIRJob& job) const {
+    Point2 first   = {x().front(), y().front()};
+    Point2 firstLL = grid().projection().lonlat(first);
+
+    std::ostringstream grid;
+    grid << "gridType=polar_stereographic;"
+         << "proj=" << proj_ << ";"
+         << "LaDInDegrees=" << LaDInDegrees_ << ";"
+         << "orientationOfTheGridInDegrees=" << orientationOfTheGridInDegrees_ << ";"
+         << "southPoleOnProjectionPlane=" << (southPoleOnProjectionPlane_ ? "1" : "0") << ";"
+         << "Ni=" << x()->size() << ";"
+         << "Nj=" << y()->size() << ";"
+         << "grid=" << std::abs(x().step()) << "/" << std::abs(y().step()) << ";"
+         << "latitudeOfFirstGridPointInDegrees=" << firstLL[LLCOORDS::LAT] << ";"
+         << "longitudeOfFirstGridPointInDegrees=" << firstLL[LLCOORDS::LON];
 
+    job.set("grid", grid.str());
 
-void PolarStereographic::fillGrib(grib_info& /*info*/) const {
-    NOTIMP;
+    // some extra keys are edition-specific, so parent call is here
+    RegularGrid::fillJob(job);
 }
 
 

src/mir/repres/regular/PolarStereographic.h

@@ -12,6 +12,8 @@
 
 #pragma once
 
+#include <string>
+
 #include "mir/repres/regular/RegularGrid.h"
 
 
@@ -25,7 +27,8 @@ class PolarStereographic final : public RegularGrid {
 
     // -- Constructors
 
-    explicit PolarStereographic(const param::MIRParametrisation&);
+    PolarStereographic(const param::MIRParametrisation&);
+    PolarStereographic(const PolarStereographic&) = delete;
 
     // -- Destructor
     // None
@@ -34,7 +37,8 @@ class PolarStereographic final : public RegularGrid {
     // None
 
     // -- Operators
-    // None
+
+    PolarStereographic& operator=(const PolarStereographic&) = delete;
 
     // -- Methods
     // None
@@ -50,14 +54,19 @@ class PolarStereographic final : public RegularGrid {
 
 private:
     // -- Members
-    // None
+
+    std::string proj_;
+    double LaDInDegrees_;
+    double orientationOfTheGridInDegrees_;
+    bool southPoleOnProjectionPlane_;
 
     // -- Methods
     // None
 
     // -- Overridden methods
 
     void fillGrib(grib_info&) const override;
+    void fillJob(api::MIRJob&) const override;
 
     // -- Class members
     // None

src/mir/repres/regular/RegularGrid.cc

@@ -23,6 +23,7 @@
 #include "eckit/utils/MD5.h"
 #include "eckit/utils/StringTools.h"
 
+#include "mir/api/MIRJob.h"
 #include "mir/param/MIRParametrisation.h"
 #include "mir/repres/Iterator.h"
 #include "mir/util/Domain.h"
@@ -217,11 +218,26 @@ void RegularGrid::fillGrib(grib_info& info) const {
     shape_.fillGrib(info, grid_.projection().spec());
 
     // scanningMode
-    info.grid.iScansNegatively = x_.back() < x_.front() ? 1L : 0L;
+    info.grid.iScansNegatively = x_.front() < x_.back() ? 0L : 1L;
     info.grid.jScansPositively = y_.front() < y_.back() ? 1L : 0L;
 }
 
 
+void RegularGrid::fillJob(api::MIRJob& job) const {
+    // shape of the reference system
+    shape_.fillJob(job, grid_.projection().spec());
+
+    // scanningMode
+    std::string grid;
+    ASSERT(job.get("grid", grid) && !grid.empty());
+
+    grid += ";iScansNegatively=" + std::to_string(x().front() < x().back() ? 0 : 1);
+    grid += ";jScansPositively=" + std::to_string(y().front() < y().back() ? 1 : 0);
+
+    job.set("grid", grid);
+}
+
+
 bool RegularGrid::includesNorthPole() const {
     return bbox_.north() == Latitude::NORTH_POLE;
 }

src/mir/repres/regular/RegularGrid.h

@@ -102,7 +102,9 @@ class RegularGrid : public Gridded {
     bool isPeriodicWestEast() const override;
 
     void fillGrib(grib_info&) const override;
+    void fillJob(api::MIRJob&) const override;
     void fillMeshGen(util::MeshGeneratorParameters&) const override;
+
     void reorder(long scanningMode, MIRValuesVector&) const override;
     void validate(const MIRValuesVector&) const override;
     void makeName(std::ostream&) const override;

src/mir/util/Shape.cc

@@ -12,6 +12,9 @@
 
 #include "mir/util/Shape.h"
 
+#include <sstream>
+
+#include "mir/api/MIRJob.h"
 #include "mir/param/MIRParametrisation.h"
 
 
@@ -50,12 +53,6 @@ Shape::Shape(const Projection::Spec& spec) : edition(0) {
 }
 
 
-Shape::Shape(const Shape&) = default;
-
-
-Shape& Shape::operator=(const Shape&) = default;
-
-
 void Shape::fillGrib(grib_info& info, const Projection::Spec& spec) const {
     // GRIB2 encoding of user-provided shape
     if (edition != 2) {
@@ -93,4 +90,48 @@ void Shape::fillGrib(grib_info& info, const Projection::Spec& spec) const {
 }
 
 
+void Shape::fillJob(api::MIRJob& job, const Projection::Spec& spec) const {
+    // GRIB2 encoding of user-provided shape
+    if (edition != 2) {
+        return;
+    }
+
+    std::ostringstream shape;
+
+    // shape given by radius or semi-major/minor axis
+    if (provided) {
+        shape << ";shapeOfTheEarth=" << code;
+        switch (code) {
+            case 1:
+                shape << ";radius", spec.getDouble("radius", a);
+                break;
+            case 3:
+                shape << ";earthMajorAxis=" << spec.getDouble("semi_major_axis", a) / 1000.;
+                shape << ";earthMinorAxis=" << spec.getDouble("semi_minor_axis", b) / 1000.;
+                break;
+            case 7:
+                shape << ";earthMajorAxis=" << spec.getDouble("semi_major_axis", a);
+                shape << ";earthMinorAxis=" << spec.getDouble("semi_minor_axis", b);
+                break;
+            default:
+                break;
+        }
+    }
+    else if (spec.has("radius")) {
+        shape << ";shapeOfTheEarth=" << 1L;
+        shape << ";radius", spec.getDouble("radius");
+    }
+    else if (spec.has("semi_major_axis") && spec.has("semi_minor_axis")) {
+        shape << ";shapeOfTheEarth=" << 7L;
+        shape << ";earthMajorAxis=" << spec.getDouble("semi_major_axis");
+        shape << ";earthMinorAxis=" << spec.getDouble("semi_minor_axis");
+    }
+
+    std::string grid;
+    ASSERT(job.get("grid", grid) && !grid.empty());
+
+    job.set("grid", grid + ";" + shape.str());
+}
+
+
 }  // namespace mir::util

src/mir/util/Shape.h

@@ -16,9 +16,14 @@
 #include "mir/util/Grib.h"
 
 
-namespace mir::param {
+namespace mir {
+namespace param {
 class MIRParametrisation;
-}  // namespace mir::param
+}
+namespace api {
+class MIRJob;
+}
+}  // namespace mir
 
 
 namespace mir::util {
@@ -30,11 +35,17 @@ struct Shape {
     Shape(const param::MIRParametrisation&);
     Shape(const Projection::Spec&);
     Shape() : code(6L), a(0.), b(0.), provided(false) {}
-    Shape(const Shape& other);
+
+    Shape(const Shape& other) = default;
+    Shape(Shape&&)            = default;
+
     virtual ~Shape() = default;
 
-    Shape& operator=(const Shape& other);
-    void fillGrib(grib_info& info, const Projection::Spec&) const;
+    Shape& operator=(const Shape&) = default;
+    Shape& operator=(Shape&&)      = default;
+
+    void fillGrib(grib_info&, const Projection::Spec&) const;
+    void fillJob(api::MIRJob&, const Projection::Spec&) const;
 
     long code;
     double a;

tests/assertions/MIR-673.001.test

@@ -0,0 +1,10 @@
+# input
+gridType=polar_stereographic.grib2
+# mir
+--grid=gridType=polar_stereographic
+# grib_get
+gridType=polar_stereographic
+Ni=2869
+Nj=2869
+Dx=2500000
+Dy=2500000

tests/assertions/MIR-673.002.test

@@ -0,0 +1,10 @@
+# input
+gridType=polar_stereographic.grib2
+# mir
+--area=66.5/-24.5/63.4/-13.5
+# grib_get
+gridType=polar_stereographic
+Ni=2869
+Nj=2869
+Dx=2500000
+Dy=2500000

tests/assertions/MIR-673.003.fail.test

@@ -0,0 +1,4 @@
+# input
+gridType=polar_stereographic.grib2
+# mir
+--area=66.5/-24.5/63.4/-13.5 --area-mode=crop

tests/assertions/gridType=polar_stereographic.grib2

@@ -0,0 +1 @@
+../data/gridType=polar_stereographic.grib2

tests/plans/MIR-673.001.test

@@ -0,0 +1,6 @@
+# input
+gridType=polar_stereographic.grib2
+# mir
+--grid=gridType=polar_stereographic
+# plan
+Copy[]

tests/plans/MIR-673.002.test

@@ -0,0 +1,6 @@
+# input
+gridType=polar_stereographic.grib2
+# mir
+--area=66.5/-24.5/63.4/-13.5
+# plan
+AreaMasker[bbox=BoundingBox[north=66.5,west=-24.5,south=63.4,east=-13.5]]|Save[output=...]

tests/plans/MIR-673.003.test

@@ -0,0 +1,6 @@
+# input
+gridType=polar_stereographic.grib2
+# mir
+--area=66.5/-24.5/63.4/-13.5 --area-mode=crop
+# plan
+AreaCropper[bbox=BoundingBox[north=66.5,west=-24.5,south=63.4,east=-13.5]]|Save[output=...]

tests/plans/MIR-673.004.test

@@ -0,0 +1,6 @@
+# input
+gridType=polar_stereographic.grib2
+# mir
+--same=data.in.MIR-673.003.test
+# plan
+Gridded2TypedGrid[grid={"LaDInDegrees":"90","Ni":"2869","Nj":"2869","grid":"2500/2500","gridType":"polar_stereographic","iScansNegatively":"0","jScansPositively":"1","latitudeOfFirstGridPointInDegrees":"42.1138","longitudeOfFirstGridPointInDegrees":"-71.6168","orientationOfTheGridInDegrees":"-30","proj":"+proj=stere +lat_ts=90.000003 +lat_0=90 +lon_0=-30.000001 +k_0=1 +x_0=0 +y_0=0 +R=6371229.000000","shapeOfTheEarth":"6","southPoleOnProjectionPlane":"0"},interpolation=linear,method=FiniteElement[name=linear,nonLinear[MissingIfHeaviestMissing[]],solver=Multiply[],cropping=none,lsmWeightAdjustment=0.2,pruneEpsilon=1e-10,poleDisplacement=0,validateMesh=0,projectionFail=missing-value]]|Save[output=...]

tests/plans/gridType=polar_stereographic.grib2

@@ -0,0 +1 @@
+../data/gridType=polar_stereographic.grib2