From 7a2da20c5d099da3eb93254e7bb05b08636e7d3d Mon Sep 17 00:00:00 2001 From: Pedro Maciel Date: Sun, 8 Oct 2023 08:23:12 +0100 Subject: [PATCH] eckit::geo::projection::Mercator --- src/eckit/geo/CMakeLists.txt | 4 + src/eckit/geo/grid/regular/Mercator.cc | 59 ++++++++++++++ src/eckit/geo/grid/regular/Mercator.h | 89 +++++++++++++++++++++ src/eckit/geo/projection/Mercator.cc | 106 +++++++++++++++++++++++++ src/eckit/geo/projection/Mercator.h | 99 +++++++++++++++++++++++ 5 files changed, 357 insertions(+) create mode 100644 src/eckit/geo/grid/regular/Mercator.cc create mode 100644 src/eckit/geo/grid/regular/Mercator.h create mode 100644 src/eckit/geo/projection/Mercator.cc create mode 100644 src/eckit/geo/projection/Mercator.h diff --git a/src/eckit/geo/CMakeLists.txt b/src/eckit/geo/CMakeLists.txt index 57d76ef2a..3fd2b2d15 100644 --- a/src/eckit/geo/CMakeLists.txt +++ b/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 diff --git a/src/eckit/geo/grid/regular/Mercator.cc b/src/eckit/geo/grid/regular/Mercator.cc new file mode 100644 index 000000000..1f7f3be80 --- /dev/null +++ b/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 __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& Mercator::longitudes() const { + NOTIMP; +} + + +const std::vector& Mercator::latitudes() const { + NOTIMP; +} + + +} // namespace eckit::geo::grid::regular diff --git a/src/eckit/geo/grid/regular/Mercator.h b/src/eckit/geo/grid/regular/Mercator.h new file mode 100644 index 000000000..42ea855de --- /dev/null +++ b/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& longitudes() const override; + const std::vector& latitudes() const override; + + // -- Class members + // None + + // -- Class methods + // None + + // -- Friends + // None +}; + + +} // namespace eckit::geo::grid::regular diff --git a/src/eckit/geo/projection/Mercator.cc b/src/eckit/geo/projection/Mercator.cc new file mode 100644 index 000000000..8236b25d4 --- /dev/null +++ b/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 + +#include "eckit/config/MappedConfiguration.h" +#include "eckit/geo/util.h" +#include "eckit/types/FloatCompare.h" + + +namespace eckit::geo::projection { + + +static ProjectionBuilder __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 diff --git a/src/eckit/geo/projection/Mercator.h b/src/eckit/geo/projection/Mercator.h new file mode 100644 index 000000000..d2479ab09 --- /dev/null +++ b/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 + +#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_; + 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(p)); } + Point inv(const Point& q) const override { return inv(std::get(q)); } + + // -- Class members + // None + + // -- Class methods + // None + + // -- Friends + // None +}; + + +} // namespace eckit::geo::projection