projections/impl/pj_apply_gridshift.hpp

0001 // Boost.Geometry
0002 // This file is manually converted from PROJ4
0003
0004 // This file was modified by Oracle on 2018, 2019.
0005 // Modifications copyright (c) 2018-2019, Oracle and/or its affiliates.
0006 // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
0007
0008 // Use, modification and distribution is subject to the Boost Software License,
0009 // Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
0010 // http://www.boost.org/LICENSE_1_0.txt)
0011
0012 // This file is converted from PROJ4, http://trac.osgeo.org/proj
0013 // PROJ4 is originally written by Gerald Evenden (then of the USGS)
0014 // PROJ4 is maintained by Frank Warmerdam
0015 // This file was converted to Geometry Library by Adam Wulkiewicz
0016
0017 // Original copyright notice:
0018 // Author:   Frank Warmerdam, warmerdam@pobox.com
0019
0020 // Copyright (c) 2000, Frank Warmerdam
0021
0022 // Permission is hereby granted, free of charge, to any person obtaining a
0023 // copy of this software and associated documentation files (the "Software"),
0024 // to deal in the Software without restriction, including without limitation
0025 // the rights to use, copy, modify, merge, publish, distribute, sublicense,
0026 // and/or sell copies of the Software, and to permit persons to whom the
0027 // Software is furnished to do so, subject to the following conditions:
0028
0029 // The above copyright notice and this permission notice shall be included
0030 // in all copies or substantial portions of the Software.
0031
0032 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
0033 // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
0034 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
0035 // THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
0036 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
0037 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
0038 // DEALINGS IN THE SOFTWARE.
0039
0040 #ifndef BOOST_GEOMETRY_SRS_PROJECTIONS_IMPL_PJ_APPLY_GRIDSHIFT_HPP
0041 #define BOOST_GEOMETRY_SRS_PROJECTIONS_IMPL_PJ_APPLY_GRIDSHIFT_HPP
0042
0043
0044 #include <boost/geometry/core/assert.hpp>
0045 #include <boost/geometry/core/radian_access.hpp>
0046
0047 #include <boost/geometry/srs/projections/impl/adjlon.hpp>
0048 #include <boost/geometry/srs/projections/impl/function_overloads.hpp>
0049 #include <boost/geometry/srs/projections/impl/pj_gridlist.hpp>
0050
0051 #include <boost/geometry/util/range.hpp>
0052
0053
0054 namespace boost { namespace geometry { namespace projections
0055 {
0056
0057 namespace detail
0058 {
0059
0060 // Originally implemented in nad_intr.c
0061 template <typename CalcT>
0062 inline void nad_intr(CalcT in_lon, CalcT in_lat,
0063                      CalcT & out_lon, CalcT & out_lat,
0064                      pj_ctable const& ct)
0065 {
0066     pj_ctable::lp_t frct;
0067     pj_ctable::ilp_t indx;
0068     boost::int32_t in;
0069
0070     indx.lam = int_floor(in_lon /= ct.del.lam);
0071     indx.phi = int_floor(in_lat /= ct.del.phi);
0072     frct.lam = in_lon - indx.lam;
0073     frct.phi = in_lat - indx.phi;
0074     // TODO: implement differently
0075     out_lon = out_lat = HUGE_VAL;
0076     if (indx.lam < 0) {
0077         if (indx.lam == -1 && frct.lam > 0.99999999999) {
0078             ++indx.lam;
0079             frct.lam = 0.;
0080         } else
0081             return;
0082     } else if ((in = indx.lam + 1) >= ct.lim.lam) {
0083         if (in == ct.lim.lam && frct.lam < 1e-11) {
0084             --indx.lam;
0085             frct.lam = 1.;
0086         } else
0087             return;
0088     }
0089     if (indx.phi < 0) {
0090         if (indx.phi == -1 && frct.phi > 0.99999999999) {
0091             ++indx.phi;
0092             frct.phi = 0.;
0093         } else
0094             return;
0095     } else if ((in = indx.phi + 1) >= ct.lim.phi) {
0096         if (in == ct.lim.phi && frct.phi < 1e-11) {
0097             --indx.phi;
0098             frct.phi = 1.;
0099         } else
0100             return;
0101     }
0102     boost::int32_t index = indx.phi * ct.lim.lam + indx.lam;
0103     pj_ctable::flp_t const& f00 = ct.cvs[index++];
0104     pj_ctable::flp_t const& f10 = ct.cvs[index];
0105     index += ct.lim.lam;
0106     pj_ctable::flp_t const& f11 = ct.cvs[index--];
0107     pj_ctable::flp_t const& f01 = ct.cvs[index];
0108     CalcT m00, m10, m01, m11;
0109     m11 = m10 = frct.lam;
0110     m00 = m01 = 1. - frct.lam;
0111     m11 *= frct.phi;
0112     m01 *= frct.phi;
0113     frct.phi = 1. - frct.phi;
0114     m00 *= frct.phi;
0115     m10 *= frct.phi;
0116     out_lon = m00 * f00.lam + m10 * f10.lam +
0117               m01 * f01.lam + m11 * f11.lam;
0118     out_lat = m00 * f00.phi + m10 * f10.phi +
0119               m01 * f01.phi + m11 * f11.phi;
0120 }
0121
0122 // Originally implemented in nad_cvt.c
0123 template <bool Inverse, typename CalcT>
0124 inline void nad_cvt(CalcT const& in_lon, CalcT const& in_lat,
0125                     CalcT & out_lon, CalcT & out_lat,
0126                     pj_gi const& gi)
0127 {
0128     static const int max_iterations = 10;
0129     static const CalcT tol = 1e-12;
0130     static const CalcT toltol = tol * tol;
0131     static const CalcT pi = math::pi<CalcT>();
0132
0133     // horizontal grid expected
0134     BOOST_GEOMETRY_ASSERT_MSG(gi.format != pj_gi::gtx,
0135         "Vertical grid cannot be used in horizontal shift.");
0136
0137     pj_ctable const& ct = gi.ct;
0138
0139     // TODO: implement differently
0140     if (in_lon == HUGE_VAL)
0141     {
0142         out_lon = HUGE_VAL;
0143         out_lat = HUGE_VAL;
0144         return;
0145     }
0146
0147     // normalize input to ll origin
0148     pj_ctable::lp_t tb;
0149     tb.lam = in_lon - ct.ll.lam;
0150     tb.phi = in_lat - ct.ll.phi;
0151     tb.lam = adjlon (tb.lam - pi) + pi;
0152
0153     pj_ctable::lp_t t;
0154     nad_intr(tb.lam, tb.phi, t.lam, t.phi, ct);
0155     if (t.lam == HUGE_VAL)
0156     {
0157         out_lon = HUGE_VAL;
0158         out_lat = HUGE_VAL;
0159         return;
0160     }
0161
0162     if (! Inverse)
0163     {
0164         out_lon = in_lon - t.lam;
0165         out_lat = in_lat - t.phi;
0166         return;
0167     }
0168
0169     t.lam = tb.lam + t.lam;
0170     t.phi = tb.phi - t.phi;
0171
0172     int i = max_iterations;
0173     pj_ctable::lp_t del, dif;
0174     do
0175     {
0176         nad_intr(t.lam, t.phi, del.lam, del.phi, ct);
0177
0178         // This case used to return failure, but I have
0179         // changed it to return the first order approximation
0180         // of the inverse shift.  This avoids cases where the
0181         // grid shift *into* this grid came from another grid.
0182         // While we aren't returning optimally correct results
0183         // I feel a close result in this case is better than
0184         // no result.  NFW
0185         // To demonstrate use -112.5839956 49.4914451 against
0186         // the NTv2 grid shift file from Canada.
0187         if (del.lam == HUGE_VAL)
0188         {
0189             // Inverse grid shift iteration failed, presumably at grid edge. Using first approximation.
0190             break;
0191         }
0192
0193         dif.lam = t.lam - del.lam - tb.lam;
0194         dif.phi = t.phi + del.phi - tb.phi;
0195         t.lam -= dif.lam;
0196         t.phi -= dif.phi;
0197
0198     }
0199     while (--i && (dif.lam*dif.lam + dif.phi*dif.phi > toltol)); // prob. slightly faster than hypot()
0200
0201     if (i==0)
0202     {
0203         // Inverse grid shift iterator failed to converge.
0204         out_lon = HUGE_VAL;
0205         out_lat = HUGE_VAL;
0206         return;
0207     }
0208
0209     out_lon = adjlon (t.lam + ct.ll.lam);
0210     out_lat = t.phi + ct.ll.phi;
0211 }
0212
0213
0214 /************************************************************************/
0215 /*                             find_grid()                              */
0216 /*                                                                      */
0217 /*    Determine which grid is the correct given an input coordinate.    */
0218 /************************************************************************/
0219
0220 // Originally find_ctable()
0221 // here divided into grid_disjoint(), find_grid() and load_grid()
0222
0223 template <typename T>
0224 inline bool grid_disjoint(T const& lam, T const& phi,
0225                           pj_ctable const& ct)
0226 {
0227     double epsilon = (fabs(ct.del.phi)+fabs(ct.del.lam))/10000.0;
0228     return ct.ll.phi - epsilon > phi
0229         || ct.ll.lam - epsilon > lam
0230         || (ct.ll.phi + (ct.lim.phi-1) * ct.del.phi + epsilon < phi)
0231         || (ct.ll.lam + (ct.lim.lam-1) * ct.del.lam + epsilon < lam);
0232 }
0233
0234 template <typename T>
0235 inline pj_gi * find_grid(T const& lam,
0236                          T const& phi,
0237                          std::vector<pj_gi>::iterator first,
0238                          std::vector<pj_gi>::iterator last)
0239 {
0240     pj_gi * gip = NULL;
0241
0242     for( ; first != last ; ++first )
0243     {
0244         // skip tables that don't match our point at all.
0245         if (! grid_disjoint(lam, phi, first->ct))
0246         {
0247             // skip vertical grids
0248             if (first->format != pj_gi::gtx)
0249             {
0250                 gip = boost::addressof(*first);
0251                 break;
0252             }
0253         }
0254     }
0255
0256     // If we didn't find a child then nothing more to do
0257     if( gip == NULL )
0258         return gip;
0259
0260     // Otherwise use the child, first checking it's children
0261     pj_gi * child = find_grid(lam, phi, first->children.begin(), first->children.end());
0262     if (child != NULL)
0263         gip = child;
0264
0265     return gip;
0266 }
0267
0268 template <typename T>
0269 inline pj_gi * find_grid(T const& lam,
0270                          T const& phi,
0271                          pj_gridinfo & grids,
0272                          std::vector<std::size_t> const& gridindexes)
0273 {
0274     pj_gi * gip = NULL;
0275
0276     // keep trying till we find a table that works
0277     for (std::size_t i = 0 ; i < gridindexes.size() ; ++i)
0278     {
0279         pj_gi & gi = grids[gridindexes[i]];
0280
0281         // skip tables that don't match our point at all.
0282         if (! grid_disjoint(lam, phi, gi.ct))
0283         {
0284             // skip vertical grids
0285             if (gi.format != pj_gi::gtx)
0286             {
0287                 gip = boost::addressof(gi);
0288                 break;
0289             }
0290         }
0291     }
0292
0293     if (gip == NULL)
0294         return gip;
0295
0296     // If we have child nodes, check to see if any of them apply.
0297     pj_gi * child = find_grid(lam, phi, gip->children.begin(), gip->children.end());
0298     if (child != NULL)
0299         gip = child;
0300
0301     // if we get this far we have found a suitable grid
0302     return gip;
0303 }
0304
0305
0306 template <typename StreamPolicy>
0307 inline bool load_grid(StreamPolicy const& stream_policy, pj_gi_load & gi)
0308 {
0309     // load the grid shift info if we don't have it.
0310     if (gi.ct.cvs.empty())
0311     {
0312         typename StreamPolicy::stream_type is;
0313         stream_policy.open(is, gi.gridname);
0314
0315         if (! pj_gridinfo_load(is, gi))
0316         {
0317             //pj_ctx_set_errno( ctx, PJD_ERR_FAILED_TO_LOAD_GRID );
0318             return false;
0319         }
0320     }
0321
0322     return true;
0323 }
0324
0325
0326 /************************************************************************/
0327 /*                        pj_apply_gridshift_3()                        */
0328 /*                                                                      */
0329 /*      This is the real workhorse, given a gridlist.                   */
0330 /************************************************************************/
0331
0332 // Generic stream policy and standard grids
0333 template <bool Inverse, typename CalcT, typename StreamPolicy, typename Range, typename Grids>
0334 inline bool pj_apply_gridshift_3(StreamPolicy const& stream_policy,
0335                                  Range & range,
0336                                  Grids & grids,
0337                                  std::vector<std::size_t> const& gridindexes,
0338                                  grids_tag)
0339 {
0340     typedef typename boost::range_size<Range>::type size_type;
0341
0342     // If the grids are empty the indexes are as well
0343     if (gridindexes.empty())
0344     {
0345         //pj_ctx_set_errno(ctx, PJD_ERR_FAILED_TO_LOAD_GRID);
0346         //return PJD_ERR_FAILED_TO_LOAD_GRID;
0347         return false;
0348     }
0349
0350     size_type point_count = boost::size(range);
0351
0352     for (size_type i = 0 ; i < point_count ; ++i)
0353     {
0354         typename boost::range_reference<Range>::type
0355             point = range::at(range, i);
0356
0357         CalcT in_lon = geometry::get_as_radian<0>(point);
0358         CalcT in_lat = geometry::get_as_radian<1>(point);
0359
0360         pj_gi * gip = find_grid(in_lon, in_lat, grids.gridinfo, gridindexes);
0361
0362         if ( gip != NULL )
0363         {
0364             // load the grid shift info if we don't have it.
0365             if (! gip->ct.cvs.empty() || load_grid(stream_policy, *gip))
0366             {
0367                 // TODO: use set_invalid_point() or similar mechanism
0368                 CalcT out_lon = HUGE_VAL;
0369                 CalcT out_lat = HUGE_VAL;
0370
0371                 nad_cvt<Inverse>(in_lon, in_lat, out_lon, out_lat, *gip);
0372
0373                 // TODO: check differently
0374                 if ( out_lon != HUGE_VAL )
0375                 {
0376                     geometry::set_from_radian<0>(point, out_lon);
0377                     geometry::set_from_radian<1>(point, out_lat);
0378                 }
0379             }
0380         }
0381     }
0382
0383     return true;
0384 }
0385
0386 // Generic stream policy and shared grids
0387 template <bool Inverse, typename CalcT, typename StreamPolicy, typename Range, typename SharedGrids>
0388 inline bool pj_apply_gridshift_3(StreamPolicy const& stream_policy,
0389                                  Range & range,
0390                                  SharedGrids & grids,
0391                                  std::vector<std::size_t> const& gridindexes,
0392                                  shared_grids_tag)
0393 {
0394     typedef typename boost::range_size<Range>::type size_type;
0395
0396     // If the grids are empty the indexes are as well
0397     if (gridindexes.empty())
0398     {
0399         //pj_ctx_set_errno(ctx, PJD_ERR_FAILED_TO_LOAD_GRID);
0400         //return PJD_ERR_FAILED_TO_LOAD_GRID;
0401         return false;
0402     }
0403
0404     size_type point_count = boost::size(range);
0405
0406     // local storage
0407     pj_gi_load local_gi;
0408
0409     for (size_type i = 0 ; i < point_count ; )
0410     {
0411         bool load_needed = false;
0412
0413         CalcT in_lon = 0;
0414         CalcT in_lat = 0;
0415
0416         {
0417             typename SharedGrids::read_locked lck_grids(grids);
0418
0419             for ( ; i < point_count ; ++i )
0420             {
0421                 typename boost::range_reference<Range>::type
0422                     point = range::at(range, i);
0423
0424                 in_lon = geometry::get_as_radian<0>(point);
0425                 in_lat = geometry::get_as_radian<1>(point);
0426
0427                 pj_gi * gip = find_grid(in_lon, in_lat, lck_grids.gridinfo, gridindexes);
0428
0429                 if (gip == NULL)
0430                 {
0431                     // do nothing
0432                 }
0433                 else if (! gip->ct.cvs.empty())
0434                 {
0435                     // TODO: use set_invalid_point() or similar mechanism
0436                     CalcT out_lon = HUGE_VAL;
0437                     CalcT out_lat = HUGE_VAL;
0438
0439                     nad_cvt<Inverse>(in_lon, in_lat, out_lon, out_lat, *gip);
0440
0441                     // TODO: check differently
0442                     if (out_lon != HUGE_VAL)
0443                     {
0444                         geometry::set_from_radian<0>(point, out_lon);
0445                         geometry::set_from_radian<1>(point, out_lat);
0446                     }
0447                 }
0448                 else
0449                 {
0450                     // loading is needed
0451                     local_gi = *gip;
0452                     load_needed = true;
0453                     break;
0454                 }
0455             }
0456         }
0457
0458         if (load_needed)
0459         {
0460             if (load_grid(stream_policy, local_gi))
0461             {
0462                 typename SharedGrids::write_locked lck_grids(grids);
0463
0464                 // check again in case other thread already loaded the grid.
0465                 pj_gi * gip = find_grid(in_lon, in_lat, lck_grids.gridinfo, gridindexes);
0466
0467                 if (gip != NULL && gip->ct.cvs.empty())
0468                 {
0469                     // swap loaded local storage with empty grid
0470                     local_gi.swap(*gip);
0471                 }
0472             }
0473             else
0474             {
0475                 ++i;
0476             }
0477         }
0478     }
0479
0480     return true;
0481 }
0482
0483
0484 /************************************************************************/
0485 /*                        pj_apply_gridshift_2()                        */
0486 /*                                                                      */
0487 /*      This implementation uses the gridlist from a coordinate         */
0488 /*      system definition.  If the gridlist has not yet been            */
0489 /*      populated in the coordinate system definition we set it up      */
0490 /*      now.                                                            */
0491 /************************************************************************/
0492
0493 template <bool Inverse, typename Par, typename Range, typename ProjGrids>
0494 inline bool pj_apply_gridshift_2(Par const& /*defn*/, Range & range, ProjGrids const& proj_grids)
0495 {
0496     /*if( defn->catalog_name != NULL )
0497         return pj_gc_apply_gridshift( defn, inverse, point_count, point_offset,
0498                                       x, y, z );*/
0499
0500     /*std::vector<std::size_t> gridindexes;
0501     pj_gridlist_from_nadgrids(pj_get_param_s(defn.params, "nadgrids"),
0502                               grids.storage_ptr->stream_policy,
0503                               grids.storage_ptr->grids,
0504                               gridindexes);*/
0505
0506     // At this point the grids should be initialized
0507     if (proj_grids.hindexes.empty())
0508         return false;
0509
0510     return pj_apply_gridshift_3
0511             <
0512                 Inverse, typename Par::type
0513             >(proj_grids.grids_storage().stream_policy,
0514               range,
0515               proj_grids.grids_storage().hgrids,
0516               proj_grids.hindexes,
0517               typename ProjGrids::grids_storage_type::grids_type::tag());
0518 }
0519
0520 template <bool Inverse, typename Par, typename Range>
0521 inline bool pj_apply_gridshift_2(Par const& , Range & , srs::detail::empty_projection_grids const& )
0522 {
0523     return false;
0524 }
0525
0526
0527 } // namespace detail
0528
0529 }}} // namespace boost::geometry::projections
0530
0531 #endif // BOOST_GEOMETRY_SRS_PROJECTIONS_IMPL_PJ_APPLY_GRIDSHIFT_HPP