← Back to branch summary

~ubuntu-branches/ubuntu/saucy/merkaartor/saucy

« back to all changes in this revision

Viewing changes to include/ggl/extensions/gis/projections/proj/tmerc.hpp

Tags: upstream-0.15.3+svn20934
Import upstream version 0.15.3+svn20934

Show diffs side-by-side

added added

removed removed

Lines of Context:
include/ggl/extensions/gis/projections/proj/tmerc.hpp
1
 
#ifndef GGL_PROJECTIONS_TMERC_HPP
2
 
#define GGL_PROJECTIONS_TMERC_HPP
3
 
 
4
 
// Generic Geometry Library - projections (based on PROJ4)
5
 
// This file is automatically generated. DO NOT EDIT.
6
 
 
7
 
// Copyright Barend Gehrels (1995-2009), Geodan Holding B.V. Amsterdam, the Netherlands.
8
 
// Copyright Bruno Lalande (2008-2009)
9
 
// Use, modification and distribution is subject to the Boost Software License,
10
 
// Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
11
 
// http://www.boost.org/LICENSE_1_0.txt)
12
 
 
13
 
// This file is converted from PROJ4, http://trac.osgeo.org/proj
14
 
// PROJ4 is originally written by Gerald Evenden (then of the USGS)
15
 
// PROJ4 is maintained by Frank Warmerdam
16
 
// PROJ4 is converted to Geometry Library by Barend Gehrels (Geodan, Amsterdam)
17
 
 
18
 
// Original copyright notice:
19
 
 
20
 
// Permission is hereby granted, free of charge, to any person obtaining a
21
 
// copy of this software and associated documentation files (the "Software"),
22
 
// to deal in the Software without restriction, including without limitation
23
 
// the rights to use, copy, modify, merge, publish, distribute, sublicense,
24
 
// and/or sell copies of the Software, and to permit persons to whom the
25
 
// Software is furnished to do so, subject to the following conditions:
26
 
 
27
 
// The above copyright notice and this permission notice shall be included
28
 
// in all copies or substantial portions of the Software.
29
 
 
30
 
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
31
 
// OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
32
 
// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
33
 
// THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
34
 
// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
35
 
// FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
36
 
// DEALINGS IN THE SOFTWARE.
37
 
 
38
 
#include <boost/concept_check.hpp>
39
 
#include <boost/math/special_functions/hypot.hpp>
40
 
 
41
 
#include <ggl/extensions/gis/projections/impl/base_static.hpp>
42
 
#include <ggl/extensions/gis/projections/impl/base_dynamic.hpp>
43
 
#include <ggl/extensions/gis/projections/impl/projects.hpp>
44
 
#include <ggl/extensions/gis/projections/impl/factory_entry.hpp>
45
 
#include <ggl/extensions/gis/projections/impl/function_overloads.hpp>
46
 
#include <ggl/extensions/gis/projections/impl/pj_mlfn.hpp>
47
 
 
48
 
#include <ggl/extensions/gis/projections/epsg_traits.hpp>
49
 
 
50
 
namespace ggl { namespace projection
51
 
{
52
 
    #ifndef DOXYGEN_NO_DETAIL
53
 
    namespace detail { namespace tmerc{ 
54
 
            static const double EPS10 = 1.e-10;
55
 
            static const double FC1 = 1.;
56
 
            static const double FC2 = .5;
57
 
            static const double FC3 = .16666666666666666666;
58
 
            static const double FC4 = .08333333333333333333;
59
 
            static const double FC5 = .05;
60
 
            static const double FC6 = .03333333333333333333;
61
 
            static const double FC7 = .02380952380952380952;
62
 
            static const double FC8 = .01785714285714285714;
63
 
 
64
 
            struct par_tmerc
65
 
            {
66
 
                double    esp;
67
 
                double    ml0;
68
 
                double    en[EN_SIZE];
69
 
            };
70
 
 
71
 
            // template class, using CRTP to implement forward/inverse
72
 
            template <typename Geographic, typename Cartesian, typename Parameters>
73
 
            struct base_tmerc_ellipsoid : public base_t_fi<base_tmerc_ellipsoid<Geographic, Cartesian, Parameters>,
74
 
                     Geographic, Cartesian, Parameters>
75
 
            {
76
 
 
77
 
                 typedef double geographic_type;
78
 
                 typedef double cartesian_type;
79
 
 
80
 
                par_tmerc m_proj_parm;
81
 
 
82
 
                inline base_tmerc_ellipsoid(const Parameters& par)
83
 
                    : base_t_fi<base_tmerc_ellipsoid<Geographic, Cartesian, Parameters>,
84
 
                     Geographic, Cartesian, Parameters>(*this, par) {}
85
 
 
86
 
                inline void fwd(geographic_type& lp_lon, geographic_type& lp_lat, cartesian_type& xy_x, cartesian_type& xy_y) const
87
 
                {
88
 
                    double al, als, n, cosphi, sinphi, t;
89
 
                
90
 
                        /*
91
 
                         * Fail if our longitude is more than 90 degrees from the 
92
 
                         * central meridian since the results are essentially garbage. 
93
 
                         * Is error -20 really an appropriate return value?
94
 
                         * 
95
 
                         *  http://trac.osgeo.org/proj/ticket/5
96
 
                         */
97
 
                        if( lp_lon < -HALFPI || lp_lon > HALFPI )
98
 
                        {
99
 
                            xy_x = HUGE_VAL;
100
 
                            xy_y = HUGE_VAL;
101
 
                            throw proj_exception(  -14);
102
 
                            return;
103
 
                        }
104
 
                
105
 
                    sinphi = sin(lp_lat); cosphi = cos(lp_lat);
106
 
                    t = fabs(cosphi) > 1e-10 ? sinphi/cosphi : 0.;
107
 
                    t *= t;
108
 
                    al = cosphi * lp_lon;
109
 
                    als = al * al;
110
 
                    al /= sqrt(1. - this->m_par.es * sinphi * sinphi);
111
 
                    n = this->m_proj_parm.esp * cosphi * cosphi;
112
 
                    xy_x = this->m_par.k0 * al * (FC1 +
113
 
                        FC3 * als * (1. - t + n +
114
 
                        FC5 * als * (5. + t * (t - 18.) + n * (14. - 58. * t)
115
 
                        + FC7 * als * (61. + t * ( t * (179. - t) - 479. ) )
116
 
                        )));
117
 
                    xy_y = this->m_par.k0 * (pj_mlfn(lp_lat, sinphi, cosphi, this->m_proj_parm.en) - this->m_proj_parm.ml0 +
118
 
                        sinphi * al * lp_lon * FC2 * ( 1. +
119
 
                        FC4 * als * (5. - t + n * (9. + 4. * n) +
120
 
                        FC6 * als * (61. + t * (t - 58.) + n * (270. - 330 * t)
121
 
                        + FC8 * als * (1385. + t * ( t * (543. - t) - 3111.) )
122
 
                        ))));
123
 
                }
124
 
 
125
 
                inline void inv(cartesian_type& xy_x, cartesian_type& xy_y, geographic_type& lp_lon, geographic_type& lp_lat) const
126
 
                {
127
 
                    double n, con, cosphi, d, ds, sinphi, t;
128
 
                
129
 
                    lp_lat = pj_inv_mlfn(this->m_proj_parm.ml0 + xy_y / this->m_par.k0, this->m_par.es, this->m_proj_parm.en);
130
 
                    if (fabs(lp_lat) >= HALFPI) {
131
 
                        lp_lat = xy_y < 0. ? -HALFPI : HALFPI;
132
 
                        lp_lon = 0.;
133
 
                    } else {
134
 
                        sinphi = sin(lp_lat);
135
 
                        cosphi = cos(lp_lat);
136
 
                        t = fabs(cosphi) > 1e-10 ? sinphi/cosphi : 0.;
137
 
                        n = this->m_proj_parm.esp * cosphi * cosphi;
138
 
                        d = xy_x * sqrt(con = 1. - this->m_par.es * sinphi * sinphi) / this->m_par.k0;
139
 
                        con *= t;
140
 
                        t *= t;
141
 
                        ds = d * d;
142
 
                        lp_lat -= (con * ds / (1.-this->m_par.es)) * FC2 * (1. -
143
 
                            ds * FC4 * (5. + t * (3. - 9. *  n) + n * (1. - 4 * n) -
144
 
                            ds * FC6 * (61. + t * (90. - 252. * n +
145
 
                                45. * t) + 46. * n
146
 
                           - ds * FC8 * (1385. + t * (3633. + t * (4095. + 1574. * t)) )
147
 
                            )));
148
 
                        lp_lon = d*(FC1 -
149
 
                            ds*FC3*( 1. + 2.*t + n -
150
 
                            ds*FC5*(5. + t*(28. + 24.*t + 8.*n) + 6.*n
151
 
                           - ds * FC7 * (61. + t * (662. + t * (1320. + 720. * t)) )
152
 
                        ))) / cosphi;
153
 
                    }
154
 
                }
155
 
            };
156
 
 
157
 
            // template class, using CRTP to implement forward/inverse
158
 
            template <typename Geographic, typename Cartesian, typename Parameters>
159
 
            struct base_tmerc_spheroid : public base_t_fi<base_tmerc_spheroid<Geographic, Cartesian, Parameters>,
160
 
                     Geographic, Cartesian, Parameters>
161
 
            {
162
 
 
163
 
                 typedef double geographic_type;
164
 
                 typedef double cartesian_type;
165
 
 
166
 
                par_tmerc m_proj_parm;
167
 
 
168
 
                inline base_tmerc_spheroid(const Parameters& par)
169
 
                    : base_t_fi<base_tmerc_spheroid<Geographic, Cartesian, Parameters>,
170
 
                     Geographic, Cartesian, Parameters>(*this, par) {}
171
 
 
172
 
                inline void fwd(geographic_type& lp_lon, geographic_type& lp_lat, cartesian_type& xy_x, cartesian_type& xy_y) const
173
 
                {
174
 
                    double b, cosphi;
175
 
                
176
 
                        /*
177
 
                         * Fail if our longitude is more than 90 degrees from the 
178
 
                         * central meridian since the results are essentially garbage. 
179
 
                         * Is error -20 really an appropriate return value?
180
 
                         * 
181
 
                         *  http://trac.osgeo.org/proj/ticket/5
182
 
                         */
183
 
                        if( lp_lon < -HALFPI || lp_lon > HALFPI )
184
 
                        {
185
 
                            xy_x = HUGE_VAL;
186
 
                            xy_y = HUGE_VAL;
187
 
                            throw proj_exception(  -14);
188
 
                            return;
189
 
                        }
190
 
                
191
 
                    b = (cosphi = cos(lp_lat)) * sin(lp_lon);
192
 
                    if (fabs(fabs(b) - 1.) <= EPS10) throw proj_exception();;
193
 
                    xy_x = this->m_proj_parm.ml0 * log((1. + b) / (1. - b));
194
 
                    if ((b = fabs( xy_y = cosphi * cos(lp_lon) / sqrt(1. - b * b) )) >= 1.) {
195
 
                        if ((b - 1.) > EPS10) throw proj_exception();
196
 
                        else xy_y = 0.;
197
 
                    } else
198
 
                        xy_y = acos(xy_y);
199
 
                    if (lp_lat < 0.) xy_y = -xy_y;
200
 
                    xy_y = this->m_proj_parm.esp * (xy_y - this->m_par.phi0);
201
 
                }
202
 
 
203
 
                inline void inv(cartesian_type& xy_x, cartesian_type& xy_y, geographic_type& lp_lon, geographic_type& lp_lat) const
204
 
                {
205
 
                    double h, g;
206
 
                
207
 
                    h = exp(xy_x / this->m_proj_parm.esp);
208
 
                    g = .5 * (h - 1. / h);
209
 
                    h = cos(this->m_par.phi0 + xy_y / this->m_proj_parm.esp);
210
 
                    lp_lat = asin(sqrt((1. - h * h) / (1. + g * g)));
211
 
                    if (xy_y < 0.) lp_lat = -lp_lat;
212
 
                    lp_lon = (g || h) ? atan2(g, h) : 0.;
213
 
                }
214
 
            };
215
 
 
216
 
            template <typename Parameters>
217
 
            void setup(Parameters& par, par_tmerc& proj_parm)  /* general initialization */
218
 
            {
219
 
                boost::ignore_unused_variable_warning(par);
220
 
                boost::ignore_unused_variable_warning(proj_parm);
221
 
                if (par.es) {
222
 
                    pj_enfn(par.es, proj_parm.en);
223
 
            
224
 
                    proj_parm.ml0 = pj_mlfn(par.phi0, sin(par.phi0), cos(par.phi0), proj_parm.en);
225
 
                    proj_parm.esp = par.es / (1. - par.es);
226
 
                // par.inv = e_inverse;
227
 
                // par.fwd = e_forward;
228
 
                } else {
229
 
                    proj_parm.esp = par.k0;
230
 
                    proj_parm.ml0 = .5 * proj_parm.esp;
231
 
                // par.inv = s_inverse;
232
 
                // par.fwd = s_forward;
233
 
                }
234
 
            }
235
 
 
236
 
 
237
 
            // Transverse Mercator
238
 
            template <typename Parameters>
239
 
            void setup_tmerc(Parameters& par, par_tmerc& proj_parm)
240
 
            {
241
 
                setup(par, proj_parm);
242
 
            }
243
 
 
244
 
            // Universal Transverse Mercator (UTM)
245
 
            template <typename Parameters>
246
 
            void setup_utm(Parameters& par, par_tmerc& proj_parm)
247
 
            {
248
 
                int zone;
249
 
                if (!par.es) throw proj_exception(-34);
250
 
                par.y0 = pj_param(par.params, "bsouth").i ? 10000000. : 0.;
251
 
                par.x0 = 500000.;
252
 
                if (pj_param(par.params, "tzone").i) /* zone input ? */
253
 
                    if ((zone = pj_param(par.params, "izone").i) > 0 && zone <= 60)
254
 
                        --zone;
255
 
                    else
256
 
                        throw proj_exception(-35);
257
 
                else /* nearest central meridian input */
258
 
                    if ((zone = int_floor((adjlon(par.lam0) + PI) * 30. / PI)) < 0)
259
 
                        zone = 0;
260
 
                    else if (zone >= 60)
261
 
                        zone = 59;
262
 
                par.lam0 = (zone + .5) * PI / 30. - PI;
263
 
                par.k0 = 0.9996;
264
 
                par.phi0 = 0.;
265
 
                setup(par, proj_parm);
266
 
            }
267
 
 
268
 
        }} // namespace detail::tmerc
269
 
    #endif // doxygen 
270
 
 
271
 
    /*!
272
 
        \brief Transverse Mercator projection
273
 
        \ingroup projections
274
 
        \tparam Geographic latlong point type
275
 
        \tparam Cartesian xy point type
276
 
        \tparam Parameters parameter type
277
 
        \par Projection characteristics
278
 
         - Cylindrical
279
 
         - Spheroid
280
 
         - Ellipsoid
281
 
        \par Example
282
 
        \image html ex_tmerc.gif
283
 
    */
284
 
    template <typename Geographic, typename Cartesian, typename Parameters = parameters>
285
 
    struct tmerc_ellipsoid : public detail::tmerc::base_tmerc_ellipsoid<Geographic, Cartesian, Parameters>
286
 
    {
287
 
        inline tmerc_ellipsoid(const Parameters& par) : detail::tmerc::base_tmerc_ellipsoid<Geographic, Cartesian, Parameters>(par)
288
 
        {
289
 
            detail::tmerc::setup_tmerc(this->m_par, this->m_proj_parm);
290
 
        }
291
 
    };
292
 
 
293
 
    /*!
294
 
        \brief Universal Transverse Mercator (UTM) projection
295
 
        \ingroup projections
296
 
        \tparam Geographic latlong point type
297
 
        \tparam Cartesian xy point type
298
 
        \tparam Parameters parameter type
299
 
        \par Projection characteristics
300
 
         - Cylindrical
301
 
         - Spheroid
302
 
         - zone= south
303
 
        \par Example
304
 
        \image html ex_utm.gif
305
 
    */
306
 
    template <typename Geographic, typename Cartesian, typename Parameters = parameters>
307
 
    struct utm_ellipsoid : public detail::tmerc::base_tmerc_ellipsoid<Geographic, Cartesian, Parameters>
308
 
    {
309
 
        inline utm_ellipsoid(const Parameters& par) : detail::tmerc::base_tmerc_ellipsoid<Geographic, Cartesian, Parameters>(par)
310
 
        {
311
 
            detail::tmerc::setup_utm(this->m_par, this->m_proj_parm);
312
 
        }
313
 
    };
314
 
 
315
 
    /*!
316
 
        \brief Transverse Mercator projection
317
 
        \ingroup projections
318
 
        \tparam Geographic latlong point type
319
 
        \tparam Cartesian xy point type
320
 
        \tparam Parameters parameter type
321
 
        \par Projection characteristics
322
 
         - Cylindrical
323
 
         - Spheroid
324
 
         - Ellipsoid
325
 
        \par Example
326
 
        \image html ex_tmerc.gif
327
 
    */
328
 
    template <typename Geographic, typename Cartesian, typename Parameters = parameters>
329
 
    struct tmerc_spheroid : public detail::tmerc::base_tmerc_spheroid<Geographic, Cartesian, Parameters>
330
 
    {
331
 
        inline tmerc_spheroid(const Parameters& par) : detail::tmerc::base_tmerc_spheroid<Geographic, Cartesian, Parameters>(par)
332
 
        {
333
 
            detail::tmerc::setup_tmerc(this->m_par, this->m_proj_parm);
334
 
        }
335
 
    };
336
 
 
337
 
    #ifndef DOXYGEN_NO_DETAIL
338
 
    namespace detail
339
 
    {
340
 
 
341
 
        // Factory entry(s)
342
 
        template <typename Geographic, typename Cartesian, typename Parameters>
343
 
        class tmerc_entry : public detail::factory_entry<Geographic, Cartesian, Parameters>
344
 
        {
345
 
            public :
346
 
                virtual projection<Geographic, Cartesian>* create_new(const Parameters& par) const
347
 
                {
348
 
                    if (par.es)
349
 
                        return new base_v_fi<tmerc_ellipsoid<Geographic, Cartesian, Parameters>, Geographic, Cartesian, Parameters>(par);
350
 
                    else
351
 
                        return new base_v_fi<tmerc_spheroid<Geographic, Cartesian, Parameters>, Geographic, Cartesian, Parameters>(par);
352
 
                }
353
 
        };
354
 
 
355
 
        template <typename Geographic, typename Cartesian, typename Parameters>
356
 
        class utm_entry : public detail::factory_entry<Geographic, Cartesian, Parameters>
357
 
        {
358
 
            public :
359
 
                virtual projection<Geographic, Cartesian>* create_new(const Parameters& par) const
360
 
                {
361
 
                    return new base_v_fi<utm_ellipsoid<Geographic, Cartesian, Parameters>, Geographic, Cartesian, Parameters>(par);
362
 
                }
363
 
        };
364
 
 
365
 
        template <typename Geographic, typename Cartesian, typename Parameters>
366
 
        inline void tmerc_init(detail::base_factory<Geographic, Cartesian, Parameters>& factory)
367
 
        {
368
 
            factory.add_to_factory("tmerc", new tmerc_entry<Geographic, Cartesian, Parameters>);
369
 
            factory.add_to_factory("utm", new utm_entry<Geographic, Cartesian, Parameters>);
370
 
        }
371
 
 
372
 
    } // namespace detail 
373
 
    // Create EPSG specializations
374
 
    // (Proof of Concept, only for some)
375
 
 
376
 
    template<typename LatLongRadian, typename Cartesian, typename Parameters>
377
 
    struct epsg_traits<2000, LatLongRadian, Cartesian, Parameters>
378
 
    {
379
 
        typedef tmerc_ellipsoid<LatLongRadian, Cartesian, Parameters> type;
380
 
        static inline std::string par()
381
 
        {
382
 
            return "+proj=tmerc +lat_0=0 +lon_0=-62 +k=0.9995000000000001 +x_0=400000 +y_0=0 +ellps=clrk80 +units=m";
383
 
        }
384
 
    };
385
 
 
386
 
 
387
 
    template<typename LatLongRadian, typename Cartesian, typename Parameters>
388
 
    struct epsg_traits<2001, LatLongRadian, Cartesian, Parameters>
389
 
    {
390
 
        typedef tmerc_ellipsoid<LatLongRadian, Cartesian, Parameters> type;
391
 
        static inline std::string par()
392
 
        {
393
 
            return "+proj=tmerc +lat_0=0 +lon_0=-62 +k=0.9995000000000001 +x_0=400000 +y_0=0 +ellps=clrk80 +units=m";
394
 
        }
395
 
    };
396
 
 
397
 
 
398
 
    template<typename LatLongRadian, typename Cartesian, typename Parameters>
399
 
    struct epsg_traits<2002, LatLongRadian, Cartesian, Parameters>
400
 
    {
401
 
        typedef tmerc_ellipsoid<LatLongRadian, Cartesian, Parameters> type;
402
 
        static inline std::string par()
403
 
        {
404
 
            return "+proj=tmerc +lat_0=0 +lon_0=-62 +k=0.9995000000000001 +x_0=400000 +y_0=0 +ellps=clrk80 +towgs84=725,685,536,0,0,0,0 +units=m";
405
 
        }
406
 
    };
407
 
 
408
 
 
409
 
    template<typename LatLongRadian, typename Cartesian, typename Parameters>
410
 
    struct epsg_traits<2003, LatLongRadian, Cartesian, Parameters>
411
 
    {
412
 
        typedef tmerc_ellipsoid<LatLongRadian, Cartesian, Parameters> type;
413
 
        static inline std::string par()
414
 
        {
415
 
            return "+proj=tmerc +lat_0=0 +lon_0=-62 +k=0.9995000000000001 +x_0=400000 +y_0=0 +ellps=clrk80 +towgs84=72,213.7,93,0,0,0,0 +units=m";
416
 
        }
417
 
    };
418
 
 
419
 
 
420
 
    template<typename LatLongRadian, typename Cartesian, typename Parameters>
421
 
    struct epsg_traits<2039, LatLongRadian, Cartesian, Parameters>
422
 
    {
423
 
        typedef tmerc_ellipsoid<LatLongRadian, Cartesian, Parameters> type;
424
 
        static inline std::string par()
425
 
        {
426
 
            return "+proj=tmerc +lat_0=31.73439361111111 +lon_0=35.20451694444445 +k=1.0000067 +x_0=219529.584 +y_0=626907.39 +ellps=GRS80 +towgs84=-48,55,52,0,0,0,0 +units=m";
427
 
        }
428
 
    };
429
 
 
430
 
 
431
 
    template<typename LatLongRadian, typename Cartesian, typename Parameters>
432
 
    struct epsg_traits<29118, LatLongRadian, Cartesian, Parameters>
433
 
    {
434
 
        typedef utm_ellipsoid<LatLongRadian, Cartesian, Parameters> type;
435
 
        static inline std::string par()
436
 
        {
437
 
            return "+proj=utm +zone=18 +ellps=GRS67 +units=m";
438
 
        }
439
 
    };
440
 
 
441
 
 
442
 
    template<typename LatLongRadian, typename Cartesian, typename Parameters>
443
 
    struct epsg_traits<29119, LatLongRadian, Cartesian, Parameters>
444
 
    {
445
 
        typedef utm_ellipsoid<LatLongRadian, Cartesian, Parameters> type;
446
 
        static inline std::string par()
447
 
        {
448
 
            return "+proj=utm +zone=19 +ellps=GRS67 +units=m";
449
 
        }
450
 
    };
451
 
 
452
 
 
453
 
    #endif // doxygen
454
 
 
455
 
}} // namespace ggl::projection
456
 
 
457
 
#endif // GGL_PROJECTIONS_TMERC_HPP
458