eckit::geo::projection::Mercator

src/eckit/geo/CMakeLists.txt

@@ -71,6 +71,8 @@ list( APPEND eckit_geo_srcs
     grid/reduced/ReducedLL.h
     grid/regular/IrregularLL.cc
     grid/regular/IrregularLL.h
+    grid/regular/Mercator.cc
+    grid/regular/Mercator.h
     grid/regular/RegularGaussian.cc
     grid/regular/RegularGaussian.h
     grid/regular/RegularLL.cc
@@ -89,6 +91,8 @@ list( APPEND eckit_geo_srcs
     polygon/Polygon.h
     projection/LonLatToXYZ.cc
     projection/LonLatToXYZ.h
+    projection/Mercator.cc
+    projection/Mercator.h
     projection/None.cc
     projection/None.h
     projection/Rotation.cc

src/eckit/geo/grid/regular/Mercator.cc

@@ -0,0 +1,59 @@
+/*
+ * (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 "eckit/geo/grid/regular/Mercator.h"
+
+#include "eckit/config/MappedConfiguration.h"
+#include "eckit/exception/Exceptions.h"
+#include "eckit/geo/iterator/Regular.h"
+#include "eckit/geo/util/regex.h"
+#include "eckit/types/Fraction.h"
+#include "eckit/utils/Translator.h"
+
+
+namespace eckit::geo::grid::regular {
+
+
+static const GridRegisterType<Mercator> __grid_type("mercator");
+
+
+Mercator::Mercator(const Configuration& config) :
+    Mercator(config.getUnsigned("ni"), config.getUnsigned("nj"), Increments{config}, area::BoundingBox{config}) {}
+
+
+Mercator::Mercator(size_t ni, size_t nj, const Increments& inc, const area::BoundingBox& bbox) :
+    Regular(bbox), inc_(inc), ni_(ni), nj_(nj) {
+    ASSERT(ni_ > 0);
+    ASSERT(nj_ > 0);
+}
+
+
+Grid::iterator Mercator::cbegin() const {
+    return iterator{new geo::iterator::Regular(*this, 0)};
+}
+
+
+Grid::iterator Mercator::cend() const {
+    return iterator{new geo::iterator::Regular(*this, size())};
+}
+
+const std::vector<double>& Mercator::longitudes() const {
+    NOTIMP;
+}
+
+
+const std::vector<double>& Mercator::latitudes() const {
+    NOTIMP;
+}
+
+
+}  // namespace eckit::geo::grid::regular

src/eckit/geo/grid/regular/Mercator.h

@@ -0,0 +1,89 @@
+/*
+ * (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 "eckit/geo/grid/Regular.h"
+
+#include "eckit/geo/Increments.h"
+
+
+namespace eckit::geo::grid::regular {
+
+
+class Mercator final : public Regular {
+public:
+    // -- Types
+    // None
+
+    // -- Exceptions
+    // None
+
+    // -- Constructors
+
+    explicit Mercator(const Configuration&);
+
+    Mercator(size_t ni, size_t nj, const Increments&, const area::BoundingBox&);
+
+    // -- Destructor
+    // None
+
+    // -- Convertors
+    // None
+
+    // -- Operators
+    // None
+
+    // -- Methods
+    // None
+
+    // -- Overridden methods
+
+    iterator cbegin() const override;
+    iterator cend() const override;
+
+    size_t ni() const override { return ni_; }
+    size_t nj() const override { return nj_; }
+
+    // -- Class members
+    // None
+
+    // -- Class methods
+    // None
+
+private:
+    // -- Members
+
+    Increments inc_;
+    size_t ni_;
+    size_t nj_;
+
+    // -- Methods
+    // None
+
+    // -- Overridden methods
+
+    const std::vector<double>& longitudes() const override;
+    const std::vector<double>& latitudes() const override;
+
+    // -- Class members
+    // None
+
+    // -- Class methods
+    // None
+
+    // -- Friends
+    // None
+};
+
+
+}  // namespace eckit::geo::grid::regular

src/eckit/geo/projection/Mercator.cc

@@ -0,0 +1,106 @@
+/*
+ * (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 "eckit/geo/projection/Mercator.h"
+
+#include <cmath>
+
+#include "eckit/config/MappedConfiguration.h"
+#include "eckit/geo/util.h"
+#include "eckit/types/FloatCompare.h"
+
+
+namespace eckit::geo::projection {
+
+
+static ProjectionBuilder<Mercator> __projection("mercator");
+
+
+Mercator::Mercator(double meridian, double parallel, Figure* figure, PointLonLat first) :
+    meridian_(PointLonLat::normalise_angle_to_minimum(meridian, -180.)),
+    parallel_(parallel),
+    figure_(figure),
+    first_(first),
+    eps_(1e-10),
+    max_iter_(15) {
+    // Map Projections: A Working Manual, John P. Snyder (1987)
+    // - Equation (7-9) to calculate phi iteratively
+    // - Equation (15-11) to calculate t
+    ASSERT(figure_);
+
+    ASSERT_MSG(!types::is_approximately_equal(first.lat, 90.) && !types::is_approximately_equal(first.lat, -90.),
+               "Mercator: projection cannot be calculated at the poles");
+
+    e_    = figure_->eccentricity();
+    lam0_ = util::degree_to_radian * meridian_;
+
+    auto phi0 = util::degree_to_radian * parallel_;
+    auto lam1 = util::degree_to_radian * first.lon;
+    auto phi1 = util::degree_to_radian * first.lat;
+
+    m_ = figure_->a() * std::cos(phi0) / (std::sqrt(1. - e_ * e_ * std::sin(phi0) * std::sin(phi0)));
+    ASSERT(!types::is_approximately_equal(m_, 0.));
+
+    w_  = 1. / m_;
+    x0_ = m_ * (lam0_ - lam1);
+    y0_ = m_ * std::log(std::tan(M_PI_4 - 0.5 * phi1) /
+                        std::pow(((1. - e_ * std::sin(phi1)) / (1. + e_ * std::sin(phi1))), 0.5 * e_));
+
+    ASSERT(types::is_approximately_equal(phi1, calculate_phi(std::exp(y0_ * w_)), eps_));
+}
+
+
+Mercator::Mercator(const Configuration& config) :
+    Mercator(config.getDouble("meridian"),
+             config.getDouble("parallel"),
+             FigureFactory::instance().get(config.getString("figure")).create(config),
+             PointLonLat{config.getDouble("lon0"), config.getDouble("lat0")}) {}
+
+
+double Mercator::calculate_phi(double t) const {
+    auto phi = M_PI_2 - 2 * std::atan(t);
+    for (size_t i = 0; i <= max_iter_; i++) {
+        auto es   = e_ * std::sin(phi);
+        auto dphi = M_PI_2 - 2 * std::atan(t * (std::pow(((1. - es * es) / (1. + es * es)), 0.5 * e_))) - phi;
+
+        phi += dphi;
+        if (types::is_approximately_equal(dphi, 0., eps_)) {
+            return phi;
+        }
+    }
+
+    throw SeriousBug("Mercator: phi calculation failed to converge");
+}
+
+
+Point2 Mercator::fwd(const PointLonLat& p) const {
+    auto phi = util::degree_to_radian * p.lat;
+    auto lam = util::degree_to_radian * p.lon;
+
+    return {x0_ + m_ * (lam - lam0_),
+            y0_ - m_ * std::log(std::tan(M_PI_4 - 0.5 * phi) /
+                                std::pow(((1. - e_ * std::sin(phi)) / (1. + e_ * std::sin(phi))), 0.5 * e_))};
+}
+
+
+PointLonLat Mercator::inv(const Point2& q) const {
+    return PointLonLat::make(util::radian_to_degree * (lam0_ + (q.X - x0_) * w_),
+                             util::radian_to_degree * calculate_phi(std::exp(-(q.Y - y0_) * w_)));
+}
+
+
+Projection::Spec Mercator::spec() const {
+    return Spec{{{"lad", parallel_}, {"orientation", meridian_}}};
+}
+
+
+}  // namespace eckit::geo::projection

src/eckit/geo/projection/Mercator.h

@@ -0,0 +1,99 @@
+/*
+ * (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 <memory>
+
+#include "eckit/geo/Figure.h"
+#include "eckit/geo/Projection.h"
+
+
+namespace eckit::geo::projection {
+
+
+/// Calculate coordinates of a point on a rotated sphere given new location of South Pole (vector) and angle
+class Mercator final : public Projection {
+public:
+    // -- Exceptions
+    // None
+
+    // -- Constructors
+
+    Mercator(double meridian, double parallel, Figure*, PointLonLat first = {0, 0});
+
+    explicit Mercator(const Configuration&);
+
+    // -- Destructor
+    // None
+
+    // -- Convertors
+    // None
+
+    // -- Operators
+    // None
+
+    // -- Methods
+
+    Point2 fwd(const PointLonLat& p) const;
+    PointLonLat inv(const Point2& q) const;
+
+    // -- Overridden methods
+
+    Spec spec() const override;
+
+    // -- Class members
+    // None
+
+    // -- Class methods
+    // None
+
+private:
+    // -- Members
+
+    const double meridian_;  // angle [degree] between Eastward direction and the Equator, range [0, 90]
+    const double parallel_;  // latitude [degree] of projection intersecting ellipsoid
+
+    const std::unique_ptr<Figure> figure_;
+    PointLonLat first_;
+
+    const double eps_;
+    const size_t max_iter_;
+
+    double lam0_;
+    double x0_;
+    double y0_;
+    double e_;
+    double m_;
+    double w_;
+
+    // -- Methods
+
+    double calculate_phi(double t) const;
+
+    // -- Overridden methods
+
+    Point fwd(const Point& p) const override { return fwd(std::get<PointLonLat>(p)); }
+    Point inv(const Point& q) const override { return inv(std::get<Point2>(q)); }
+
+    // -- Class members
+    // None
+
+    // -- Class methods
+    // None
+
+    // -- Friends
+    // None
+};
+
+
+}  // namespace eckit::geo::projection