Back to home page

EIC code displayed by LXR
 
 

    

File indexing completed on 2025-01-18 09:35:39 

0001 // Boost.Geometry
0002 // This file is manually converted from PROJ4
0003 
0004 // Copyright (c) 2023 Adam Wulkiewicz, Lodz, Poland.
0005 
0006 // This file was modified by Oracle on 2017, 2018.
0007 // Modifications copyright (c) 2017-2018, Oracle and/or its affiliates.
0008 // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
0009 
0010 // Use, modification and distribution is subject to the Boost Software License,
0011 // Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
0012 // http://www.boost.org/LICENSE_1_0.txt)
0013 
0014 // This file is converted from PROJ4, http://trac.osgeo.org/proj
0015 // PROJ4 is originally written by Gerald Evenden (then of the USGS)
0016 // PROJ4 is maintained by Frank Warmerdam
0017 // This file was converted to Geometry Library by Adam Wulkiewicz
0018 
0019 // Original copyright notice:
0020 
0021 // Copyright (c) 2000, Frank Warmerdam
0022 
0023 // Permission is hereby granted, free of charge, to any person obtaining a
0024 // copy of this software and associated documentation files (the "Software"),
0025 // to deal in the Software without restriction, including without limitation
0026 // the rights to use, copy, modify, merge, publish, distribute, sublicense,
0027 // and/or sell copies of the Software, and to permit persons to whom the
0028 // Software is furnished to do so, subject to the following conditions:
0029 
0030 // The above copyright notice and this permission notice shall be included
0031 // in all copies or substantial portions of the Software.
0032 
0033 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
0034 // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
0035 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
0036 // THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
0037 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
0038 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
0039 // DEALINGS IN THE SOFTWARE.
0040 
0041 #ifndef BOOST_GEOMETRY_SRS_PROJECTIONS_IMPL_PJ_TRANSFORM_HPP
0042 #define BOOST_GEOMETRY_SRS_PROJECTIONS_IMPL_PJ_TRANSFORM_HPP
0043 
0044 
0045 #include <boost/geometry/core/access.hpp>
0046 #include <boost/geometry/core/coordinate_dimension.hpp>
0047 #include <boost/geometry/core/radian_access.hpp>
0048 
0049 #include <boost/geometry/srs/projections/impl/geocent.hpp>
0050 #include <boost/geometry/srs/projections/impl/pj_apply_gridshift.hpp>
0051 #include <boost/geometry/srs/projections/impl/projects.hpp>
0052 #include <boost/geometry/srs/projections/invalid_point.hpp>
0053 
0054 #include <boost/geometry/util/condition.hpp>
0055 #include <boost/geometry/util/range.hpp>
0056 
0057 #include <cstring>
0058 #include <cmath>
0059 
0060 
0061 namespace boost { namespace geometry { namespace projections
0062 {
0063 
0064 namespace detail
0065 {
0066 
0067 // -----------------------------------------------------------
0068 // Boost.Geometry helpers begin
0069 // -----------------------------------------------------------
0070 
0071 template
0072 <
0073     typename Point,
0074     bool HasCoord2 = (dimension<Point>::value > 2)
0075 >
0076 struct z_access
0077 {
0078     typedef typename coordinate_type<Point>::type type;
0079     static inline type get(Point const& point)
0080     {
0081         return geometry::get<2>(point);
0082     }
0083     static inline void set(Point & point, type const& h)
0084     {
0085         return geometry::set<2>(point, h);
0086     }
0087 };
0088 
0089 template <typename Point>
0090 struct z_access<Point, false>
0091 {
0092     typedef typename coordinate_type<Point>::type type;
0093     static inline type get(Point const& )
0094     {
0095         return type(0);
0096     }
0097     static inline void set(Point & , type const& )
0098     {}
0099 };
0100 
0101 template <typename XYorXYZ>
0102 inline typename z_access<XYorXYZ>::type
0103 get_z(XYorXYZ const& xy_or_xyz)
0104 {
0105     return z_access<XYorXYZ>::get(xy_or_xyz);
0106 }
0107 
0108 template <typename XYorXYZ>
0109 inline void set_z(XYorXYZ & xy_or_xyz,
0110                   typename z_access<XYorXYZ>::type const& z)
0111 {
0112     return z_access<XYorXYZ>::set(xy_or_xyz, z);
0113 }
0114 
0115 template
0116 <
0117     typename Range,
0118     bool AddZ = (dimension<typename boost::range_value<Range>::type>::value < 3)
0119 >
0120 struct range_wrapper
0121 {
0122     typedef Range range_type;
0123     typedef typename boost::range_value<Range>::type point_type;
0124     typedef typename coordinate_type<point_type>::type coord_t;
0125 
0126     range_wrapper(Range & range)
0127         : m_range(range)
0128     {}
0129 
0130     range_type & get_range() { return m_range; }
0131 
0132     coord_t get_z(std::size_t i) { return detail::get_z(range::at(m_range, i)); }
0133     void set_z(std::size_t i, coord_t const& v) { return detail::set_z(range::at(m_range, i), v); }
0134 
0135 private:
0136     Range & m_range;
0137 };
0138 
0139 template <typename Range>
0140 struct range_wrapper<Range, true>
0141 {
0142     typedef Range range_type;
0143     typedef typename boost::range_value<Range>::type point_type;
0144     typedef typename coordinate_type<point_type>::type coord_t;
0145 
0146     range_wrapper(Range & range)
0147         : m_range(range)
0148         , m_zs(boost::size(range), coord_t(0))
0149     {}
0150 
0151     range_type & get_range() { return m_range; }
0152 
0153     coord_t get_z(std::size_t i) { return m_zs[i]; }
0154     void set_z(std::size_t i, coord_t const& v) { m_zs[i] = v; }
0155 
0156 private:
0157     Range & m_range;
0158     std::vector<coord_t> m_zs;
0159 };
0160 
0161 // -----------------------------------------------------------
0162 // Boost.Geometry helpers end
0163 // -----------------------------------------------------------
0164 
0165 template <typename Par>
0166 inline typename Par::type Dx_BF(Par const& defn) { return defn.datum_params[0]; }
0167 template <typename Par>
0168 inline typename Par::type Dy_BF(Par const& defn) { return defn.datum_params[1]; }
0169 template <typename Par>
0170 inline typename Par::type Dz_BF(Par const& defn) { return defn.datum_params[2]; }
0171 template <typename Par>
0172 inline typename Par::type Rx_BF(Par const& defn) { return defn.datum_params[3]; }
0173 template <typename Par>
0174 inline typename Par::type Ry_BF(Par const& defn) { return defn.datum_params[4]; }
0175 template <typename Par>
0176 inline typename Par::type Rz_BF(Par const& defn) { return defn.datum_params[5]; }
0177 template <typename Par>
0178 inline typename Par::type M_BF(Par const& defn) { return defn.datum_params[6]; }
0179 
0180 /*
0181 ** This table is intended to indicate for any given error code in
0182 ** the range 0 to -56, whether that error will occur for all locations (ie.
0183 ** it is a problem with the coordinate system as a whole) in which case the
0184 ** value would be 0, or if the problem is with the point being transformed
0185 ** in which case the value is 1.
0186 **
0187 ** At some point we might want to move this array in with the error message
0188 ** list or something, but while experimenting with it this should be fine.
0189 **
0190 **
0191 ** NOTE (2017-10-01): Non-transient errors really should have resulted in a
0192 ** PJ==0 during initialization, and hence should be handled at the level
0193 ** before calling pj_transform. The only obvious example of the contrary
0194 ** appears to be the PJD_ERR_GRID_AREA case, which may also be taken to
0195 ** mean "no grids available"
0196 **
0197 **
0198 */
0199 
0200 static const int transient_error[60] = {
0201     /*             0  1  2  3  4  5  6  7  8  9   */
0202     /* 0 to 9 */   0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0203     /* 10 to 19 */ 0, 0, 0, 0, 1, 1, 0, 1, 1, 1,
0204     /* 20 to 29 */ 1, 0, 0, 0, 0, 0, 0, 1, 0, 0,
0205     /* 30 to 39 */ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0206     /* 40 to 49 */ 0, 0, 0, 0, 0, 0, 0, 0, 1, 0,
0207     /* 50 to 59 */ 1, 0, 1, 0, 1, 1, 1, 1, 0, 0 };
0208 
0209 
0210 template <typename T, typename Range>
0211 inline int pj_geocentric_to_geodetic( T const& a, T const& es,
0212                                       Range & range );
0213 template <typename T, typename Range>
0214 inline int pj_geodetic_to_geocentric( T const& a, T const& es,
0215                                       Range & range );
0216 
0217 /************************************************************************/
0218 /*                            pj_transform()                            */
0219 /*                                                                      */
0220 /*      Currently this function doesn't recognise if two projections    */
0221 /*      are identical (to short circuit reprojection) because it is     */
0222 /*      difficult to compare PJ structures (since there are some        */
0223 /*      projection specific components).                                */
0224 /************************************************************************/
0225 
0226 template <
0227     typename SrcPrj,
0228     typename DstPrj2,
0229     typename Par,
0230     typename Range,
0231     typename Grids
0232 >
0233 inline bool pj_transform(SrcPrj const& srcprj, Par const& srcdefn,
0234                          DstPrj2 const& dstprj, Par const& dstdefn,
0235                          Range & range,
0236                          Grids const& srcgrids,
0237                          Grids const& dstgrids)
0238 
0239 {
0240     typedef typename boost::range_value<Range>::type point_type;
0241     typedef typename coordinate_type<point_type>::type coord_t;
0242     static const std::size_t dimension = geometry::dimension<point_type>::value;
0243     std::size_t point_count = boost::size(range);
0244     bool result = true;
0245 
0246 /* -------------------------------------------------------------------- */
0247 /*      Transform unusual input coordinate axis orientation to          */
0248 /*      standard form if needed.                                        */
0249 /* -------------------------------------------------------------------- */
0250     // Ignored
0251 
0252 /* -------------------------------------------------------------------- */
0253 /*      Transform Z to meters if it isn't already.                      */
0254 /* -------------------------------------------------------------------- */
0255     if( BOOST_GEOMETRY_CONDITION(srcdefn.vto_meter != 1.0 && dimension > 2) )
0256     {
0257         for( std::size_t i = 0; i < point_count; i++ )
0258         {
0259             point_type & point = geometry::range::at(range, i);
0260             set_z(point, get_z(point) * srcdefn.vto_meter);
0261         }
0262     }
0263 
0264 /* -------------------------------------------------------------------- */
0265 /*      Transform geocentric source coordinates to lat/long.            */
0266 /* -------------------------------------------------------------------- */
0267     if( BOOST_GEOMETRY_CONDITION(srcdefn.is_geocent) )
0268     {
0269         // Point should be cartesian 3D (ECEF)
0270         if ( BOOST_GEOMETRY_CONDITION(dimension < 3) )
0271             BOOST_THROW_EXCEPTION( projection_exception(error_geocentric) );
0272             //return PJD_ERR_GEOCENTRIC;
0273 
0274         if( srcdefn.to_meter != 1.0 )
0275         {
0276             for(std::size_t i = 0; i < point_count; i++ )
0277             {
0278                 point_type & point = range::at(range, i);
0279                 if( ! is_invalid_point(point) )
0280                 {
0281                     set<0>(point, get<0>(point) * srcdefn.to_meter);
0282                     set<1>(point, get<1>(point) * srcdefn.to_meter);
0283                 }
0284             }
0285         }
0286 
0287         range_wrapper<Range, false> rng(range);
0288         int err = pj_geocentric_to_geodetic( srcdefn.a_orig, srcdefn.es_orig,
0289                                              rng );
0290         if( err != 0 )
0291             BOOST_THROW_EXCEPTION( projection_exception(err) );
0292             //return err;
0293 
0294         // NOTE: here 3D cartesian ECEF is converted into 3D geodetic LLH
0295     }
0296 
0297 /* -------------------------------------------------------------------- */
0298 /*      Transform source points to lat/long, if they aren't             */
0299 /*      already.                                                        */
0300 /* -------------------------------------------------------------------- */
0301     else if( !srcdefn.is_latlong )
0302     {
0303         // Point should be cartesian 2D or 3D (map projection)
0304 
0305         /* Check first if projection is invertible. */
0306         /*if( (srcdefn->inv3d == NULL) && (srcdefn->inv == NULL))
0307         {
0308             pj_ctx_set_errno( pj_get_ctx(srcdefn), -17 );
0309             pj_log( pj_get_ctx(srcdefn), PJ_LOG_ERROR,
0310                     "pj_transform(): source projection not invertible" );
0311             return -17;
0312         }*/
0313 
0314         /* If invertible - First try inv3d if defined */
0315         //if (srcdefn->inv3d != NULL)
0316         //{
0317         //    /* Three dimensions must be defined */
0318         //    if ( z == NULL)
0319         //    {
0320         //        pj_ctx_set_errno( pj_get_ctx(srcdefn), PJD_ERR_GEOCENTRIC);
0321         //        return PJD_ERR_GEOCENTRIC;
0322         //    }
0323 
0324         //    for (i=0; i < point_count; i++)
0325         //    {
0326         //        XYZ projected_loc;
0327         //        XYZ geodetic_loc;
0328 
0329         //        projected_loc.u = x[point_offset*i];
0330         //        projected_loc.v = y[point_offset*i];
0331         //        projected_loc.w = z[point_offset*i];
0332 
0333         //        if (projected_loc.u == HUGE_VAL)
0334         //            continue;
0335 
0336         //        geodetic_loc = pj_inv3d(projected_loc, srcdefn);
0337         //        if( srcdefn->ctx->last_errno != 0 )
0338         //        {
0339         //            if( (srcdefn->ctx->last_errno != 33 /*EDOM*/
0340         //                 && srcdefn->ctx->last_errno != 34 /*ERANGE*/ )
0341         //                && (srcdefn->ctx->last_errno > 0
0342         //                    || srcdefn->ctx->last_errno < -44 || point_count == 1
0343         //                    || transient_error[-srcdefn->ctx->last_errno] == 0 ) )
0344         //                return srcdefn->ctx->last_errno;
0345         //            else
0346         //            {
0347         //                geodetic_loc.u = HUGE_VAL;
0348         //                geodetic_loc.v = HUGE_VAL;
0349         //                geodetic_loc.w = HUGE_VAL;
0350         //            }
0351         //        }
0352 
0353         //        x[point_offset*i] = geodetic_loc.u;
0354         //        y[point_offset*i] = geodetic_loc.v;
0355         //        z[point_offset*i] = geodetic_loc.w;
0356 
0357         //    }
0358 
0359         //}
0360         //else
0361         {
0362             /* Fallback to the original PROJ.4 API 2d inversion - inv */
0363             for( std::size_t i = 0; i < point_count; i++ )
0364             {
0365                 point_type & point = range::at(range, i);
0366 
0367                 if( is_invalid_point(point) )
0368                     continue;
0369 
0370                 try
0371                 {
0372                     pj_inv(srcprj, srcdefn, point, point);
0373                 }
0374                 catch(projection_exception const& e)
0375                 {
0376                     if( (e.code() != 33 /*EDOM*/
0377                         && e.code() != 34 /*ERANGE*/ )
0378                         && (e.code() > 0
0379                             || e.code() < -44 /*|| point_count == 1*/
0380                             || transient_error[-e.code()] == 0) ) {
0381                         BOOST_RETHROW
0382                     } else {
0383                         set_invalid_point(point);
0384                         result = false;
0385                         if (point_count == 1)
0386                             return result;
0387                     }
0388                 }
0389             }
0390         }
0391     }
0392 
0393 /* -------------------------------------------------------------------- */
0394 /*      But if they are already lat long, adjust for the prime          */
0395 /*      meridian if there is one in effect.                             */
0396 /* -------------------------------------------------------------------- */
0397     if( srcdefn.from_greenwich != 0.0 )
0398     {
0399         for( std::size_t i = 0; i < point_count; i++ )
0400         {
0401             point_type & point = range::at(range, i);
0402 
0403             if( ! is_invalid_point(point) )
0404                 set<0>(point, get<0>(point) + srcdefn.from_greenwich);
0405         }
0406     }
0407 
0408 /* -------------------------------------------------------------------- */
0409 /*      Do we need to translate from geoid to ellipsoidal vertical      */
0410 /*      datum?                                                          */
0411 /* -------------------------------------------------------------------- */
0412     /*if( srcdefn->has_geoid_vgrids && z != NULL )
0413     {
0414         if( pj_apply_vgridshift( srcdefn, "sgeoidgrids",
0415                                  &(srcdefn->vgridlist_geoid),
0416                                  &(srcdefn->vgridlist_geoid_count),
0417                                  0, point_count, point_offset, x, y, z ) != 0 )
0418             return pj_ctx_get_errno(srcdefn->ctx);
0419     }*/
0420 
0421 /* -------------------------------------------------------------------- */
0422 /*      Convert datums if needed, and possible.                         */
0423 /* -------------------------------------------------------------------- */
0424     if ( ! pj_datum_transform( srcdefn, dstdefn, range, srcgrids, dstgrids ) )
0425     {
0426         result = false;
0427     }
0428 
0429 /* -------------------------------------------------------------------- */
0430 /*      Do we need to translate from ellipsoidal to geoid vertical      */
0431 /*      datum?                                                          */
0432 /* -------------------------------------------------------------------- */
0433     /*if( dstdefn->has_geoid_vgrids && z != NULL )
0434     {
0435         if( pj_apply_vgridshift( dstdefn, "sgeoidgrids",
0436                                  &(dstdefn->vgridlist_geoid),
0437                                  &(dstdefn->vgridlist_geoid_count),
0438                                  1, point_count, point_offset, x, y, z ) != 0 )
0439             return dstdefn->ctx->last_errno;
0440     }*/
0441 
0442 /* -------------------------------------------------------------------- */
0443 /*      But if they are staying lat long, adjust for the prime          */
0444 /*      meridian if there is one in effect.                             */
0445 /* -------------------------------------------------------------------- */
0446     if( dstdefn.from_greenwich != 0.0 )
0447     {
0448         for( std::size_t i = 0; i < point_count; i++ )
0449         {
0450             point_type & point = range::at(range, i);
0451 
0452             if( ! is_invalid_point(point) )
0453                 set<0>(point, get<0>(point) - dstdefn.from_greenwich);
0454         }
0455     }
0456 
0457 /* -------------------------------------------------------------------- */
0458 /*      Transform destination latlong to geocentric if required.        */
0459 /* -------------------------------------------------------------------- */
0460     if( BOOST_GEOMETRY_CONDITION(dstdefn.is_geocent) )
0461     {
0462         // Point should be cartesian 3D (ECEF)
0463         if ( BOOST_GEOMETRY_CONDITION(dimension < 3) )
0464             BOOST_THROW_EXCEPTION( projection_exception(error_geocentric) );
0465             //return PJD_ERR_GEOCENTRIC;
0466 
0467         // NOTE: In the original code the return value of the following
0468         // function is not checked
0469         range_wrapper<Range, false> rng(range);
0470         int err = pj_geodetic_to_geocentric( dstdefn.a_orig, dstdefn.es_orig,
0471                                              rng );
0472         if( err == -14 )
0473             result = false;
0474         else
0475             BOOST_THROW_EXCEPTION( projection_exception(err) );
0476 
0477         if( dstdefn.fr_meter != 1.0 )
0478         {
0479             for( std::size_t i = 0; i < point_count; i++ )
0480             {
0481                 point_type & point = range::at(range, i);
0482                 if( ! is_invalid_point(point) )
0483                 {
0484                     set<0>(point, get<0>(point) * dstdefn.fr_meter);
0485                     set<1>(point, get<1>(point) * dstdefn.fr_meter);
0486                 }
0487             }
0488         }
0489     }
0490 
0491 /* -------------------------------------------------------------------- */
0492 /*      Transform destination points to projection coordinates, if      */
0493 /*      desired.                                                        */
0494 /* -------------------------------------------------------------------- */
0495     else if( !dstdefn.is_latlong )
0496     {
0497 
0498         //if( dstdefn->fwd3d != NULL)
0499         //{
0500         //    for( i = 0; i < point_count; i++ )
0501         //    {
0502         //        XYZ projected_loc;
0503         //        LPZ geodetic_loc;
0504 
0505         //        geodetic_loc.u = x[point_offset*i];
0506         //        geodetic_loc.v = y[point_offset*i];
0507         //        geodetic_loc.w = z[point_offset*i];
0508 
0509         //        if (geodetic_loc.u == HUGE_VAL)
0510         //            continue;
0511 
0512         //        projected_loc = pj_fwd3d( geodetic_loc, dstdefn);
0513         //        if( dstdefn->ctx->last_errno != 0 )
0514         //        {
0515         //            if( (dstdefn->ctx->last_errno != 33 /*EDOM*/
0516         //                 && dstdefn->ctx->last_errno != 34 /*ERANGE*/ )
0517         //                && (dstdefn->ctx->last_errno > 0
0518         //                    || dstdefn->ctx->last_errno < -44 || point_count == 1
0519         //                    || transient_error[-dstdefn->ctx->last_errno] == 0 ) )
0520         //                return dstdefn->ctx->last_errno;
0521         //            else
0522         //            {
0523         //                projected_loc.u = HUGE_VAL;
0524         //                projected_loc.v = HUGE_VAL;
0525         //                projected_loc.w = HUGE_VAL;
0526         //            }
0527         //        }
0528 
0529         //        x[point_offset*i] = projected_loc.u;
0530         //        y[point_offset*i] = projected_loc.v;
0531         //        z[point_offset*i] = projected_loc.w;
0532         //    }
0533 
0534         //}
0535         //else
0536         {
0537             for(std::size_t i = 0; i < point_count; i++ )
0538             {
0539                 point_type & point = range::at(range, i);
0540 
0541                 if( is_invalid_point(point) )
0542                     continue;
0543 
0544                 try {
0545                     pj_fwd(dstprj, dstdefn, point, point);
0546                 } catch (projection_exception const& e) {
0547 
0548                     if( (e.code() != 33 /*EDOM*/
0549                          && e.code() != 34 /*ERANGE*/ )
0550                         && (e.code() > 0
0551                             || e.code() < -44 /*|| point_count == 1*/
0552                             || transient_error[-e.code()] == 0) ) {
0553                         BOOST_RETHROW
0554                     } else {
0555                         set_invalid_point(point);
0556                         result = false;
0557                         if (point_count == 1)
0558                             return result;
0559                     }
0560                 }
0561             }
0562         }
0563     }
0564 
0565 /* -------------------------------------------------------------------- */
0566 /*      If a wrapping center other than 0 is provided, rewrap around    */
0567 /*      the suggested center (for latlong coordinate systems only).     */
0568 /* -------------------------------------------------------------------- */
0569     else if( dstdefn.is_latlong && dstdefn.is_long_wrap_set )
0570     {
0571         for( std::size_t i = 0; i < point_count; i++ )
0572         {
0573             point_type & point = range::at(range, i);
0574             coord_t x = get_as_radian<0>(point);
0575 
0576             if( is_invalid_point(point) )
0577                 continue;
0578 
0579             // TODO - units-dependant constants could be used instead
0580             while( x < dstdefn.long_wrap_center - math::pi<coord_t>() )
0581                 x += math::two_pi<coord_t>();
0582             while( x > dstdefn.long_wrap_center + math::pi<coord_t>() )
0583                 x -= math::two_pi<coord_t>();
0584 
0585             set_from_radian<0>(point, x);
0586         }
0587     }
0588 
0589 /* -------------------------------------------------------------------- */
0590 /*      Transform Z from meters if needed.                              */
0591 /* -------------------------------------------------------------------- */
0592     if( dstdefn.vto_meter != 1.0 && dimension > 2 )
0593     {
0594         for( std::size_t i = 0; i < point_count; i++ )
0595         {
0596             point_type & point = geometry::range::at(range, i);
0597             set_z(point, get_z(point) * dstdefn.vfr_meter);
0598         }
0599     }
0600 
0601 /* -------------------------------------------------------------------- */
0602 /*      Transform normalized axes into unusual output coordinate axis   */
0603 /*      orientation if needed.                                          */
0604 /* -------------------------------------------------------------------- */
0605     // Ignored
0606 
0607     return result;
0608 }
0609 
0610 /************************************************************************/
0611 /*                     pj_geodetic_to_geocentric()                      */
0612 /************************************************************************/
0613 
0614 template <typename T, typename Range, bool AddZ>
0615 inline int pj_geodetic_to_geocentric( T const& a, T const& es,
0616                                       range_wrapper<Range, AddZ> & range_wrapper )
0617 
0618 {
0619     //typedef typename boost::range_iterator<Range>::type iterator;
0620     typedef typename boost::range_value<Range>::type point_type;
0621     //typedef typename coordinate_type<point_type>::type coord_t;
0622 
0623     Range & rng = range_wrapper.get_range();
0624     std::size_t point_count = boost::size(rng);
0625 
0626     int ret_errno = 0;
0627 
0628     T const b = (es == 0.0) ? a : a * sqrt(1-es);
0629 
0630     GeocentricInfo<T> gi;
0631     if( pj_Set_Geocentric_Parameters( gi, a, b ) != 0 )
0632     {
0633         return error_geocentric;
0634     }
0635 
0636     for( std::size_t i = 0 ; i < point_count ; ++i )
0637     {
0638         point_type & point = range::at(rng, i);
0639 
0640         if( is_invalid_point(point) )
0641             continue;
0642 
0643         T X = 0, Y = 0, Z = 0;
0644         if( pj_Convert_Geodetic_To_Geocentric( gi,
0645                                                get_as_radian<0>(point),
0646                                                get_as_radian<1>(point),
0647                                                range_wrapper.get_z(i), // Height
0648                                                X, Y, Z ) != 0 )
0649         {
0650             ret_errno = error_lat_or_lon_exceed_limit;
0651             set_invalid_point(point);
0652             /* but keep processing points! */
0653         }
0654         else
0655         {
0656             set<0>(point, X);
0657             set<1>(point, Y);
0658             range_wrapper.set_z(i, Z);
0659         }
0660     }
0661 
0662     return ret_errno;
0663 }
0664 
0665 /************************************************************************/
0666 /*                     pj_geodetic_to_geocentric()                      */
0667 /************************************************************************/
0668 
0669 template <typename T, typename Range, bool AddZ>
0670 inline int pj_geocentric_to_geodetic( T const& a, T const& es,
0671                                       range_wrapper<Range, AddZ> & range_wrapper )
0672 
0673 {
0674     //typedef typename boost::range_iterator<Range>::type iterator;
0675     typedef typename boost::range_value<Range>::type point_type;
0676     //typedef typename coordinate_type<point_type>::type coord_t;
0677 
0678     Range & rng = range_wrapper.get_range();
0679     std::size_t point_count = boost::size(rng);
0680 
0681     T const b = (es == 0.0) ? a : a * sqrt(1-es);
0682 
0683     GeocentricInfo<T> gi;
0684     if( pj_Set_Geocentric_Parameters( gi, a, b ) != 0 )
0685     {
0686         return error_geocentric;
0687     }
0688 
0689     for( std::size_t i = 0 ; i < point_count ; ++i )
0690     {
0691         point_type & point = range::at(rng, i);
0692 
0693         if( is_invalid_point(point) )
0694             continue;
0695 
0696         T Longitude = 0, Latitude = 0, Height = 0;
0697         pj_Convert_Geocentric_To_Geodetic( gi,
0698                                            get<0>(point),
0699                                            get<1>(point),
0700                                            range_wrapper.get_z(i), // z
0701                                            Longitude, Latitude, Height );
0702 
0703         set_from_radian<0>(point, Longitude);
0704         set_from_radian<1>(point, Latitude);
0705         range_wrapper.set_z(i, Height); // Height
0706     }
0707 
0708     return 0;
0709 }
0710 
0711 /************************************************************************/
0712 /*                         pj_compare_datums()                          */
0713 /*                                                                      */
0714 /*      Returns TRUE if the two datums are identical, otherwise         */
0715 /*      FALSE.                                                          */
0716 /************************************************************************/
0717 
0718 template <typename Par>
0719 inline bool pj_compare_datums( Par & srcdefn, Par & dstdefn )
0720 {
0721     if( srcdefn.datum_type != dstdefn.datum_type )
0722     {
0723         return false;
0724     }
0725     else if( srcdefn.a_orig != dstdefn.a_orig
0726              || math::abs(srcdefn.es_orig - dstdefn.es_orig) > 0.000000000050 )
0727     {
0728         /* the tolerance for es is to ensure that GRS80 and WGS84 are
0729            considered identical */
0730         return false;
0731     }
0732     else if( srcdefn.datum_type == datum_3param )
0733     {
0734         return (srcdefn.datum_params[0] == dstdefn.datum_params[0]
0735                 && srcdefn.datum_params[1] == dstdefn.datum_params[1]
0736                 && srcdefn.datum_params[2] == dstdefn.datum_params[2]);
0737     }
0738     else if( srcdefn.datum_type == datum_7param )
0739     {
0740         return (srcdefn.datum_params[0] == dstdefn.datum_params[0]
0741                 && srcdefn.datum_params[1] == dstdefn.datum_params[1]
0742                 && srcdefn.datum_params[2] == dstdefn.datum_params[2]
0743                 && srcdefn.datum_params[3] == dstdefn.datum_params[3]
0744                 && srcdefn.datum_params[4] == dstdefn.datum_params[4]
0745                 && srcdefn.datum_params[5] == dstdefn.datum_params[5]
0746                 && srcdefn.datum_params[6] == dstdefn.datum_params[6]);
0747     }
0748     else if( srcdefn.datum_type == datum_gridshift )
0749     {
0750         return srcdefn.nadgrids == dstdefn.nadgrids;
0751     }
0752     else
0753         return true;
0754 }
0755 
0756 /************************************************************************/
0757 /*                       pj_geocentic_to_wgs84()                        */
0758 /************************************************************************/
0759 
0760 template <typename Par, typename Range, bool AddZ>
0761 inline int pj_geocentric_to_wgs84( Par const& defn,
0762                                    range_wrapper<Range, AddZ> & range_wrapper )
0763 
0764 {
0765     typedef typename boost::range_value<Range>::type point_type;
0766     typedef typename coordinate_type<point_type>::type coord_t;
0767 
0768     Range & rng = range_wrapper.get_range();
0769     std::size_t point_count = boost::size(rng);
0770 
0771     if( defn.datum_type == datum_3param )
0772     {
0773         for(std::size_t i = 0; i < point_count; i++ )
0774         {
0775             point_type & point = range::at(rng, i);
0776 
0777             if( is_invalid_point(point) )
0778                 continue;
0779 
0780             set<0>(point,                   get<0>(point) + Dx_BF(defn));
0781             set<1>(point,                   get<1>(point) + Dy_BF(defn));
0782             range_wrapper.set_z(i, range_wrapper.get_z(i) + Dz_BF(defn));
0783         }
0784     }
0785     else if( defn.datum_type == datum_7param )
0786     {
0787         for(std::size_t i = 0; i < point_count; i++ )
0788         {
0789             point_type & point = range::at(rng, i);
0790 
0791             if( is_invalid_point(point) )
0792                 continue;
0793 
0794             coord_t x = get<0>(point);
0795             coord_t y = get<1>(point);
0796             coord_t z = range_wrapper.get_z(i);
0797 
0798             coord_t x_out, y_out, z_out;
0799 
0800             x_out = M_BF(defn)*(             x - Rz_BF(defn)*y + Ry_BF(defn)*z) + Dx_BF(defn);
0801             y_out = M_BF(defn)*( Rz_BF(defn)*x +             y - Rx_BF(defn)*z) + Dy_BF(defn);
0802             z_out = M_BF(defn)*(-Ry_BF(defn)*x + Rx_BF(defn)*y +             z) + Dz_BF(defn);
0803 
0804             set<0>(point, x_out);
0805             set<1>(point, y_out);
0806             range_wrapper.set_z(i, z_out);
0807         }
0808     }
0809 
0810     return 0;
0811 }
0812 
0813 /************************************************************************/
0814 /*                      pj_geocentic_from_wgs84()                       */
0815 /************************************************************************/
0816 
0817 template <typename Par, typename Range, bool AddZ>
0818 inline int pj_geocentric_from_wgs84( Par const& defn,
0819                                      range_wrapper<Range, AddZ> & range_wrapper )
0820 
0821 {
0822     typedef typename boost::range_value<Range>::type point_type;
0823     typedef typename coordinate_type<point_type>::type coord_t;
0824 
0825     Range & rng = range_wrapper.get_range();
0826     std::size_t point_count = boost::size(rng);
0827 
0828     if( defn.datum_type == datum_3param )
0829     {
0830         for(std::size_t i = 0; i < point_count; i++ )
0831         {
0832             point_type & point = range::at(rng, i);
0833 
0834             if( is_invalid_point(point) )
0835                 continue;
0836 
0837             set<0>(point,                   get<0>(point) - Dx_BF(defn));
0838             set<1>(point,                   get<1>(point) - Dy_BF(defn));
0839             range_wrapper.set_z(i, range_wrapper.get_z(i) - Dz_BF(defn));
0840         }
0841     }
0842     else if( defn.datum_type == datum_7param )
0843     {
0844         for(std::size_t i = 0; i < point_count; i++ )
0845         {
0846             point_type & point = range::at(rng, i);
0847 
0848             if( is_invalid_point(point) )
0849                 continue;
0850 
0851             coord_t x = get<0>(point);
0852             coord_t y = get<1>(point);
0853             coord_t z = range_wrapper.get_z(i);
0854 
0855             coord_t x_tmp = (x - Dx_BF(defn)) / M_BF(defn);
0856             coord_t y_tmp = (y - Dy_BF(defn)) / M_BF(defn);
0857             coord_t z_tmp = (z - Dz_BF(defn)) / M_BF(defn);
0858 
0859             x =              x_tmp + Rz_BF(defn)*y_tmp - Ry_BF(defn)*z_tmp;
0860             y = -Rz_BF(defn)*x_tmp +             y_tmp + Rx_BF(defn)*z_tmp;
0861             z =  Ry_BF(defn)*x_tmp - Rx_BF(defn)*y_tmp +             z_tmp;
0862 
0863             set<0>(point, x);
0864             set<1>(point, y);
0865             range_wrapper.set_z(i, z);
0866         }
0867     }
0868 
0869     return 0;
0870 }
0871 
0872 
0873 inline bool pj_datum_check_error(int err)
0874 {
0875     return err != 0 && (err > 0 || transient_error[-err] == 0);
0876 }
0877 
0878 /************************************************************************/
0879 /*                         pj_datum_transform()                         */
0880 /*                                                                      */
0881 /*      The input should be long/lat/z coordinates in radians in the    */
0882 /*      source datum, and the output should be long/lat/z               */
0883 /*      coordinates in radians in the destination datum.                */
0884 /************************************************************************/
0885 
0886 template <typename Par, typename Range, typename Grids>
0887 inline bool pj_datum_transform(Par const& srcdefn,
0888                                Par const& dstdefn,
0889                                Range & range,
0890                                Grids const& srcgrids,
0891                                Grids const& dstgrids)
0892 
0893 {
0894     typedef typename Par::type calc_t;
0895 
0896     // This has to be consistent with default spheroid and pj_ellps
0897     // TODO: Define in one place
0898     static const calc_t wgs84_a = 6378137.0;
0899     static const calc_t wgs84_b = 6356752.3142451793;
0900     static const calc_t wgs84_es = 1. - (wgs84_b * wgs84_b) / (wgs84_a * wgs84_a);
0901 
0902     bool result = true;
0903 
0904     calc_t      src_a, src_es, dst_a, dst_es;
0905 
0906 /* -------------------------------------------------------------------- */
0907 /*      We cannot do any meaningful datum transformation if either      */
0908 /*      the source or destination are of an unknown datum type          */
0909 /*      (ie. only a +ellps declaration, no +datum).  This is new        */
0910 /*      behavior for PROJ 4.6.0.                                        */
0911 /* -------------------------------------------------------------------- */
0912     if( srcdefn.datum_type == datum_unknown
0913         || dstdefn.datum_type == datum_unknown )
0914         return result;
0915 
0916 /* -------------------------------------------------------------------- */
0917 /*      Short cut if the datums are identical.                          */
0918 /* -------------------------------------------------------------------- */
0919     if( pj_compare_datums( srcdefn, dstdefn ) )
0920         return result;
0921 
0922     src_a = srcdefn.a_orig;
0923     src_es = srcdefn.es_orig;
0924 
0925     dst_a = dstdefn.a_orig;
0926     dst_es = dstdefn.es_orig;
0927 
0928 /* -------------------------------------------------------------------- */
0929 /*      Create a temporary Z array if one is not provided.              */
0930 /* -------------------------------------------------------------------- */
0931 
0932     range_wrapper<Range> z_range(range);
0933 
0934 /* -------------------------------------------------------------------- */
0935 /*      If this datum requires grid shifts, then apply it to geodetic   */
0936 /*      coordinates.                                                    */
0937 /* -------------------------------------------------------------------- */
0938     if( srcdefn.datum_type == datum_gridshift )
0939     {
0940         try {
0941             pj_apply_gridshift_2<false>( srcdefn, range, srcgrids );
0942         } catch (projection_exception const& e) {
0943             if (pj_datum_check_error(e.code())) {
0944                 BOOST_RETHROW
0945             }
0946         }
0947 
0948         src_a = wgs84_a;
0949         src_es = wgs84_es;
0950     }
0951 
0952     if( dstdefn.datum_type == datum_gridshift )
0953     {
0954         dst_a = wgs84_a;
0955         dst_es = wgs84_es;
0956     }
0957 
0958 /* ==================================================================== */
0959 /*      Do we need to go through geocentric coordinates?                */
0960 /* ==================================================================== */
0961     if( src_es != dst_es || src_a != dst_a
0962         || srcdefn.datum_type == datum_3param
0963         || srcdefn.datum_type == datum_7param
0964         || dstdefn.datum_type == datum_3param
0965         || dstdefn.datum_type == datum_7param)
0966     {
0967 /* -------------------------------------------------------------------- */
0968 /*      Convert to geocentric coordinates.                              */
0969 /* -------------------------------------------------------------------- */
0970         int err = pj_geodetic_to_geocentric( src_a, src_es, z_range );
0971         if (pj_datum_check_error(err))
0972             BOOST_THROW_EXCEPTION( projection_exception(err) );
0973         else if (err != 0)
0974             result = false;
0975 
0976 /* -------------------------------------------------------------------- */
0977 /*      Convert between datums.                                         */
0978 /* -------------------------------------------------------------------- */
0979         if( srcdefn.datum_type == datum_3param
0980             || srcdefn.datum_type == datum_7param )
0981         {
0982             try {
0983                 pj_geocentric_to_wgs84( srcdefn, z_range );
0984             } catch (projection_exception const& e) {
0985                 if (pj_datum_check_error(e.code())) {
0986                     BOOST_RETHROW
0987                 }
0988             }
0989         }
0990 
0991         if( dstdefn.datum_type == datum_3param
0992             || dstdefn.datum_type == datum_7param )
0993         {
0994             try {
0995                 pj_geocentric_from_wgs84( dstdefn, z_range );
0996             } catch (projection_exception const& e) {
0997                 if (pj_datum_check_error(e.code())) {
0998                     BOOST_RETHROW
0999                 }
1000             }
1001         }
1002 
1003 /* -------------------------------------------------------------------- */
1004 /*      Convert back to geodetic coordinates.                           */
1005 /* -------------------------------------------------------------------- */
1006         err = pj_geocentric_to_geodetic( dst_a, dst_es, z_range );
1007         if (pj_datum_check_error(err))
1008             BOOST_THROW_EXCEPTION( projection_exception(err) );
1009         else if (err != 0)
1010             result = false;
1011     }
1012 
1013 /* -------------------------------------------------------------------- */
1014 /*      Apply grid shift to destination if required.                    */
1015 /* -------------------------------------------------------------------- */
1016     if( dstdefn.datum_type == datum_gridshift )
1017     {
1018         try {
1019             pj_apply_gridshift_2<true>( dstdefn, range, dstgrids );
1020         } catch (projection_exception const& e) {
1021             if (pj_datum_check_error(e.code()))
1022                 BOOST_RETHROW
1023         }
1024     }
1025 
1026     return result;
1027 }
1028 
1029 } // namespace detail
1030 
1031 }}} // namespace boost::geometry::projections
1032 
1033 #endif // BOOST_GEOMETRY_SRS_PROJECTIONS_IMPL_PJ_TRANSFORM_HPP
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!