projections/impl/pj_gridinfo.hpp

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 2018, 2019.
0007 // Modifications copyright (c) 2018-2019, 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 // Author:   Frank Warmerdam, warmerdam@pobox.com
0021
0022 // Copyright (c) 2000, Frank Warmerdam
0023
0024 // Permission is hereby granted, free of charge, to any person obtaining a
0025 // copy of this software and associated documentation files (the "Software"),
0026 // to deal in the Software without restriction, including without limitation
0027 // the rights to use, copy, modify, merge, publish, distribute, sublicense,
0028 // and/or sell copies of the Software, and to permit persons to whom the
0029 // Software is furnished to do so, subject to the following conditions:
0030
0031 // The above copyright notice and this permission notice shall be included
0032 // in all copies or substantial portions of the Software.
0033
0034 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
0035 // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
0036 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
0037 // THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
0038 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
0039 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
0040 // DEALINGS IN THE SOFTWARE.
0041
0042 #ifndef BOOST_GEOMETRY_SRS_PROJECTIONS_IMPL_PJ_GRIDINFO_HPP
0043 #define BOOST_GEOMETRY_SRS_PROJECTIONS_IMPL_PJ_GRIDINFO_HPP
0044
0045
0046 #include <boost/geometry/core/assert.hpp>
0047 #include <boost/geometry/util/math.hpp>
0048
0049 #include <boost/algorithm/string/predicate.hpp>
0050 #include <boost/algorithm/string/trim.hpp>
0051 #include <boost/cstdint.hpp>
0052
0053 #include <algorithm>
0054 #include <string>
0055 #include <vector>
0056
0057
0058 namespace boost { namespace geometry { namespace projections
0059 {
0060
0061 namespace detail
0062 {
0063
0064 /************************************************************************/
0065 /*                             swap_words()                             */
0066 /*                                                                      */
0067 /*      Convert the byte order of the given word(s) in place.           */
0068 /************************************************************************/
0069
0070 inline bool is_lsb()
0071 {
0072     static const int byte_order_test = 1;
0073     static bool result = (1 == ((const unsigned char *) (&byte_order_test))[0]);
0074     return result;
0075 }
0076
0077 inline void swap_words( char *data, int word_size, int word_count )
0078 {
0079     for (int word = 0; word < word_count; word++)
0080     {
0081         for (int i = 0; i < word_size/2; i++)
0082         {
0083             std::swap(data[i], data[word_size-i-1]);
0084         }
0085
0086         data += word_size;
0087     }
0088 }
0089
0090 inline bool cstr_equal(const char * s1, const char * s2, std::size_t n)
0091 {
0092     return std::equal(s1, s1 + n, s2);
0093 }
0094
0095 struct is_trimmable_char
0096 {
0097     inline bool operator()(char ch)
0098     {
0099         return ch == '\n' || ch == ' ';
0100     }
0101 };
0102
0103 // structs originally defined in projects.h
0104
0105 struct pj_ctable
0106 {
0107     struct lp_t  { double lam, phi; };
0108     struct flp_t { float lam, phi; };
0109     struct ilp_t { boost::int32_t lam, phi; };
0110
0111     std::string id;         // ascii info
0112     lp_t ll;                // lower left corner coordinates
0113     lp_t del;               // size of cells
0114     ilp_t lim;              // limits of conversion matrix
0115     std::vector<flp_t> cvs; // conversion matrix
0116
0117     inline void swap(pj_ctable & r)
0118     {
0119         id.swap(r.id);
0120         std::swap(ll, r.ll);
0121         std::swap(del, r.del);
0122         std::swap(lim, r.lim);
0123         cvs.swap(r.cvs);
0124     }
0125 };
0126
0127 struct pj_gi_load
0128 {
0129     enum format_t { missing = 0, ntv1, ntv2, gtx, ctable, ctable2 };
0130     typedef boost::long_long_type offset_t;
0131
0132     explicit pj_gi_load(std::string const& gname = "",
0133                         format_t f = missing,
0134                         offset_t off = 0,
0135                         bool swap = false)
0136         : gridname(gname)
0137         , format(f)
0138         , grid_offset(off)
0139         , must_swap(swap)
0140     {}
0141
0142     std::string gridname; // identifying name of grid, eg "conus" or ntv2_0.gsb
0143
0144     format_t format;      // format of this grid, ie "ctable", "ntv1", "ntv2" or "missing".
0145
0146     offset_t grid_offset; // offset in file, for delayed loading
0147     bool must_swap;       // only for NTv2
0148
0149     pj_ctable ct;
0150
0151     inline void swap(pj_gi_load & r)
0152     {
0153         gridname.swap(r.gridname);
0154         std::swap(format, r.format);
0155         std::swap(grid_offset, r.grid_offset);
0156         std::swap(must_swap, r.must_swap);
0157         ct.swap(r.ct);
0158     }
0159
0160 };
0161
0162 struct pj_gi
0163     : pj_gi_load
0164 {
0165     explicit pj_gi(std::string const& gname = "",
0166                    pj_gi_load::format_t f = missing,
0167                    pj_gi_load::offset_t off = 0,
0168                    bool swap = false)
0169         : pj_gi_load(gname, f, off, swap)
0170     {}
0171
0172     std::vector<pj_gi> children;
0173
0174     inline void swap(pj_gi & r)
0175     {
0176         pj_gi_load::swap(r);
0177         children.swap(r.children);
0178     }
0179 };
0180
0181 typedef std::vector<pj_gi> pj_gridinfo;
0182
0183
0184 /************************************************************************/
0185 /*                   pj_gridinfo_load_ctable()                          */
0186 /*                                                                      */
0187 /*      Load the data portion of a ctable formatted grid.               */
0188 /************************************************************************/
0189
0190 // Originally nad_ctable_load() defined in nad_init.c
0191 template <typename IStream>
0192 bool pj_gridinfo_load_ctable(IStream & is, pj_gi_load & gi)
0193 {
0194     pj_ctable & ct = gi.ct;
0195
0196     // Move the input stream by the size of the proj4 original CTABLE
0197     std::size_t header_size = 80
0198                             + 2 * sizeof(pj_ctable::lp_t)
0199                             + sizeof(pj_ctable::ilp_t)
0200                             + sizeof(pj_ctable::flp_t*);
0201     is.seekg(header_size);
0202
0203     // read all the actual shift values
0204     std::size_t a_size = ct.lim.lam * ct.lim.phi;
0205     ct.cvs.resize(a_size);
0206
0207     std::size_t ch_size = sizeof(pj_ctable::flp_t) * a_size;
0208     is.read(reinterpret_cast<char*>(&ct.cvs[0]), ch_size);
0209
0210     if (is.fail() || std::size_t(is.gcount()) != ch_size)
0211     {
0212         ct.cvs.clear();
0213         //ctable loading failed on fread() - binary incompatible?
0214         return false;
0215     }
0216
0217     return true;
0218 }
0219
0220 /************************************************************************/
0221 /*                  pj_gridinfo_load_ctable2()                          */
0222 /*                                                                      */
0223 /*      Load the data portion of a ctable2 formatted grid.              */
0224 /************************************************************************/
0225
0226 // Originally nad_ctable2_load() defined in nad_init.c
0227 template <typename IStream>
0228 bool pj_gridinfo_load_ctable2(IStream & is, pj_gi_load & gi)
0229 {
0230     pj_ctable & ct = gi.ct;
0231
0232     is.seekg(160);
0233
0234     // read all the actual shift values
0235     std::size_t a_size = ct.lim.lam * ct.lim.phi;
0236     ct.cvs.resize(a_size);
0237
0238     std::size_t ch_size = sizeof(pj_ctable::flp_t) * a_size;
0239     is.read(reinterpret_cast<char*>(&ct.cvs[0]), ch_size);
0240
0241     if (is.fail() || std::size_t(is.gcount()) != ch_size)
0242     {
0243         //ctable2 loading failed on fread() - binary incompatible?
0244         ct.cvs.clear();
0245         return false;
0246     }
0247
0248     if (! is_lsb())
0249     {
0250         swap_words(reinterpret_cast<char*>(&ct.cvs[0]), 4, (int)a_size * 2);
0251     }
0252
0253     return true;
0254 }
0255
0256 /************************************************************************/
0257 /*                      pj_gridinfo_load_ntv1()                         */
0258 /*                                                                      */
0259 /*      NTv1 format.                                                    */
0260 /*      We process one line at a time.  Note that the array storage     */
0261 /*      direction (e-w) is different in the NTv1 file and what          */
0262 /*      the CTABLE is supposed to have.  The phi/lam are also           */
0263 /*      reversed, and we have to be aware of byte swapping.             */
0264 /************************************************************************/
0265
0266 // originally in pj_gridinfo_load() function
0267 template <typename IStream>
0268 inline bool pj_gridinfo_load_ntv1(IStream & is, pj_gi_load & gi)
0269 {
0270     static const double s2r = math::d2r<double>() / 3600.0;
0271
0272     std::size_t const r_size = gi.ct.lim.lam * 2;
0273     std::size_t const ch_size = sizeof(double) * r_size;
0274
0275     is.seekg(gi.grid_offset);
0276
0277     std::vector<double> row_buf(r_size);
0278     gi.ct.cvs.resize(gi.ct.lim.lam * gi.ct.lim.phi);
0279
0280     for (boost::int32_t row = 0; row < gi.ct.lim.phi; row++ )
0281     {
0282         is.read(reinterpret_cast<char*>(&row_buf[0]), ch_size);
0283
0284         if (is.fail() || std::size_t(is.gcount()) != ch_size)
0285         {
0286             gi.ct.cvs.clear();
0287             return false;
0288         }
0289
0290         if (is_lsb())
0291             swap_words(reinterpret_cast<char*>(&row_buf[0]), 8, (int)r_size);
0292
0293         // convert seconds to radians
0294         for (boost::int32_t i = 0; i < gi.ct.lim.lam; i++ )
0295         {
0296             pj_ctable::flp_t & cvs = gi.ct.cvs[row * gi.ct.lim.lam + (gi.ct.lim.lam - i - 1)];
0297
0298             cvs.phi = (float) (row_buf[i*2] * s2r);
0299             cvs.lam = (float) (row_buf[i*2+1] * s2r);
0300         }
0301     }
0302
0303     return true;
0304 }
0305
0306 /* -------------------------------------------------------------------- */
0307 /*                         pj_gridinfo_load_ntv2()                      */
0308 /*                                                                      */
0309 /*      NTv2 format.                                                    */
0310 /*      We process one line at a time.  Note that the array storage     */
0311 /*      direction (e-w) is different in the NTv2 file and what          */
0312 /*      the CTABLE is supposed to have.  The phi/lam are also           */
0313 /*      reversed, and we have to be aware of byte swapping.             */
0314 /* -------------------------------------------------------------------- */
0315
0316 // originally in pj_gridinfo_load() function
0317 template <typename IStream>
0318 inline bool pj_gridinfo_load_ntv2(IStream & is, pj_gi_load & gi)
0319 {
0320     static const double s2r = math::d2r<double>() / 3600.0;
0321
0322     std::size_t const r_size = gi.ct.lim.lam * 4;
0323     std::size_t const ch_size = sizeof(float) * r_size;
0324
0325     is.seekg(gi.grid_offset);
0326
0327     std::vector<float> row_buf(r_size);
0328     gi.ct.cvs.resize(gi.ct.lim.lam * gi.ct.lim.phi);
0329
0330     for (boost::int32_t row = 0; row < gi.ct.lim.phi; row++ )
0331     {
0332         is.read(reinterpret_cast<char*>(&row_buf[0]), ch_size);
0333
0334         if (is.fail() || std::size_t(is.gcount()) != ch_size)
0335         {
0336             gi.ct.cvs.clear();
0337             return false;
0338         }
0339
0340         if (gi.must_swap)
0341         {
0342             swap_words(reinterpret_cast<char*>(&row_buf[0]), 4, (int)r_size);
0343         }
0344
0345         // convert seconds to radians
0346         for (boost::int32_t i = 0; i < gi.ct.lim.lam; i++ )
0347         {
0348             pj_ctable::flp_t & cvs = gi.ct.cvs[row * gi.ct.lim.lam + (gi.ct.lim.lam - i - 1)];
0349
0350             // skip accuracy values
0351             cvs.phi = (float) (row_buf[i*4] * s2r);
0352             cvs.lam = (float) (row_buf[i*4+1] * s2r);
0353         }
0354     }
0355
0356     return true;
0357 }
0358
0359 /************************************************************************/
0360 /*                   pj_gridinfo_load_gtx()                             */
0361 /*                                                                      */
0362 /*      GTX format.                                                     */
0363 /************************************************************************/
0364
0365 // originally in pj_gridinfo_load() function
0366 template <typename IStream>
0367 inline bool pj_gridinfo_load_gtx(IStream & is, pj_gi_load & gi)
0368 {
0369     boost::int32_t words = gi.ct.lim.lam * gi.ct.lim.phi;
0370     std::size_t const ch_size = sizeof(float) * words;
0371
0372     is.seekg(gi.grid_offset);
0373
0374     // TODO: Consider changing this unintuitive code
0375     // NOTE: Vertical shift data (one float per point) is stored in a container
0376     // holding horizontal shift data (two floats per point).
0377     gi.ct.cvs.resize((words + 1) / 2);
0378
0379     is.read(reinterpret_cast<char*>(&gi.ct.cvs[0]), ch_size);
0380
0381     if (is.fail() || std::size_t(is.gcount()) != ch_size)
0382     {
0383         gi.ct.cvs.clear();
0384         return false;
0385     }
0386
0387     if (is_lsb())
0388     {
0389         swap_words(reinterpret_cast<char*>(&gi.ct.cvs[0]), 4, words);
0390     }
0391
0392     return true;
0393 }
0394
0395 /************************************************************************/
0396 /*                          pj_gridinfo_load()                          */
0397 /*                                                                      */
0398 /*      This function is intended to implement delayed loading of       */
0399 /*      the data contents of a grid file.  The header and related       */
0400 /*      stuff are loaded by pj_gridinfo_init().                         */
0401 /************************************************************************/
0402
0403 template <typename IStream>
0404 inline bool pj_gridinfo_load(IStream & is, pj_gi_load & gi)
0405 {
0406     if (! gi.ct.cvs.empty())
0407     {
0408         return true;
0409     }
0410
0411     if (! is.is_open())
0412     {
0413         return false;
0414     }
0415
0416     // Original platform specific CTable format.
0417     if (gi.format == pj_gi::ctable)
0418     {
0419         return pj_gridinfo_load_ctable(is, gi);
0420     }
0421     // CTable2 format.
0422     else if (gi.format == pj_gi::ctable2)
0423     {
0424         return pj_gridinfo_load_ctable2(is, gi);
0425     }
0426     // NTv1 format.
0427     else if (gi.format == pj_gi::ntv1)
0428     {
0429         return pj_gridinfo_load_ntv1(is, gi);
0430     }
0431     // NTv2 format.
0432     else if (gi.format == pj_gi::ntv2)
0433     {
0434         return pj_gridinfo_load_ntv2(is, gi);
0435     }
0436     // GTX format.
0437     else if (gi.format == pj_gi::gtx)
0438     {
0439         return pj_gridinfo_load_gtx(is, gi);
0440     }
0441     else
0442     {
0443         return false;
0444     }
0445 }
0446
0447 /************************************************************************/
0448 /*                        pj_gridinfo_parent()                          */
0449 /*                                                                      */
0450 /*      Seek a parent grid file by name from a grid list                */
0451 /************************************************************************/
0452
0453 template <typename It>
0454 inline It pj_gridinfo_parent(It first, It last, std::string const& name)
0455 {
0456     for ( ; first != last ; ++first)
0457     {
0458         if (first->ct.id == name)
0459             return first;
0460
0461         It parent = pj_gridinfo_parent(first->children.begin(), first->children.end(), name);
0462         if( parent != first->children.end() )
0463             return parent;
0464     }
0465
0466     return last;
0467 }
0468
0469 /************************************************************************/
0470 /*                       pj_gridinfo_init_ntv2()                        */
0471 /*                                                                      */
0472 /*      Load a ntv2 (.gsb) file.                                        */
0473 /************************************************************************/
0474
0475 template <typename IStream>
0476 inline bool pj_gridinfo_init_ntv2(std::string const& gridname,
0477                                   IStream & is,
0478                                   pj_gridinfo & gridinfo)
0479 {
0480     BOOST_STATIC_ASSERT( sizeof(boost::int32_t) == 4 );
0481     BOOST_STATIC_ASSERT( sizeof(double) == 8 );
0482
0483     static const double s2r = math::d2r<double>() / 3600.0;
0484
0485     std::size_t gridinfo_orig_size = gridinfo.size();
0486
0487     // Read the overview header.
0488     char header[11*16];
0489
0490     is.read(header, sizeof(header));
0491     if( is.fail() )
0492     {
0493         return false;
0494     }
0495
0496     bool must_swap = (header[8] == 11)
0497                    ? !is_lsb()
0498                    : is_lsb();
0499
0500     // NOTE: This check is not implemented in proj4
0501     if (! cstr_equal(header + 56, "SECONDS", 7))
0502     {
0503         return false;
0504     }
0505
0506     // Byte swap interesting fields if needed.
0507     if( must_swap )
0508     {
0509         swap_words( header+8, 4, 1 );
0510         swap_words( header+8+16, 4, 1 );
0511         swap_words( header+8+32, 4, 1 );
0512         swap_words( header+8+7*16, 8, 1 );
0513         swap_words( header+8+8*16, 8, 1 );
0514         swap_words( header+8+9*16, 8, 1 );
0515         swap_words( header+8+10*16, 8, 1 );
0516     }
0517
0518     // Get the subfile count out ... all we really use for now.
0519     boost::int32_t num_subfiles;
0520     memcpy( &num_subfiles, header+8+32, 4 );
0521
0522     // Step through the subfiles, creating a PJ_GRIDINFO for each.
0523     for( boost::int32_t subfile = 0; subfile < num_subfiles; subfile++ )
0524     {
0525         // Read header.
0526         is.read(header, sizeof(header));
0527         if( is.fail() )
0528         {
0529             return false;
0530         }
0531
0532         if(! cstr_equal(header, "SUB_NAME", 8))
0533         {
0534             return false;
0535         }
0536
0537         // Byte swap interesting fields if needed.
0538         if( must_swap )
0539         {
0540             swap_words( header+8+16*4, 8, 1 );
0541             swap_words( header+8+16*5, 8, 1 );
0542             swap_words( header+8+16*6, 8, 1 );
0543             swap_words( header+8+16*7, 8, 1 );
0544             swap_words( header+8+16*8, 8, 1 );
0545             swap_words( header+8+16*9, 8, 1 );
0546             swap_words( header+8+16*10, 4, 1 );
0547         }
0548
0549         // Initialize a corresponding "ct" structure.
0550         pj_ctable ct;
0551         pj_ctable::lp_t ur;
0552
0553         ct.id = std::string(header + 8, 8);
0554
0555         ct.ll.lam = - *((double *) (header+7*16+8)); /* W_LONG */
0556         ct.ll.phi = *((double *) (header+4*16+8));   /* S_LAT */
0557
0558         ur.lam = - *((double *) (header+6*16+8));     /* E_LONG */
0559         ur.phi = *((double *) (header+5*16+8));       /* N_LAT */
0560
0561         ct.del.lam = *((double *) (header+9*16+8));
0562         ct.del.phi = *((double *) (header+8*16+8));
0563
0564         ct.lim.lam = (boost::int32_t) (fabs(ur.lam-ct.ll.lam)/ct.del.lam + 0.5) + 1;
0565         ct.lim.phi = (boost::int32_t) (fabs(ur.phi-ct.ll.phi)/ct.del.phi + 0.5) + 1;
0566
0567         ct.ll.lam *= s2r;
0568         ct.ll.phi *= s2r;
0569         ct.del.lam *= s2r;
0570         ct.del.phi *= s2r;
0571
0572         boost::int32_t gs_count;
0573         memcpy( &gs_count, header + 8 + 16*10, 4 );
0574         if( gs_count != ct.lim.lam * ct.lim.phi )
0575         {
0576             return false;
0577         }
0578
0579         //ct.cvs.clear();
0580
0581         // Create a new gridinfo for this if we aren't processing the
0582         // 1st subfile, and initialize our grid info.
0583
0584         // Attach to the correct list or sublist.
0585
0586         // TODO is offset needed?
0587         pj_gi gi(gridname, pj_gi::ntv2, is.tellg(), must_swap);
0588         gi.ct = ct;
0589
0590         if( subfile == 0 )
0591         {
0592             gridinfo.push_back(gi);
0593         }
0594         else if( cstr_equal(header+24, "NONE", 4) )
0595         {
0596             gridinfo.push_back(gi);
0597         }
0598         else
0599         {
0600             pj_gridinfo::iterator git = pj_gridinfo_parent(gridinfo.begin() + gridinfo_orig_size,
0601                                                            gridinfo.end(),
0602                                                            std::string((const char*)header+24, 8));
0603
0604             if( git == gridinfo.end() )
0605             {
0606                 gridinfo.push_back(gi);
0607             }
0608             else
0609             {
0610                 git->children.push_back(gi);
0611             }
0612         }
0613
0614         // Seek past the data.
0615         is.seekg(gs_count * 16, std::ios::cur);
0616     }
0617
0618     return true;
0619 }
0620
0621 /************************************************************************/
0622 /*                       pj_gridinfo_init_ntv1()                        */
0623 /*                                                                      */
0624 /*      Load an NTv1 style Canadian grid shift file.                    */
0625 /************************************************************************/
0626
0627 template <typename IStream>
0628 inline bool pj_gridinfo_init_ntv1(std::string const& gridname,
0629                                   IStream & is,
0630                                   pj_gridinfo & gridinfo)
0631 {
0632     BOOST_STATIC_ASSERT( sizeof(boost::int32_t) == 4 );
0633     BOOST_STATIC_ASSERT( sizeof(double) == 8 );
0634
0635     static const double d2r = math::d2r<double>();
0636
0637     // Read the header.
0638     char header[176];
0639
0640     is.read(header, sizeof(header));
0641     if( is.fail() )
0642     {
0643         return false;
0644     }
0645
0646     // Regularize fields of interest.
0647     if( is_lsb() )
0648     {
0649         swap_words( header+8, 4, 1 );
0650         swap_words( header+24, 8, 1 );
0651         swap_words( header+40, 8, 1 );
0652         swap_words( header+56, 8, 1 );
0653         swap_words( header+72, 8, 1 );
0654         swap_words( header+88, 8, 1 );
0655         swap_words( header+104, 8, 1 );
0656     }
0657
0658     // NTv1 grid shift file has wrong record count, corrupt?
0659     if( *((boost::int32_t *) (header+8)) != 12 )
0660     {
0661         return false;
0662     }
0663
0664     // NOTE: This check is not implemented in proj4
0665     if (! cstr_equal(header + 120, "SECONDS", 7))
0666     {
0667         return false;
0668     }
0669
0670     // Fill in CTABLE structure.
0671     pj_ctable ct;
0672     pj_ctable::lp_t ur;
0673
0674     ct.id = "NTv1 Grid Shift File";
0675
0676     ct.ll.lam = - *((double *) (header+72));
0677     ct.ll.phi = *((double *) (header+24));
0678     ur.lam = - *((double *) (header+56));
0679     ur.phi = *((double *) (header+40));
0680     ct.del.lam = *((double *) (header+104));
0681     ct.del.phi = *((double *) (header+88));
0682     ct.lim.lam = (boost::int32_t) (fabs(ur.lam-ct.ll.lam)/ct.del.lam + 0.5) + 1;
0683     ct.lim.phi = (boost::int32_t) (fabs(ur.phi-ct.ll.phi)/ct.del.phi + 0.5) + 1;
0684
0685     ct.ll.lam *= d2r;
0686     ct.ll.phi *= d2r;
0687     ct.del.lam *= d2r;
0688     ct.del.phi *= d2r;
0689     //ct.cvs.clear();
0690
0691     // is offset needed?
0692     gridinfo.push_back(pj_gi(gridname, pj_gi::ntv1, is.tellg()));
0693     gridinfo.back().ct = ct;
0694
0695     return true;
0696 }
0697
0698 /************************************************************************/
0699 /*                       pj_gridinfo_init_gtx()                         */
0700 /*                                                                      */
0701 /*      Load a NOAA .gtx vertical datum shift file.                     */
0702 /************************************************************************/
0703
0704 template <typename IStream>
0705 inline bool pj_gridinfo_init_gtx(std::string const& gridname,
0706                                  IStream & is,
0707                                  pj_gridinfo & gridinfo)
0708 {
0709     BOOST_STATIC_ASSERT( sizeof(boost::int32_t) == 4 );
0710     BOOST_STATIC_ASSERT( sizeof(double) == 8 );
0711
0712     static const double d2r = math::d2r<double>();
0713
0714     // Read the header.
0715     char header[40];
0716
0717     is.read(header, sizeof(header));
0718     if( is.fail() )
0719     {
0720         return false;
0721     }
0722
0723     // Regularize fields of interest and extract.
0724     double         xorigin, yorigin, xstep, ystep;
0725     boost::int32_t rows, columns;
0726
0727     if( is_lsb() )
0728     {
0729         swap_words( header+0, 8, 4 );
0730         swap_words( header+32, 4, 2 );
0731     }
0732
0733     memcpy( &yorigin, header+0, 8 );
0734     memcpy( &xorigin, header+8, 8 );
0735     memcpy( &ystep, header+16, 8 );
0736     memcpy( &xstep, header+24, 8 );
0737
0738     memcpy( &rows, header+32, 4 );
0739     memcpy( &columns, header+36, 4 );
0740
0741     // gtx file header has invalid extents, corrupt?
0742     if( xorigin < -360 || xorigin > 360
0743         || yorigin < -90 || yorigin > 90 )
0744     {
0745         return false;
0746     }
0747
0748     // Fill in CTABLE structure.
0749     pj_ctable ct;
0750
0751     ct.id = "GTX Vertical Grid Shift File";
0752
0753     ct.ll.lam = xorigin;
0754     ct.ll.phi = yorigin;
0755     ct.del.lam = xstep;
0756     ct.del.phi = ystep;
0757     ct.lim.lam = columns;
0758     ct.lim.phi = rows;
0759
0760     // some GTX files come in 0-360 and we shift them back into the
0761     // expected -180 to 180 range if possible. This does not solve
0762     // problems with grids spanning the dateline.
0763     if( ct.ll.lam >= 180.0 )
0764         ct.ll.lam -= 360.0;
0765
0766     if( ct.ll.lam >= 0.0 && ct.ll.lam + ct.del.lam * ct.lim.lam > 180.0 )
0767     {
0768         //"This GTX spans the dateline!  This will cause problems." );
0769     }
0770
0771     ct.ll.lam *= d2r;
0772     ct.ll.phi *= d2r;
0773     ct.del.lam *= d2r;
0774     ct.del.phi *= d2r;
0775     //ct.cvs.clear();
0776
0777     // is offset needed?
0778     gridinfo.push_back(pj_gi(gridname, pj_gi::gtx, 40));
0779     gridinfo.back().ct = ct;
0780
0781     return true;
0782 }
0783
0784 /************************************************************************/
0785 /*                   pj_gridinfo_init_ctable2()                         */
0786 /*                                                                      */
0787 /*      Read the header portion of a "ctable2" format grid.             */
0788 /************************************************************************/
0789
0790 // Originally nad_ctable2_init() defined in nad_init.c
0791 template <typename IStream>
0792 inline bool pj_gridinfo_init_ctable2(std::string const& gridname,
0793                                      IStream & is,
0794                                      pj_gridinfo & gridinfo)
0795 {
0796     BOOST_STATIC_ASSERT( sizeof(boost::int32_t) == 4 );
0797     BOOST_STATIC_ASSERT( sizeof(double) == 8 );
0798
0799     char header[160];
0800
0801     is.read(header, sizeof(header));
0802     if( is.fail() )
0803     {
0804         return false;
0805     }
0806
0807     if( !is_lsb() )
0808     {
0809         swap_words( header +  96, 8, 4 );
0810         swap_words( header + 128, 4, 2 );
0811     }
0812
0813     // ctable2 - wrong header!
0814     if (! cstr_equal(header, "CTABLE V2", 9))
0815     {
0816         return false;
0817     }
0818
0819     // read the table header
0820     pj_ctable ct;
0821
0822     ct.id = std::string(header + 16, std::find(header + 16, header + 16 + 80, '\0'));
0823     //memcpy( &ct.ll.lam,  header +  96, 8 );
0824     //memcpy( &ct.ll.phi,  header + 104, 8 );
0825     //memcpy( &ct.del.lam, header + 112, 8 );
0826     //memcpy( &ct.del.phi, header + 120, 8 );
0827     //memcpy( &ct.lim.lam, header + 128, 4 );
0828     //memcpy( &ct.lim.phi, header + 132, 4 );
0829     memcpy( &ct.ll,  header +  96, 40 );
0830
0831     // do some minimal validation to ensure the structure isn't corrupt
0832     if ( (ct.lim.lam < 1) || (ct.lim.lam > 100000)
0833       || (ct.lim.phi < 1) || (ct.lim.phi > 100000) )
0834     {
0835         return false;
0836     }
0837
0838     // trim white space and newlines off id
0839     boost::trim_right_if(ct.id, is_trimmable_char());
0840
0841     //ct.cvs.clear();
0842
0843     gridinfo.push_back(pj_gi(gridname, pj_gi::ctable2));
0844     gridinfo.back().ct = ct;
0845
0846     return true;
0847 }
0848
0849 /************************************************************************/
0850 /*                    pj_gridinfo_init_ctable()                         */
0851 /*                                                                      */
0852 /*      Read the header portion of a "ctable" format grid.              */
0853 /************************************************************************/
0854
0855 // Originally nad_ctable_init() defined in nad_init.c
0856 template <typename IStream>
0857 inline bool pj_gridinfo_init_ctable(std::string const& gridname,
0858                                     IStream & is,
0859                                     pj_gridinfo & gridinfo)
0860 {
0861     BOOST_STATIC_ASSERT( sizeof(boost::int32_t) == 4 );
0862     BOOST_STATIC_ASSERT( sizeof(double) == 8 );
0863
0864     // 80 + 2*8 + 2*8 + 2*4
0865     char header[120];
0866
0867     // NOTE: in proj4 data is loaded directly into CTABLE
0868
0869     is.read(header, sizeof(header));
0870     if( is.fail() )
0871     {
0872         return false;
0873     }
0874
0875     // NOTE: in proj4 LSB is not checked here
0876
0877     // read the table header
0878     pj_ctable ct;
0879
0880     ct.id = std::string(header, std::find(header, header + 80, '\0'));
0881     memcpy( &ct.ll, header + 80, 40 );
0882
0883     // do some minimal validation to ensure the structure isn't corrupt
0884     if ( (ct.lim.lam < 1) || (ct.lim.lam > 100000)
0885       || (ct.lim.phi < 1) || (ct.lim.phi > 100000) )
0886     {
0887         return false;
0888     }
0889
0890     // trim white space and newlines off id
0891     boost::trim_right_if(ct.id, is_trimmable_char());
0892
0893     //ct.cvs.clear();
0894
0895     gridinfo.push_back(pj_gi(gridname, pj_gi::ctable));
0896     gridinfo.back().ct = ct;
0897
0898     return true;
0899 }
0900
0901 /************************************************************************/
0902 /*                          pj_gridinfo_init()                          */
0903 /*                                                                      */
0904 /*      Open and parse header details from a datum gridshift file       */
0905 /*      returning a list of PJ_GRIDINFOs for the grids in that          */
0906 /*      file.  This superceeds use of nad_init() for modern             */
0907 /*      applications.                                                   */
0908 /************************************************************************/
0909
0910 template <typename IStream>
0911 inline bool pj_gridinfo_init(std::string const& gridname,
0912                              IStream & is,
0913                              pj_gridinfo & gridinfo)
0914 {
0915     char header[160];
0916
0917     // Check if the stream is opened.
0918     if (! is.is_open()) {
0919         return false;
0920     }
0921
0922     // Load a header, to determine the file type.
0923     is.read(header, sizeof(header));
0924
0925     if ( is.fail() ) {
0926         return false;
0927     }
0928
0929     is.seekg(0);
0930
0931     // Determine file type.
0932     if ( cstr_equal(header + 0, "HEADER", 6)
0933       && cstr_equal(header + 96, "W GRID", 6)
0934       && cstr_equal(header + 144, "TO      NAD83   ", 16) )
0935     {
0936         return pj_gridinfo_init_ntv1(gridname, is, gridinfo);
0937     }
0938     else if( cstr_equal(header + 0, "NUM_OREC", 8)
0939           && cstr_equal(header + 48, "GS_TYPE", 7) )
0940     {
0941         return pj_gridinfo_init_ntv2(gridname, is, gridinfo);
0942     }
0943     else if( boost::algorithm::ends_with(gridname, "gtx")
0944           || boost::algorithm::ends_with(gridname, "GTX") )
0945     {
0946         return pj_gridinfo_init_gtx(gridname, is, gridinfo);
0947     }
0948     else if( cstr_equal(header + 0, "CTABLE V2", 9) )
0949     {
0950         return pj_gridinfo_init_ctable2(gridname, is, gridinfo);
0951     }
0952     else
0953     {
0954         return pj_gridinfo_init_ctable(gridname, is, gridinfo);
0955     }
0956 }
0957
0958
0959 } // namespace detail
0960
0961 }}} // namespace boost::geometry::projections
0962
0963 #endif // BOOST_GEOMETRY_SRS_PROJECTIONS_IMPL_PJ_GRIDINFO_HPP