boost/geometry/srs/projections/impl/pj_apply_gridshift.hpp

   1 // Boost.Geometry
   2 // This file is manually converted from PROJ4
   3
   4 // This file was modified by Oracle on 2018, 2019.
   5 // Modifications copyright (c) 2018-2019, Oracle and/or its affiliates.
   6 // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
   7
   8 // Use, modification and distribution is subject to the Boost Software License,
   9 // Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
  10 // http://www.boost.org/LICENSE_1_0.txt)
  11
  12 // This file is converted from PROJ4, http://trac.osgeo.org/proj
  13 // PROJ4 is originally written by Gerald Evenden (then of the USGS)
  14 // PROJ4 is maintained by Frank Warmerdam
  15 // This file was converted to Geometry Library by Adam Wulkiewicz
  16
  17 // Original copyright notice:
  18 // Author:   Frank Warmerdam, warmerdam@pobox.com
  19
  20 // Copyright (c) 2000, Frank Warmerdam
  21
  22 // Permission is hereby granted, free of charge, to any person obtaining a
  23 // copy of this software and associated documentation files (the "Software"),
  24 // to deal in the Software without restriction, including without limitation
  25 // the rights to use, copy, modify, merge, publish, distribute, sublicense,
  26 // and/or sell copies of the Software, and to permit persons to whom the
  27 // Software is furnished to do so, subject to the following conditions:
  28
  29 // The above copyright notice and this permission notice shall be included
  30 // in all copies or substantial portions of the Software.
  31
  32 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
  33 // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
  34 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
  35 // THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
  36 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
  37 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
  38 // DEALINGS IN THE SOFTWARE.
  39
  40 #ifndef BOOST_GEOMETRY_SRS_PROJECTIONS_IMPL_PJ_APPLY_GRIDSHIFT_HPP
  41 #define BOOST_GEOMETRY_SRS_PROJECTIONS_IMPL_PJ_APPLY_GRIDSHIFT_HPP
  42
  43
  44 #include <boost/geometry/core/assert.hpp>
  45 #include <boost/geometry/core/radian_access.hpp>
  46
  47 #include <boost/geometry/srs/projections/impl/adjlon.hpp>
  48 #include <boost/geometry/srs/projections/impl/function_overloads.hpp>
  49 #include <boost/geometry/srs/projections/impl/pj_gridlist.hpp>
  50
  51 #include <boost/geometry/util/range.hpp>
  52
  53
  54 namespace boost { namespace geometry { namespace projections
  55 {
  56
  57 namespace detail
  58 {
  59
  60 // Originally implemented in nad_intr.c
  61 template <typename CalcT>
  62 inline void nad_intr(CalcT in_lon, CalcT in_lat,
  63                      CalcT & out_lon, CalcT & out_lat,
  64                      pj_ctable const& ct)
  65 {
  66         pj_ctable::lp_t frct;
  67         pj_ctable::ilp_t indx;
  68         boost::int32_t in;
  69
  70         indx.lam = int_floor(in_lon /= ct.del.lam);
  71         indx.phi = int_floor(in_lat /= ct.del.phi);
  72         frct.lam = in_lon - indx.lam;
  73         frct.phi = in_lat - indx.phi;
  74     // TODO: implement differently
  75         out_lon = out_lat = HUGE_VAL;
  76         if (indx.lam < 0) {
  77                 if (indx.lam == -1 && frct.lam > 0.99999999999) {
  78                         ++indx.lam;
  79                         frct.lam = 0.;
  80                 } else
  81                         return;
  82         } else if ((in = indx.lam + 1) >= ct.lim.lam) {
  83                 if (in == ct.lim.lam && frct.lam < 1e-11) {
  84                         --indx.lam;
  85                         frct.lam = 1.;
  86                 } else
  87                         return;
  88         }
  89         if (indx.phi < 0) {
  90                 if (indx.phi == -1 && frct.phi > 0.99999999999) {
  91                         ++indx.phi;
  92                         frct.phi = 0.;
  93                 } else
  94                         return;
  95         } else if ((in = indx.phi + 1) >= ct.lim.phi) {
  96                 if (in == ct.lim.phi && frct.phi < 1e-11) {
  97                         --indx.phi;
  98                         frct.phi = 1.;
  99                 } else
 100                         return;
 101         }
 102         boost::int32_t index = indx.phi * ct.lim.lam + indx.lam;
 103         pj_ctable::flp_t const& f00 = ct.cvs[index++];
 104         pj_ctable::flp_t const& f10 = ct.cvs[index];
 105         index += ct.lim.lam;
 106         pj_ctable::flp_t const& f11 = ct.cvs[index--];
 107         pj_ctable::flp_t const& f01 = ct.cvs[index];
 108     CalcT m00, m10, m01, m11;
 109         m11 = m10 = frct.lam;
 110         m00 = m01 = 1. - frct.lam;
 111         m11 *= frct.phi;
 112         m01 *= frct.phi;
 113         frct.phi = 1. - frct.phi;
 114         m00 *= frct.phi;
 115         m10 *= frct.phi;
 116         out_lon = m00 * f00.lam + m10 * f10.lam +
 117                           m01 * f01.lam + m11 * f11.lam;
 118         out_lat = m00 * f00.phi + m10 * f10.phi +
 119                           m01 * f01.phi + m11 * f11.phi;
 120 }
 121
 122 // Originally implemented in nad_cvt.c
 123 template <bool Inverse, typename CalcT>
 124 inline void nad_cvt(CalcT const& in_lon, CalcT const& in_lat,
 125                     CalcT & out_lon, CalcT & out_lat,
 126                     pj_gi const& gi)
 127 {
 128     static const int max_iterations = 10;
 129     static const CalcT tol = 1e-12;
 130     static const CalcT toltol = tol * tol;
 131     static const CalcT pi = math::pi<CalcT>();
 132
 133     // horizontal grid expected
 134     BOOST_GEOMETRY_ASSERT_MSG(gi.format != pj_gi::gtx,
 135         "Vertical grid cannot be used in horizontal shift.");
 136
 137     pj_ctable const& ct = gi.ct;
 138
 139     // TODO: implement differently
 140     if (in_lon == HUGE_VAL)
 141     {
 142         out_lon = HUGE_VAL;
 143         out_lat = HUGE_VAL;
 144         return;
 145     }
 146
 147     // normalize input to ll origin
 148     pj_ctable::lp_t tb;
 149     tb.lam = in_lon - ct.ll.lam;
 150     tb.phi = in_lat - ct.ll.phi;
 151     tb.lam = adjlon (tb.lam - pi) + pi;
 152
 153     pj_ctable::lp_t t;
 154     nad_intr(tb.lam, tb.phi, t.lam, t.phi, ct);
 155     if (t.lam == HUGE_VAL)
 156     {
 157         out_lon = HUGE_VAL;
 158         out_lat = HUGE_VAL;
 159         return;
 160     }
 161
 162     if (! Inverse)
 163     {
 164         out_lon = in_lon - t.lam;
 165         out_lat = in_lat - t.phi;
 166         return;
 167     }
 168
 169     t.lam = tb.lam + t.lam;
 170     t.phi = tb.phi - t.phi;
 171
 172     int i = max_iterations;
 173     pj_ctable::lp_t del, dif;
 174     do
 175     {
 176         nad_intr(t.lam, t.phi, del.lam, del.phi, ct);
 177
 178         // This case used to return failure, but I have
 179         // changed it to return the first order approximation
 180         // of the inverse shift.  This avoids cases where the
 181         // grid shift *into* this grid came from another grid.
 182         // While we aren't returning optimally correct results
 183         // I feel a close result in this case is better than
 184         // no result.  NFW
 185         // To demonstrate use -112.5839956 49.4914451 against
 186         // the NTv2 grid shift file from Canada.
 187         if (del.lam == HUGE_VAL)
 188         {
 189             // Inverse grid shift iteration failed, presumably at grid edge. Using first approximation.
 190             break;
 191         }
 192
 193         dif.lam = t.lam - del.lam - tb.lam;
 194         dif.phi = t.phi + del.phi - tb.phi;
 195         t.lam -= dif.lam;
 196         t.phi -= dif.phi;
 197
 198     }
 199     while (--i && (dif.lam*dif.lam + dif.phi*dif.phi > toltol)); // prob. slightly faster than hypot()
 200
 201     if (i==0)
 202     {
 203         // Inverse grid shift iterator failed to converge.
 204         out_lon = HUGE_VAL;
 205         out_lat = HUGE_VAL;
 206         return;
 207     }
 208
 209     out_lon = adjlon (t.lam + ct.ll.lam);
 210     out_lat = t.phi + ct.ll.phi;
 211 }
 212
 213
 214 /************************************************************************/
 215 /*                             find_grid()                              */
 216 /*                                                                      */
 217 /*    Determine which grid is the correct given an input coordinate.    */
 218 /************************************************************************/
 219
 220 // Originally find_ctable()
 221 // here divided into grid_disjoint(), find_grid() and load_grid()
 222
 223 template <typename T>
 224 inline bool grid_disjoint(T const& lam, T const& phi,
 225                           pj_ctable const& ct)
 226 {
 227     double epsilon = (fabs(ct.del.phi)+fabs(ct.del.lam))/10000.0;
 228     return ct.ll.phi - epsilon > phi
 229         || ct.ll.lam - epsilon > lam
 230         || (ct.ll.phi + (ct.lim.phi-1) * ct.del.phi + epsilon < phi)
 231         || (ct.ll.lam + (ct.lim.lam-1) * ct.del.lam + epsilon < lam);
 232 }
 233
 234 template <typename T>
 235 inline pj_gi * find_grid(T const& lam,
 236                          T const& phi,
 237                          std::vector<pj_gi>::iterator first,
 238                          std::vector<pj_gi>::iterator last)
 239 {
 240     pj_gi * gip = NULL;
 241
 242     for( ; first != last ; ++first )
 243     {
 244         // skip tables that don't match our point at all.
 245         if (! grid_disjoint(lam, phi, first->ct))
 246         {
 247             // skip vertical grids
 248             if (first->format != pj_gi::gtx)
 249             {
 250                 gip = boost::addressof(*first);
 251                 break;
 252             }
 253         }
 254     }
 255
 256     // If we didn't find a child then nothing more to do
 257     if( gip == NULL )
 258         return gip;
 259
 260     // Otherwise use the child, first checking it's children
 261     pj_gi * child = find_grid(lam, phi, first->children.begin(), first->children.end());
 262     if (child != NULL)
 263         gip = child;
 264
 265     return gip;
 266 }
 267
 268 template <typename T>
 269 inline pj_gi * find_grid(T const& lam,
 270                          T const& phi,
 271                          pj_gridinfo & grids,
 272                          std::vector<std::size_t> const& gridindexes)
 273 {
 274     pj_gi * gip = NULL;
 275
 276     // keep trying till we find a table that works
 277     for (std::size_t i = 0 ; i < gridindexes.size() ; ++i)
 278     {
 279         pj_gi & gi = grids[gridindexes[i]];
 280
 281         // skip tables that don't match our point at all.
 282         if (! grid_disjoint(lam, phi, gi.ct))
 283         {
 284             // skip vertical grids
 285             if (gi.format != pj_gi::gtx)
 286             {
 287                 gip = boost::addressof(gi);
 288                 break;
 289             }
 290         }
 291     }
 292
 293     if (gip == NULL)
 294         return gip;
 295
 296     // If we have child nodes, check to see if any of them apply.
 297     pj_gi * child = find_grid(lam, phi, gip->children.begin(), gip->children.end());
 298     if (child != NULL)
 299         gip = child;
 300
 301     // if we get this far we have found a suitable grid
 302     return gip;
 303 }
 304
 305
 306 template <typename StreamPolicy>
 307 inline bool load_grid(StreamPolicy const& stream_policy, pj_gi_load & gi)
 308 {
 309     // load the grid shift info if we don't have it.
 310     if (gi.ct.cvs.empty())
 311     {
 312         typename StreamPolicy::stream_type is;
 313         stream_policy.open(is, gi.gridname);
 314
 315         if (! pj_gridinfo_load(is, gi))
 316         {
 317             //pj_ctx_set_errno( ctx, PJD_ERR_FAILED_TO_LOAD_GRID );
 318             return false;
 319         }
 320     }
 321
 322     return true;
 323 }
 324
 325
 326 /************************************************************************/
 327 /*                        pj_apply_gridshift_3()                        */
 328 /*                                                                      */
 329 /*      This is the real workhorse, given a gridlist.                   */
 330 /************************************************************************/
 331
 332 // Generic stream policy and standard grids
 333 template <bool Inverse, typename CalcT, typename StreamPolicy, typename Range, typename Grids>
 334 inline bool pj_apply_gridshift_3(StreamPolicy const& stream_policy,
 335                                  Range & range,
 336                                  Grids & grids,
 337                                  std::vector<std::size_t> const& gridindexes,
 338                                  grids_tag)
 339 {
 340     typedef typename boost::range_size<Range>::type size_type;
 341
 342     // If the grids are empty the indexes are as well
 343     if (gridindexes.empty())
 344     {
 345         //pj_ctx_set_errno(ctx, PJD_ERR_FAILED_TO_LOAD_GRID);
 346         //return PJD_ERR_FAILED_TO_LOAD_GRID;
 347         return false;
 348     }
 349
 350     size_type point_count = boost::size(range);
 351
 352     for (size_type i = 0 ; i < point_count ; ++i)
 353     {
 354         typename boost::range_reference<Range>::type
 355             point = range::at(range, i);
 356
 357         CalcT in_lon = geometry::get_as_radian<0>(point);
 358         CalcT in_lat = geometry::get_as_radian<1>(point);
 359
 360         pj_gi * gip = find_grid(in_lon, in_lat, grids.gridinfo, gridindexes);
 361
 362         if ( gip != NULL )
 363         {
 364             // load the grid shift info if we don't have it.
 365             if (! gip->ct.cvs.empty() || load_grid(stream_policy, *gip))
 366             {
 367                 // TODO: use set_invalid_point() or similar mechanism
 368                 CalcT out_lon = HUGE_VAL;
 369                 CalcT out_lat = HUGE_VAL;
 370
 371                 nad_cvt<Inverse>(in_lon, in_lat, out_lon, out_lat, *gip);
 372
 373                 // TODO: check differently
 374                 if ( out_lon != HUGE_VAL )
 375                 {
 376                     geometry::set_from_radian<0>(point, out_lon);
 377                     geometry::set_from_radian<1>(point, out_lat);
 378                 }
 379             }
 380         }
 381     }
 382
 383     return true;
 384 }
 385
 386 // Generic stream policy and shared grids
 387 template <bool Inverse, typename CalcT, typename StreamPolicy, typename Range, typename SharedGrids>
 388 inline bool pj_apply_gridshift_3(StreamPolicy const& stream_policy,
 389                                  Range & range,
 390                                  SharedGrids & grids,
 391                                  std::vector<std::size_t> const& gridindexes,
 392                                  shared_grids_tag)
 393 {
 394     typedef typename boost::range_size<Range>::type size_type;
 395
 396     // If the grids are empty the indexes are as well
 397     if (gridindexes.empty())
 398     {
 399         //pj_ctx_set_errno(ctx, PJD_ERR_FAILED_TO_LOAD_GRID);
 400         //return PJD_ERR_FAILED_TO_LOAD_GRID;
 401         return false;
 402     }
 403
 404     size_type point_count = boost::size(range);
 405
 406     // local storage
 407     pj_gi_load local_gi;
 408
 409     for (size_type i = 0 ; i < point_count ; )
 410     {
 411         bool load_needed = false;
 412
 413         CalcT in_lon = 0;
 414         CalcT in_lat = 0;
 415
 416         {
 417             typename SharedGrids::read_locked lck_grids(grids);
 418
 419             for ( ; i < point_count ; ++i )
 420             {
 421                 typename boost::range_reference<Range>::type
 422                     point = range::at(range, i);
 423
 424                 in_lon = geometry::get_as_radian<0>(point);
 425                 in_lat = geometry::get_as_radian<1>(point);
 426
 427                 pj_gi * gip = find_grid(in_lon, in_lat, lck_grids.gridinfo, gridindexes);
 428
 429                 if (gip == NULL)
 430                 {
 431                     // do nothing
 432                 }
 433                 else if (! gip->ct.cvs.empty())
 434                 {
 435                     // TODO: use set_invalid_point() or similar mechanism
 436                     CalcT out_lon = HUGE_VAL;
 437                     CalcT out_lat = HUGE_VAL;
 438
 439                     nad_cvt<Inverse>(in_lon, in_lat, out_lon, out_lat, *gip);
 440
 441                     // TODO: check differently
 442                     if (out_lon != HUGE_VAL)
 443                     {
 444                         geometry::set_from_radian<0>(point, out_lon);
 445                         geometry::set_from_radian<1>(point, out_lat);
 446                     }
 447                 }
 448                 else
 449                 {
 450                     // loading is needed
 451                     local_gi = *gip;
 452                     load_needed = true;
 453                     break;
 454                 }
 455             }
 456         }
 457
 458         if (load_needed)
 459         {
 460             if (load_grid(stream_policy, local_gi))
 461             {
 462                 typename SharedGrids::write_locked lck_grids(grids);
 463
 464                 // check again in case other thread already loaded the grid.
 465                 pj_gi * gip = find_grid(in_lon, in_lat, lck_grids.gridinfo, gridindexes);
 466
 467                 if (gip != NULL && gip->ct.cvs.empty())
 468                 {
 469                     // swap loaded local storage with empty grid
 470                     local_gi.swap(*gip);
 471                 }
 472             }
 473             else
 474             {
 475                 ++i;
 476             }
 477         }
 478     }
 479
 480     return true;
 481 }
 482
 483
 484 /************************************************************************/
 485 /*                        pj_apply_gridshift_2()                        */
 486 /*                                                                      */
 487 /*      This implementation uses the gridlist from a coordinate         */
 488 /*      system definition.  If the gridlist has not yet been            */
 489 /*      populated in the coordinate system definition we set it up      */
 490 /*      now.                                                            */
 491 /************************************************************************/
 492
 493 template <bool Inverse, typename Par, typename Range, typename ProjGrids>
 494 inline bool pj_apply_gridshift_2(Par const& /*defn*/, Range & range, ProjGrids const& proj_grids)
 495 {
 496     /*if( defn->catalog_name != NULL )
 497         return pj_gc_apply_gridshift( defn, inverse, point_count, point_offset,
 498                                       x, y, z );*/
 499
 500     /*std::vector<std::size_t> gridindexes;
 501     pj_gridlist_from_nadgrids(pj_get_param_s(defn.params, "nadgrids"),
 502                               grids.storage_ptr->stream_policy,
 503                               grids.storage_ptr->grids,
 504                               gridindexes);*/
 505
 506     // At this point the grids should be initialized
 507     if (proj_grids.hindexes.empty())
 508         return false;
 509
 510     return pj_apply_gridshift_3
 511             <
 512                 Inverse, typename Par::type
 513             >(proj_grids.grids_storage().stream_policy,
 514               range,
 515               proj_grids.grids_storage().hgrids,
 516               proj_grids.hindexes,
 517               typename ProjGrids::grids_storage_type::grids_type::tag());
 518 }
 519
 520 template <bool Inverse, typename Par, typename Range>
 521 inline bool pj_apply_gridshift_2(Par const& , Range & , srs::detail::empty_projection_grids const& )
 522 {
 523     return false;
 524 }
 525
 526
 527 } // namespace detail
 528
 529 }}} // namespace boost::geometry::projections
 530
 531 #endif // BOOST_GEOMETRY_SRS_PROJECTIONS_IMPL_PJ_APPLY_GRIDSHIFT_HPP