← 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/krovak.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/krovak.hpp
1
 
#ifndef GGL_PROJECTIONS_KROVAK_HPP
2
 
#define GGL_PROJECTIONS_KROVAK_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/math/special_functions/hypot.hpp>
39
 
 
40
 
#include <ggl/extensions/gis/projections/impl/base_static.hpp>
41
 
#include <ggl/extensions/gis/projections/impl/base_dynamic.hpp>
42
 
#include <ggl/extensions/gis/projections/impl/projects.hpp>
43
 
#include <ggl/extensions/gis/projections/impl/factory_entry.hpp>
44
 
 
45
 
namespace ggl { namespace projection
46
 
{
47
 
    #ifndef DOXYGEN_NO_DETAIL
48
 
    namespace detail { namespace krovak{ 
49
 
 
50
 
            struct par_krovak
51
 
            {
52
 
                double    C_x;
53
 
            };
54
 
            
55
 
            
56
 
            
57
 
            
58
 
            
59
 
            /**
60
 
               NOTES: According to EPSG the full Krovak projection method should have
61
 
                      the following parameters.  Within PROJ.4 the azimuth, and pseudo
62
 
                      standard parallel are hardcoded in the algorithm and can't be 
63
 
                      altered from outside.  The others all have defaults to match the
64
 
                      common usage with Krovak projection.
65
 
            
66
 
              lat_0 = latitude of centre of the projection
67
 
                     
68
 
              lon_0 = longitude of centre of the projection
69
 
              
70
 
              ** = azimuth (true) of the centre line passing through the centre of the projection
71
 
            
72
 
              ** = latitude of pseudo standard parallel
73
 
               
74
 
              k  = scale factor on the pseudo standard parallel
75
 
              
76
 
              x_0 = False Easting of the centre of the projection at the apex of the cone
77
 
              
78
 
              y_0 = False Northing of the centre of the projection at the apex of the cone
79
 
            
80
 
             **/
81
 
            
82
 
            
83
 
            
84
 
 
85
 
            // template class, using CRTP to implement forward/inverse
86
 
            template <typename Geographic, typename Cartesian, typename Parameters>
87
 
            struct base_krovak_ellipsoid : public base_t_fi<base_krovak_ellipsoid<Geographic, Cartesian, Parameters>,
88
 
                     Geographic, Cartesian, Parameters>
89
 
            {
90
 
 
91
 
                 typedef double geographic_type;
92
 
                 typedef double cartesian_type;
93
 
 
94
 
                par_krovak m_proj_parm;
95
 
 
96
 
                inline base_krovak_ellipsoid(const Parameters& par)
97
 
                    : base_t_fi<base_krovak_ellipsoid<Geographic, Cartesian, Parameters>,
98
 
                     Geographic, Cartesian, Parameters>(*this, par) {}
99
 
 
100
 
                inline void fwd(geographic_type& lp_lon, geographic_type& lp_lat, cartesian_type& xy_x, cartesian_type& xy_y) const
101
 
                {
102
 
                /* calculate xy from lat/lon */
103
 
                
104
 
                    
105
 
                    
106
 
                
107
 
                /* Constants, identical to inverse transform function */
108
 
                    double s45, s90, e2, e, alfa, uq, u0, g, k, k1, n0, ro0, ad, a, s0, n;
109
 
                    double gfi, u, fi0, deltav, s, d, eps, ro;
110
 
                
111
 
                
112
 
                    s45 = 0.785398163397448;    /* 45 DEG */
113
 
                    s90 = 2 * s45;
114
 
                    fi0 = this->m_par.phi0;    /* Latitude of projection centre 49 DEG 30' */
115
 
                
116
 
                   /* Ellipsoid is used as Parameter in for.c and inv.c, therefore a must 
117
 
                      be set to 1 here.
118
 
                      Ellipsoid Bessel 1841 a = 6377397.155m 1/f = 299.1528128,
119
 
                      e2=0.006674372230614;
120
 
                   */
121
 
                    a =  1; /* 6377397.155; */
122
 
                    /* e2 = this->m_par.es;*/       /* 0.006674372230614; */
123
 
                    e2 = 0.006674372230614;
124
 
                    e = sqrt(e2);
125
 
                
126
 
                    alfa = sqrt(1. + (e2 * pow(cos(fi0), 4)) / (1. - e2));
127
 
                
128
 
                    uq = 1.04216856380474;      /* DU(2, 59, 42, 42.69689) */
129
 
                    u0 = asin(sin(fi0) / alfa);
130
 
                    g = pow(   (1. + e * sin(fi0)) / (1. - e * sin(fi0)) , alfa * e / 2.  );
131
 
                
132
 
                    k = tan( u0 / 2. + s45) / pow  (tan(fi0 / 2. + s45) , alfa) * g;
133
 
                
134
 
                    k1 = this->m_par.k0;
135
 
                    n0 = a * sqrt(1. - e2) / (1. - e2 * pow(sin(fi0), 2));
136
 
                    s0 = 1.37008346281555;       /* Latitude of pseudo standard parallel 78 DEG 30'00" N */
137
 
                    n = sin(s0);
138
 
                    ro0 = k1 * n0 / tan(s0);
139
 
                    ad = s90 - uq;
140
 
                
141
 
                /* Transformation */
142
 
                
143
 
                    gfi =pow ( ((1. + e * sin(lp_lat)) /
144
 
                               (1. - e * sin(lp_lat))) , (alfa * e / 2.));
145
 
                
146
 
                    u= 2. * (atan(k * pow( tan(lp_lat / 2. + s45), alfa) / gfi)-s45);
147
 
                
148
 
                    deltav = - lp_lon * alfa;
149
 
                
150
 
                    s = asin(cos(ad) * sin(u) + sin(ad) * cos(u) * cos(deltav));
151
 
                    d = asin(cos(u) * sin(deltav) / cos(s));
152
 
                    eps = n * d;
153
 
                    ro = ro0 * pow(tan(s0 / 2. + s45) , n) / pow(tan(s / 2. + s45) , n)   ;
154
 
                
155
 
                   /* x and y are reverted! */
156
 
                    xy_y = ro * cos(eps) / a;
157
 
                    xy_x = ro * sin(eps) / a;
158
 
                
159
 
                        if( !pj_param(this->m_par.params, "tczech").i )
160
 
                      {
161
 
                        xy_y *= -1.0;
162
 
                        xy_x *= -1.0;
163
 
                      }
164
 
                
165
 
                            return;
166
 
                }
167
 
                
168
 
                
169
 
                
170
 
 
171
 
                inline void inv(cartesian_type& xy_x, cartesian_type& xy_y, geographic_type& lp_lon, geographic_type& lp_lat) const
172
 
                {
173
 
                    /* calculate lat/lon from xy */
174
 
                
175
 
                /* Constants, identisch wie in der Umkehrfunktion */
176
 
                    double s45, s90, fi0, e2, e, alfa, uq, u0, g, k, k1, n0, ro0, ad, a, s0, n;
177
 
                    double u, deltav, s, d, eps, ro, fi1, xy0;
178
 
                    int ok;
179
 
                
180
 
                    s45 = 0.785398163397448;    /* 45 DEG */
181
 
                    s90 = 2 * s45;
182
 
                    fi0 = this->m_par.phi0;    /* Latitude of projection centre 49 DEG 30' */
183
 
                
184
 
                
185
 
                   /* Ellipsoid is used as Parameter in for.c and inv.c, therefore a must 
186
 
                      be set to 1 here.
187
 
                      Ellipsoid Bessel 1841 a = 6377397.155m 1/f = 299.1528128,
188
 
                      e2=0.006674372230614;
189
 
                   */
190
 
                    a = 1; /* 6377397.155; */
191
 
                    /* e2 = this->m_par.es; */      /* 0.006674372230614; */
192
 
                    e2 = 0.006674372230614;
193
 
                    e = sqrt(e2);
194
 
                
195
 
                    alfa = sqrt(1. + (e2 * pow(cos(fi0), 4)) / (1. - e2));
196
 
                    uq = 1.04216856380474;      /* DU(2, 59, 42, 42.69689) */
197
 
                    u0 = asin(sin(fi0) / alfa);
198
 
                    g = pow(   (1. + e * sin(fi0)) / (1. - e * sin(fi0)) , alfa * e / 2.  );
199
 
                
200
 
                    k = tan( u0 / 2. + s45) / pow  (tan(fi0 / 2. + s45) , alfa) * g;
201
 
                
202
 
                    k1 = this->m_par.k0;
203
 
                    n0 = a * sqrt(1. - e2) / (1. - e2 * pow(sin(fi0), 2));
204
 
                    s0 = 1.37008346281555;       /* Latitude of pseudo standard parallel 78 DEG 30'00" N */
205
 
                    n = sin(s0);
206
 
                    ro0 = k1 * n0 / tan(s0);
207
 
                    ad = s90 - uq;
208
 
                
209
 
                
210
 
                /* Transformation */
211
 
                   /* revert y, x*/
212
 
                    xy0=xy_x;
213
 
                    xy_x=xy_y;
214
 
                    xy_y=xy0;
215
 
                
216
 
                        if( !pj_param(this->m_par.params, "tczech").i )
217
 
                      {
218
 
                        xy_x *= -1.0;
219
 
                        xy_y *= -1.0;
220
 
                      }
221
 
                
222
 
                    ro = sqrt(xy_x * xy_x + xy_y * xy_y);
223
 
                    eps = atan2(xy_y, xy_x);
224
 
                    d = eps / sin(s0);
225
 
                    s = 2. * (atan(  pow(ro0 / ro, 1. / n) * tan(s0 / 2. + s45)) - s45);
226
 
                
227
 
                    u = asin(cos(ad) * sin(s) - sin(ad) * cos(s) * cos(d));
228
 
                    deltav = asin(cos(s) * sin(d) / cos(u));
229
 
                
230
 
                    lp_lon = this->m_par.lam0 - deltav / alfa;
231
 
                
232
 
                /* ITERATION FOR lp_lat */
233
 
                   fi1 = u;
234
 
                
235
 
                   ok = 0;
236
 
                   do
237
 
                   {
238
 
                       lp_lat = 2. * ( atan( pow( k, -1. / alfa)  *
239
 
                                            pow( tan(u / 2. + s45) , 1. / alfa)  *
240
 
                                            pow( (1. + e * sin(fi1)) / (1. - e * sin(fi1)) , e / 2.)
241
 
                                           )  - s45);
242
 
                
243
 
                      if (fabs(fi1 - lp_lat) < 0.000000000000001) ok=1;
244
 
                      fi1 = lp_lat;
245
 
                
246
 
                   }
247
 
                   while (ok==0);
248
 
                
249
 
                   lp_lon -= this->m_par.lam0;
250
 
                
251
 
                            return;
252
 
                }
253
 
                
254
 
            };
255
 
 
256
 
            // Krovak
257
 
            template <typename Parameters>
258
 
            void setup_krovak(Parameters& par, par_krovak& proj_parm)
259
 
            {
260
 
                double ts;
261
 
                /* read some Parameters,
262
 
                 * here Latitude Truescale */
263
 
                ts = pj_param(par.params, "rlat_ts").f;
264
 
                proj_parm.C_x = ts;
265
 
                
266
 
                /* we want Bessel as fixed ellipsoid */
267
 
                par.a = 6377397.155;
268
 
                par.e = sqrt(par.es = 0.006674372230614);
269
 
                    /* if latitude of projection center is not set, use 49d30'N */
270
 
                if (!pj_param(par.params, "tlat_0").i)
271
 
                        par.phi0 = 0.863937979737193;
272
 
             
273
 
                    /* if center long is not set use 42d30'E of Ferro - 17d40' for Ferro */
274
 
                    /* that will correspond to using longitudes relative to greenwich    */
275
 
                    /* as input and output, instead of lat/long relative to Ferro */
276
 
                if (!pj_param(par.params, "tlon_0").i)
277
 
                        par.lam0 = 0.7417649320975901 - 0.308341501185665;
278
 
                    /* if scale not set default to 0.9999 */
279
 
                if (!pj_param(par.params, "tk").i)
280
 
                        par.k0 = 0.9999;
281
 
                /* always the same */
282
 
                // par.inv = e_inverse;
283
 
             
284
 
                // par.fwd = e_forward;
285
 
            }
286
 
 
287
 
        }} // namespace detail::krovak
288
 
    #endif // doxygen 
289
 
 
290
 
    /*!
291
 
        \brief Krovak projection
292
 
        \ingroup projections
293
 
        \tparam Geographic latlong point type
294
 
        \tparam Cartesian xy point type
295
 
        \tparam Parameters parameter type
296
 
        \par Projection characteristics
297
 
         - Pseudocylindrical
298
 
         - Ellps
299
 
        \par Example
300
 
        \image html ex_krovak.gif
301
 
    */
302
 
    template <typename Geographic, typename Cartesian, typename Parameters = parameters>
303
 
    struct krovak_ellipsoid : public detail::krovak::base_krovak_ellipsoid<Geographic, Cartesian, Parameters>
304
 
    {
305
 
        inline krovak_ellipsoid(const Parameters& par) : detail::krovak::base_krovak_ellipsoid<Geographic, Cartesian, Parameters>(par)
306
 
        {
307
 
            detail::krovak::setup_krovak(this->m_par, this->m_proj_parm);
308
 
        }
309
 
    };
310
 
 
311
 
    #ifndef DOXYGEN_NO_DETAIL
312
 
    namespace detail
313
 
    {
314
 
 
315
 
        // Factory entry(s)
316
 
        template <typename Geographic, typename Cartesian, typename Parameters>
317
 
        class krovak_entry : public detail::factory_entry<Geographic, Cartesian, Parameters>
318
 
        {
319
 
            public :
320
 
                virtual projection<Geographic, Cartesian>* create_new(const Parameters& par) const
321
 
                {
322
 
                    return new base_v_fi<krovak_ellipsoid<Geographic, Cartesian, Parameters>, Geographic, Cartesian, Parameters>(par);
323
 
                }
324
 
        };
325
 
 
326
 
        template <typename Geographic, typename Cartesian, typename Parameters>
327
 
        inline void krovak_init(detail::base_factory<Geographic, Cartesian, Parameters>& factory)
328
 
        {
329
 
            factory.add_to_factory("krovak", new krovak_entry<Geographic, Cartesian, Parameters>);
330
 
        }
331
 
 
332
 
    } // namespace detail 
333
 
    #endif // doxygen
334
 
 
335
 
}} // namespace ggl::projection
336
 
 
337
 
#endif // GGL_PROJECTIONS_KROVAK_HPP
338