projections/proj/nzmg.hpp

0001 // Boost.Geometry - gis-projections (based on PROJ4)
0002
0003 // Copyright (c) 2008-2015 Barend Gehrels, Amsterdam, the Netherlands.
0004
0005 // This file was modified by Oracle on 2017, 2018, 2019.
0006 // Modifications copyright (c) 2017-2019, Oracle and/or its affiliates.
0007 // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle.
0008
0009 // Use, modification and distribution is subject to the Boost Software License,
0010 // Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
0011 // http://www.boost.org/LICENSE_1_0.txt)
0012
0013 // This file is converted from PROJ4, http://trac.osgeo.org/proj
0014 // PROJ4 is originally written by Gerald Evenden (then of the USGS)
0015 // PROJ4 is maintained by Frank Warmerdam
0016 // PROJ4 is converted to Boost.Geometry by Barend Gehrels
0017
0018 // Last updated version of proj: 5.0.0
0019
0020 // Original copyright notice:
0021
0022 // Purpose:  Implementation of the nzmg (New Zealand Map Grid) projection.
0023 //           Very loosely based upon DMA code by Bradford W. Drew
0024 // Author:   Gerald Evenden
0025 // Copyright (c) 1995, Gerald Evenden
0026
0027 // Permission is hereby granted, free of charge, to any person obtaining a
0028 // copy of this software and associated documentation files (the "Software"),
0029 // to deal in the Software without restriction, including without limitation
0030 // the rights to use, copy, modify, merge, publish, distribute, sublicense,
0031 // and/or sell copies of the Software, and to permit persons to whom the
0032 // Software is furnished to do so, subject to the following conditions:
0033
0034 // The above copyright notice and this permission notice shall be included
0035 // in all copies or substantial portions of the Software.
0036
0037 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
0038 // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
0039 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
0040 // THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
0041 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
0042 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
0043 // DEALINGS IN THE SOFTWARE.
0044
0045 #ifndef BOOST_GEOMETRY_PROJECTIONS_NZMG_HPP
0046 #define BOOST_GEOMETRY_PROJECTIONS_NZMG_HPP
0047
0048 #include <boost/geometry/util/math.hpp>
0049
0050 #include <boost/geometry/srs/projections/impl/base_static.hpp>
0051 #include <boost/geometry/srs/projections/impl/base_dynamic.hpp>
0052 #include <boost/geometry/srs/projections/impl/projects.hpp>
0053 #include <boost/geometry/srs/projections/impl/factory_entry.hpp>
0054 #include <boost/geometry/srs/projections/impl/pj_zpoly1.hpp>
0055
0056 namespace boost { namespace geometry
0057 {
0058
0059 namespace projections
0060 {
0061     #ifndef DOXYGEN_NO_DETAIL
0062     namespace detail { namespace nzmg
0063     {
0064
0065             static const double epsilon = 1e-10;
0066             static const int Nbf = 5;
0067             static const int Ntpsi = 9;
0068             static const int Ntphi = 8;
0069
0070             template <typename T>
0071             inline T sec5_to_rad() { return 0.4848136811095359935899141023; }
0072             template <typename T>
0073             inline T rad_to_sec5() { return 2.062648062470963551564733573; }
0074
0075             template <typename T>
0076             inline const pj_complex<T> * bf()
0077             {
0078                 static const pj_complex<T> result[] = {
0079                     {.7557853228,    0.0},
0080                     {.249204646,    .003371507},
0081                     {-.001541739,    .041058560},
0082                     {-.10162907,    .01727609},
0083                     {-.26623489,    -.36249218},
0084                     {-.6870983,    -1.1651967}
0085                 };
0086                 return result;
0087             }
0088
0089             template <typename T>
0090             inline const T * tphi()
0091             {
0092                 static const T result[] = { 1.5627014243, .5185406398, -.03333098,
0093                                             -.1052906,   -.0368594,     .007317,
0094                                              .01220,      .00394,      -.0013 };
0095                 return result;
0096             }
0097             template <typename T>
0098             inline const T * tpsi()
0099             {
0100                 static const T result[] = { .6399175073, -.1358797613, .063294409, -.02526853, .0117879,
0101                                            -.0055161,     .0026906,   -.001333,     .00067,   -.00034 };
0102                 return result;
0103             }
0104
0105             template <typename T, typename Parameters>
0106             struct base_nzmg_ellipsoid
0107             {
0108                 // FORWARD(e_forward)  ellipsoid
0109                 // Project coordinates from geographic (lon, lat) to cartesian (x, y)
0110                 inline void fwd(Parameters const& par, T const& lp_lon, T lp_lat, T& xy_x, T& xy_y) const
0111                 {
0112                     static const T rad_to_sec5 = nzmg::rad_to_sec5<T>();
0113
0114                     pj_complex<T> p;
0115                     const T * C;
0116                     int i;
0117
0118                     lp_lat = (lp_lat - par.phi0) * rad_to_sec5;
0119                     for (p.r = *(C = tpsi<T>() + (i = Ntpsi)); i ; --i)
0120                         p.r = *--C + lp_lat * p.r;
0121                     p.r *= lp_lat;
0122                     p.i = lp_lon;
0123                     p = pj_zpoly1(p, bf<T>(), Nbf);
0124                     xy_x = p.i;
0125                     xy_y = p.r;
0126                 }
0127
0128                 // INVERSE(e_inverse)  ellipsoid
0129                 // Project coordinates from cartesian (x, y) to geographic (lon, lat)
0130                 inline void inv(Parameters const& par, T const& xy_x, T const& xy_y, T& lp_lon, T& lp_lat) const
0131                 {
0132                     static const T sec5_to_rad = nzmg::sec5_to_rad<T>();
0133
0134                     int nn, i;
0135                     pj_complex<T> p, f, fp, dp;
0136                     T den;
0137                     const T* C;
0138
0139                     p.r = xy_y;
0140                     p.i = xy_x;
0141                     for (nn = 20; nn ;--nn) {
0142                         f = pj_zpolyd1(p, bf<T>(), Nbf, &fp);
0143                         f.r -= xy_y;
0144                         f.i -= xy_x;
0145                         den = fp.r * fp.r + fp.i * fp.i;
0146                         p.r += dp.r = -(f.r * fp.r + f.i * fp.i) / den;
0147                         p.i += dp.i = -(f.i * fp.r - f.r * fp.i) / den;
0148                         if ((fabs(dp.r) + fabs(dp.i)) <= epsilon)
0149                             break;
0150                     }
0151                     if (nn) {
0152                         lp_lon = p.i;
0153                         for (lp_lat = *(C = tphi<T>() + (i = Ntphi)); i ; --i)
0154                             lp_lat = *--C + p.r * lp_lat;
0155                         lp_lat = par.phi0 + p.r * lp_lat * sec5_to_rad;
0156                     } else
0157                         lp_lon = lp_lat = HUGE_VAL;
0158                 }
0159
0160                 static inline std::string get_name()
0161                 {
0162                     return "nzmg_ellipsoid";
0163                 }
0164
0165             };
0166
0167             // New Zealand Map Grid
0168             template <typename Parameters>
0169             inline void setup_nzmg(Parameters& par)
0170             {
0171                 typedef typename Parameters::type calc_t;
0172                 static const calc_t d2r = geometry::math::d2r<calc_t>();
0173
0174                 /* force to International major axis */
0175                 par.a = 6378388.0;
0176                 par.ra = 1. / par.a;
0177                 par.lam0 = 173. * d2r;
0178                 par.phi0 = -41. * d2r;
0179                 par.x0 = 2510000.;
0180                 par.y0 = 6023150.;
0181             }
0182
0183     }} // namespace detail::nzmg
0184     #endif // doxygen
0185
0186     /*!
0187         \brief New Zealand Map Grid projection
0188         \ingroup projections
0189         \tparam Geographic latlong point type
0190         \tparam Cartesian xy point type
0191         \tparam Parameters parameter type
0192         \par Projection characteristics
0193          - Fixed Earth
0194         \par Example
0195         \image html ex_nzmg.gif
0196     */
0197     template <typename T, typename Parameters>
0198     struct nzmg_ellipsoid : public detail::nzmg::base_nzmg_ellipsoid<T, Parameters>
0199     {
0200         template <typename Params>
0201         inline nzmg_ellipsoid(Params const& , Parameters & par)
0202         {
0203             detail::nzmg::setup_nzmg(par);
0204         }
0205     };
0206
0207     #ifndef DOXYGEN_NO_DETAIL
0208     namespace detail
0209     {
0210
0211         // Static projection
0212         BOOST_GEOMETRY_PROJECTIONS_DETAIL_STATIC_PROJECTION_FI(srs::spar::proj_nzmg, nzmg_ellipsoid)
0213
0214         // Factory entry(s)
0215         BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_ENTRY_FI(nzmg_entry, nzmg_ellipsoid)
0216
0217         BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_BEGIN(nzmg_init)
0218         {
0219             BOOST_GEOMETRY_PROJECTIONS_DETAIL_FACTORY_INIT_ENTRY(nzmg, nzmg_entry)
0220         }
0221
0222     } // namespace detail
0223     #endif // doxygen
0224
0225 } // namespace projections
0226
0227 }} // namespace boost::geometry
0228
0229 #endif // BOOST_GEOMETRY_PROJECTIONS_NZMG_HPP
0230