pj_gridinfo.hpp 29 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711712713714715716717718719720721722723724725726727728729730731732733734735736737738739740741742743744745746747748749750751752753754755756757758759760761762763764765766767768769770771772773774775776777778779780781782783784785786787788789790791792793794795796797798799800801802803804805806807808809810811812813814815816817818819820821822823824825826827828829830831832833834835836837838839840841842843844845846847848849850851852853854855856857858859860861862863864865866867868869870871872873874875876877878879880881882883884885886887888889890891892893894895896897898899900901902903904905906907908909910911912913914915916917918919920921922923924925926927928929930931932933934935936937938939940941942943944945946947948949950951952953954955956957958959960961
  1. // Boost.Geometry
  2. // This file is manually converted from PROJ4
  3. // This file was modified by Oracle on 2018.
  4. // Modifications copyright (c) 2018, Oracle and/or its affiliates.
  5. // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
  6. // Use, modification and distribution is subject to the Boost Software License,
  7. // Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
  8. // http://www.boost.org/LICENSE_1_0.txt)
  9. // This file is converted from PROJ4, http://trac.osgeo.org/proj
  10. // PROJ4 is originally written by Gerald Evenden (then of the USGS)
  11. // PROJ4 is maintained by Frank Warmerdam
  12. // This file was converted to Geometry Library by Adam Wulkiewicz
  13. // Original copyright notice:
  14. // Author: Frank Warmerdam, warmerdam@pobox.com
  15. // Copyright (c) 2000, Frank Warmerdam
  16. // Permission is hereby granted, free of charge, to any person obtaining a
  17. // copy of this software and associated documentation files (the "Software"),
  18. // to deal in the Software without restriction, including without limitation
  19. // the rights to use, copy, modify, merge, publish, distribute, sublicense,
  20. // and/or sell copies of the Software, and to permit persons to whom the
  21. // Software is furnished to do so, subject to the following conditions:
  22. // The above copyright notice and this permission notice shall be included
  23. // in all copies or substantial portions of the Software.
  24. // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
  25. // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
  26. // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
  27. // THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
  28. // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
  29. // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
  30. // DEALINGS IN THE SOFTWARE.
  31. #ifndef BOOST_GEOMETRY_SRS_PROJECTIONS_IMPL_PJ_GRIDINFO_HPP
  32. #define BOOST_GEOMETRY_SRS_PROJECTIONS_IMPL_PJ_GRIDINFO_HPP
  33. #include <boost/algorithm/string.hpp>
  34. #include <boost/geometry/core/assert.hpp>
  35. #include <boost/cstdint.hpp>
  36. #include <algorithm>
  37. #include <string>
  38. #include <vector>
  39. namespace boost { namespace geometry { namespace projections
  40. {
  41. namespace detail
  42. {
  43. /************************************************************************/
  44. /* swap_words() */
  45. /* */
  46. /* Convert the byte order of the given word(s) in place. */
  47. /************************************************************************/
  48. inline bool is_lsb()
  49. {
  50. static const int byte_order_test = 1;
  51. static bool result = (1 == ((const unsigned char *) (&byte_order_test))[0]);
  52. return result;
  53. }
  54. inline void swap_words( char *data, int word_size, int word_count )
  55. {
  56. for (int word = 0; word < word_count; word++)
  57. {
  58. for (int i = 0; i < word_size/2; i++)
  59. {
  60. std::swap(data[i], data[word_size-i-1]);
  61. }
  62. data += word_size;
  63. }
  64. }
  65. inline bool cstr_equal(const char * s1, const char * s2, std::size_t n)
  66. {
  67. return std::equal(s1, s1 + n, s2);
  68. }
  69. struct is_trimmable_char
  70. {
  71. inline bool operator()(char ch)
  72. {
  73. return ch == '\n' || ch == ' ';
  74. }
  75. };
  76. // structs originally defined in projects.h
  77. struct pj_ctable
  78. {
  79. struct lp_t { double lam, phi; };
  80. struct flp_t { float lam, phi; };
  81. struct ilp_t { boost::int32_t lam, phi; };
  82. std::string id; // ascii info
  83. lp_t ll; // lower left corner coordinates
  84. lp_t del; // size of cells
  85. ilp_t lim; // limits of conversion matrix
  86. std::vector<flp_t> cvs; // conversion matrix
  87. inline void swap(pj_ctable & r)
  88. {
  89. id.swap(r.id);
  90. std::swap(ll, r.ll);
  91. std::swap(del, r.del);
  92. std::swap(lim, r.lim);
  93. cvs.swap(r.cvs);
  94. }
  95. };
  96. struct pj_gi_load
  97. {
  98. enum format_t { missing = 0, ntv1, ntv2, gtx, ctable, ctable2 };
  99. typedef boost::long_long_type offset_t;
  100. explicit pj_gi_load(std::string const& gname = "",
  101. format_t f = missing,
  102. offset_t off = 0,
  103. bool swap = false)
  104. : gridname(gname)
  105. , format(f)
  106. , grid_offset(off)
  107. , must_swap(swap)
  108. {}
  109. std::string gridname; // identifying name of grid, eg "conus" or ntv2_0.gsb
  110. format_t format; // format of this grid, ie "ctable", "ntv1", "ntv2" or "missing".
  111. offset_t grid_offset; // offset in file, for delayed loading
  112. bool must_swap; // only for NTv2
  113. pj_ctable ct;
  114. inline void swap(pj_gi_load & r)
  115. {
  116. gridname.swap(r.gridname);
  117. std::swap(format, r.format);
  118. std::swap(grid_offset, r.grid_offset);
  119. std::swap(must_swap, r.must_swap);
  120. ct.swap(r.ct);
  121. }
  122. };
  123. struct pj_gi
  124. : pj_gi_load
  125. {
  126. explicit pj_gi(std::string const& gname = "",
  127. pj_gi_load::format_t f = missing,
  128. pj_gi_load::offset_t off = 0,
  129. bool swap = false)
  130. : pj_gi_load(gname, f, off, swap)
  131. {}
  132. std::vector<pj_gi> children;
  133. inline void swap(pj_gi & r)
  134. {
  135. pj_gi_load::swap(r);
  136. children.swap(r.children);
  137. }
  138. };
  139. typedef std::vector<pj_gi> pj_gridinfo;
  140. /************************************************************************/
  141. /* pj_gridinfo_load_ctable() */
  142. /* */
  143. /* Load the data portion of a ctable formatted grid. */
  144. /************************************************************************/
  145. // Originally nad_ctable_load() defined in nad_init.c
  146. template <typename IStream>
  147. bool pj_gridinfo_load_ctable(IStream & is, pj_gi_load & gi)
  148. {
  149. pj_ctable & ct = gi.ct;
  150. // Move the input stream by the size of the proj4 original CTABLE
  151. std::size_t header_size = 80
  152. + 2 * sizeof(pj_ctable::lp_t)
  153. + sizeof(pj_ctable::ilp_t)
  154. + sizeof(pj_ctable::flp_t*);
  155. is.seekg(header_size);
  156. // read all the actual shift values
  157. std::size_t a_size = ct.lim.lam * ct.lim.phi;
  158. ct.cvs.resize(a_size);
  159. std::size_t ch_size = sizeof(pj_ctable::flp_t) * a_size;
  160. is.read(reinterpret_cast<char*>(&ct.cvs[0]), ch_size);
  161. if (is.fail() || is.gcount() != ch_size)
  162. {
  163. ct.cvs.clear();
  164. //ctable loading failed on fread() - binary incompatible?
  165. return false;
  166. }
  167. return true;
  168. }
  169. /************************************************************************/
  170. /* pj_gridinfo_load_ctable2() */
  171. /* */
  172. /* Load the data portion of a ctable2 formatted grid. */
  173. /************************************************************************/
  174. // Originally nad_ctable2_load() defined in nad_init.c
  175. template <typename IStream>
  176. bool pj_gridinfo_load_ctable2(IStream & is, pj_gi_load & gi)
  177. {
  178. pj_ctable & ct = gi.ct;
  179. is.seekg(160);
  180. // read all the actual shift values
  181. std::size_t a_size = ct.lim.lam * ct.lim.phi;
  182. ct.cvs.resize(a_size);
  183. std::size_t ch_size = sizeof(pj_ctable::flp_t) * a_size;
  184. is.read(reinterpret_cast<char*>(&ct.cvs[0]), ch_size);
  185. if (is.fail() || is.gcount() != ch_size)
  186. {
  187. //ctable2 loading failed on fread() - binary incompatible?
  188. ct.cvs.clear();
  189. return false;
  190. }
  191. if (! is_lsb())
  192. {
  193. swap_words(reinterpret_cast<char*>(&ct.cvs[0]), 4, (int)a_size * 2);
  194. }
  195. return true;
  196. }
  197. /************************************************************************/
  198. /* pj_gridinfo_load_ntv1() */
  199. /* */
  200. /* NTv1 format. */
  201. /* We process one line at a time. Note that the array storage */
  202. /* direction (e-w) is different in the NTv1 file and what */
  203. /* the CTABLE is supposed to have. The phi/lam are also */
  204. /* reversed, and we have to be aware of byte swapping. */
  205. /************************************************************************/
  206. // originally in pj_gridinfo_load() function
  207. template <typename IStream>
  208. inline bool pj_gridinfo_load_ntv1(IStream & is, pj_gi_load & gi)
  209. {
  210. static const double s2r = math::d2r<double>() / 3600.0;
  211. std::size_t const r_size = gi.ct.lim.lam * 2;
  212. std::size_t const ch_size = sizeof(double) * r_size;
  213. is.seekg(gi.grid_offset);
  214. std::vector<double> row_buf(r_size);
  215. gi.ct.cvs.resize(gi.ct.lim.lam * gi.ct.lim.phi);
  216. for (boost::int32_t row = 0; row < gi.ct.lim.phi; row++ )
  217. {
  218. is.read(reinterpret_cast<char*>(&row_buf[0]), ch_size);
  219. if (is.fail() || is.gcount() != ch_size)
  220. {
  221. gi.ct.cvs.clear();
  222. return false;
  223. }
  224. if (is_lsb())
  225. swap_words(reinterpret_cast<char*>(&row_buf[0]), 8, (int)r_size);
  226. // convert seconds to radians
  227. for (boost::int32_t i = 0; i < gi.ct.lim.lam; i++ )
  228. {
  229. pj_ctable::flp_t & cvs = gi.ct.cvs[row * gi.ct.lim.lam + (gi.ct.lim.lam - i - 1)];
  230. cvs.phi = (float) (row_buf[i*2] * s2r);
  231. cvs.lam = (float) (row_buf[i*2+1] * s2r);
  232. }
  233. }
  234. return true;
  235. }
  236. /* -------------------------------------------------------------------- */
  237. /* pj_gridinfo_load_ntv2() */
  238. /* */
  239. /* NTv2 format. */
  240. /* We process one line at a time. Note that the array storage */
  241. /* direction (e-w) is different in the NTv2 file and what */
  242. /* the CTABLE is supposed to have. The phi/lam are also */
  243. /* reversed, and we have to be aware of byte swapping. */
  244. /* -------------------------------------------------------------------- */
  245. // originally in pj_gridinfo_load() function
  246. template <typename IStream>
  247. inline bool pj_gridinfo_load_ntv2(IStream & is, pj_gi_load & gi)
  248. {
  249. static const double s2r = math::d2r<double>() / 3600.0;
  250. std::size_t const r_size = gi.ct.lim.lam * 4;
  251. std::size_t const ch_size = sizeof(float) * r_size;
  252. is.seekg(gi.grid_offset);
  253. std::vector<float> row_buf(r_size);
  254. gi.ct.cvs.resize(gi.ct.lim.lam * gi.ct.lim.phi);
  255. for (boost::int32_t row = 0; row < gi.ct.lim.phi; row++ )
  256. {
  257. is.read(reinterpret_cast<char*>(&row_buf[0]), ch_size);
  258. if (is.fail() || is.gcount() != ch_size)
  259. {
  260. gi.ct.cvs.clear();
  261. return false;
  262. }
  263. if (gi.must_swap)
  264. {
  265. swap_words(reinterpret_cast<char*>(&row_buf[0]), 4, (int)r_size);
  266. }
  267. // convert seconds to radians
  268. for (boost::int32_t i = 0; i < gi.ct.lim.lam; i++ )
  269. {
  270. pj_ctable::flp_t & cvs = gi.ct.cvs[row * gi.ct.lim.lam + (gi.ct.lim.lam - i - 1)];
  271. // skip accuracy values
  272. cvs.phi = (float) (row_buf[i*4] * s2r);
  273. cvs.lam = (float) (row_buf[i*4+1] * s2r);
  274. }
  275. }
  276. return true;
  277. }
  278. /************************************************************************/
  279. /* pj_gridinfo_load_gtx() */
  280. /* */
  281. /* GTX format. */
  282. /************************************************************************/
  283. // originally in pj_gridinfo_load() function
  284. template <typename IStream>
  285. inline bool pj_gridinfo_load_gtx(IStream & is, pj_gi_load & gi)
  286. {
  287. boost::int32_t words = gi.ct.lim.lam * gi.ct.lim.phi;
  288. std::size_t const ch_size = sizeof(float) * words;
  289. is.seekg(gi.grid_offset);
  290. // TODO: Consider changing this unintuitive code
  291. // NOTE: Vertical shift data (one float per point) is stored in a container
  292. // holding horizontal shift data (two floats per point).
  293. gi.ct.cvs.resize((words + 1) / 2);
  294. is.read(reinterpret_cast<char*>(&gi.ct.cvs[0]), ch_size);
  295. if (is.fail() || is.gcount() != ch_size)
  296. {
  297. gi.ct.cvs.clear();
  298. return false;
  299. }
  300. if (is_lsb())
  301. {
  302. swap_words(reinterpret_cast<char*>(&gi.ct.cvs[0]), 4, words);
  303. }
  304. return true;
  305. }
  306. /************************************************************************/
  307. /* pj_gridinfo_load() */
  308. /* */
  309. /* This function is intended to implement delayed loading of */
  310. /* the data contents of a grid file. The header and related */
  311. /* stuff are loaded by pj_gridinfo_init(). */
  312. /************************************************************************/
  313. template <typename IStream>
  314. inline bool pj_gridinfo_load(IStream & is, pj_gi_load & gi)
  315. {
  316. if (! gi.ct.cvs.empty())
  317. {
  318. return true;
  319. }
  320. if (! is.is_open())
  321. {
  322. return false;
  323. }
  324. // Original platform specific CTable format.
  325. if (gi.format == pj_gi::ctable)
  326. {
  327. return pj_gridinfo_load_ctable(is, gi);
  328. }
  329. // CTable2 format.
  330. else if (gi.format == pj_gi::ctable2)
  331. {
  332. return pj_gridinfo_load_ctable2(is, gi);
  333. }
  334. // NTv1 format.
  335. else if (gi.format == pj_gi::ntv1)
  336. {
  337. return pj_gridinfo_load_ntv1(is, gi);
  338. }
  339. // NTv2 format.
  340. else if (gi.format == pj_gi::ntv2)
  341. {
  342. return pj_gridinfo_load_ntv2(is, gi);
  343. }
  344. // GTX format.
  345. else if (gi.format == pj_gi::gtx)
  346. {
  347. return pj_gridinfo_load_gtx(is, gi);
  348. }
  349. else
  350. {
  351. return false;
  352. }
  353. }
  354. /************************************************************************/
  355. /* pj_gridinfo_parent() */
  356. /* */
  357. /* Seek a parent grid file by name from a grid list */
  358. /************************************************************************/
  359. template <typename It>
  360. inline It pj_gridinfo_parent(It first, It last, std::string const& name)
  361. {
  362. for ( ; first != last ; ++first)
  363. {
  364. if (first->ct.id == name)
  365. return first;
  366. It parent = pj_gridinfo_parent(first->children.begin(), first->children.end(), name);
  367. if( parent != first->children.end() )
  368. return parent;
  369. }
  370. return last;
  371. }
  372. /************************************************************************/
  373. /* pj_gridinfo_init_ntv2() */
  374. /* */
  375. /* Load a ntv2 (.gsb) file. */
  376. /************************************************************************/
  377. template <typename IStream>
  378. inline bool pj_gridinfo_init_ntv2(std::string const& gridname,
  379. IStream & is,
  380. pj_gridinfo & gridinfo)
  381. {
  382. BOOST_STATIC_ASSERT( sizeof(boost::int32_t) == 4 );
  383. BOOST_STATIC_ASSERT( sizeof(double) == 8 );
  384. static const double s2r = math::d2r<double>() / 3600.0;
  385. std::size_t gridinfo_orig_size = gridinfo.size();
  386. // Read the overview header.
  387. char header[11*16];
  388. is.read(header, sizeof(header));
  389. if( is.fail() )
  390. {
  391. return false;
  392. }
  393. bool must_swap = (header[8] == 11)
  394. ? !is_lsb()
  395. : is_lsb();
  396. // NOTE: This check is not implemented in proj4
  397. if (! cstr_equal(header + 56, "SECONDS", 7))
  398. {
  399. return false;
  400. }
  401. // Byte swap interesting fields if needed.
  402. if( must_swap )
  403. {
  404. swap_words( header+8, 4, 1 );
  405. swap_words( header+8+16, 4, 1 );
  406. swap_words( header+8+32, 4, 1 );
  407. swap_words( header+8+7*16, 8, 1 );
  408. swap_words( header+8+8*16, 8, 1 );
  409. swap_words( header+8+9*16, 8, 1 );
  410. swap_words( header+8+10*16, 8, 1 );
  411. }
  412. // Get the subfile count out ... all we really use for now.
  413. boost::int32_t num_subfiles;
  414. memcpy( &num_subfiles, header+8+32, 4 );
  415. // Step through the subfiles, creating a PJ_GRIDINFO for each.
  416. for( boost::int32_t subfile = 0; subfile < num_subfiles; subfile++ )
  417. {
  418. // Read header.
  419. is.read(header, sizeof(header));
  420. if( is.fail() )
  421. {
  422. return false;
  423. }
  424. if(! cstr_equal(header, "SUB_NAME", 8))
  425. {
  426. return false;
  427. }
  428. // Byte swap interesting fields if needed.
  429. if( must_swap )
  430. {
  431. swap_words( header+8+16*4, 8, 1 );
  432. swap_words( header+8+16*5, 8, 1 );
  433. swap_words( header+8+16*6, 8, 1 );
  434. swap_words( header+8+16*7, 8, 1 );
  435. swap_words( header+8+16*8, 8, 1 );
  436. swap_words( header+8+16*9, 8, 1 );
  437. swap_words( header+8+16*10, 4, 1 );
  438. }
  439. // Initialize a corresponding "ct" structure.
  440. pj_ctable ct;
  441. pj_ctable::lp_t ur;
  442. ct.id = std::string(header + 8, 8);
  443. ct.ll.lam = - *((double *) (header+7*16+8)); /* W_LONG */
  444. ct.ll.phi = *((double *) (header+4*16+8)); /* S_LAT */
  445. ur.lam = - *((double *) (header+6*16+8)); /* E_LONG */
  446. ur.phi = *((double *) (header+5*16+8)); /* N_LAT */
  447. ct.del.lam = *((double *) (header+9*16+8));
  448. ct.del.phi = *((double *) (header+8*16+8));
  449. ct.lim.lam = (boost::int32_t) (fabs(ur.lam-ct.ll.lam)/ct.del.lam + 0.5) + 1;
  450. ct.lim.phi = (boost::int32_t) (fabs(ur.phi-ct.ll.phi)/ct.del.phi + 0.5) + 1;
  451. ct.ll.lam *= s2r;
  452. ct.ll.phi *= s2r;
  453. ct.del.lam *= s2r;
  454. ct.del.phi *= s2r;
  455. boost::int32_t gs_count;
  456. memcpy( &gs_count, header + 8 + 16*10, 4 );
  457. if( gs_count != ct.lim.lam * ct.lim.phi )
  458. {
  459. return false;
  460. }
  461. //ct.cvs.clear();
  462. // Create a new gridinfo for this if we aren't processing the
  463. // 1st subfile, and initialize our grid info.
  464. // Attach to the correct list or sublist.
  465. // TODO is offset needed?
  466. pj_gi gi(gridname, pj_gi::ntv2, is.tellg(), must_swap);
  467. gi.ct = ct;
  468. if( subfile == 0 )
  469. {
  470. gridinfo.push_back(gi);
  471. }
  472. else if( cstr_equal(header+24, "NONE", 4) )
  473. {
  474. gridinfo.push_back(gi);
  475. }
  476. else
  477. {
  478. pj_gridinfo::iterator git = pj_gridinfo_parent(gridinfo.begin() + gridinfo_orig_size,
  479. gridinfo.end(),
  480. std::string((const char*)header+24, 8));
  481. if( git == gridinfo.end() )
  482. {
  483. gridinfo.push_back(gi);
  484. }
  485. else
  486. {
  487. git->children.push_back(gi);
  488. }
  489. }
  490. // Seek past the data.
  491. is.seekg(gs_count * 16, std::ios::cur);
  492. }
  493. return true;
  494. }
  495. /************************************************************************/
  496. /* pj_gridinfo_init_ntv1() */
  497. /* */
  498. /* Load an NTv1 style Canadian grid shift file. */
  499. /************************************************************************/
  500. template <typename IStream>
  501. inline bool pj_gridinfo_init_ntv1(std::string const& gridname,
  502. IStream & is,
  503. pj_gridinfo & gridinfo)
  504. {
  505. BOOST_STATIC_ASSERT( sizeof(boost::int32_t) == 4 );
  506. BOOST_STATIC_ASSERT( sizeof(double) == 8 );
  507. static const double d2r = math::d2r<double>();
  508. // Read the header.
  509. char header[176];
  510. is.read(header, sizeof(header));
  511. if( is.fail() )
  512. {
  513. return false;
  514. }
  515. // Regularize fields of interest.
  516. if( is_lsb() )
  517. {
  518. swap_words( header+8, 4, 1 );
  519. swap_words( header+24, 8, 1 );
  520. swap_words( header+40, 8, 1 );
  521. swap_words( header+56, 8, 1 );
  522. swap_words( header+72, 8, 1 );
  523. swap_words( header+88, 8, 1 );
  524. swap_words( header+104, 8, 1 );
  525. }
  526. // NTv1 grid shift file has wrong record count, corrupt?
  527. if( *((boost::int32_t *) (header+8)) != 12 )
  528. {
  529. return false;
  530. }
  531. // NOTE: This check is not implemented in proj4
  532. if (! cstr_equal(header + 120, "SECONDS", 7))
  533. {
  534. return false;
  535. }
  536. // Fill in CTABLE structure.
  537. pj_ctable ct;
  538. pj_ctable::lp_t ur;
  539. ct.id = "NTv1 Grid Shift File";
  540. ct.ll.lam = - *((double *) (header+72));
  541. ct.ll.phi = *((double *) (header+24));
  542. ur.lam = - *((double *) (header+56));
  543. ur.phi = *((double *) (header+40));
  544. ct.del.lam = *((double *) (header+104));
  545. ct.del.phi = *((double *) (header+88));
  546. ct.lim.lam = (boost::int32_t) (fabs(ur.lam-ct.ll.lam)/ct.del.lam + 0.5) + 1;
  547. ct.lim.phi = (boost::int32_t) (fabs(ur.phi-ct.ll.phi)/ct.del.phi + 0.5) + 1;
  548. ct.ll.lam *= d2r;
  549. ct.ll.phi *= d2r;
  550. ct.del.lam *= d2r;
  551. ct.del.phi *= d2r;
  552. //ct.cvs.clear();
  553. // is offset needed?
  554. gridinfo.push_back(pj_gi(gridname, pj_gi::ntv1, is.tellg()));
  555. gridinfo.back().ct = ct;
  556. return true;
  557. }
  558. /************************************************************************/
  559. /* pj_gridinfo_init_gtx() */
  560. /* */
  561. /* Load a NOAA .gtx vertical datum shift file. */
  562. /************************************************************************/
  563. template <typename IStream>
  564. inline bool pj_gridinfo_init_gtx(std::string const& gridname,
  565. IStream & is,
  566. pj_gridinfo & gridinfo)
  567. {
  568. BOOST_STATIC_ASSERT( sizeof(boost::int32_t) == 4 );
  569. BOOST_STATIC_ASSERT( sizeof(double) == 8 );
  570. static const double d2r = math::d2r<double>();
  571. // Read the header.
  572. char header[40];
  573. is.read(header, sizeof(header));
  574. if( is.fail() )
  575. {
  576. return false;
  577. }
  578. // Regularize fields of interest and extract.
  579. double xorigin, yorigin, xstep, ystep;
  580. boost::int32_t rows, columns;
  581. if( is_lsb() )
  582. {
  583. swap_words( header+0, 8, 4 );
  584. swap_words( header+32, 4, 2 );
  585. }
  586. memcpy( &yorigin, header+0, 8 );
  587. memcpy( &xorigin, header+8, 8 );
  588. memcpy( &ystep, header+16, 8 );
  589. memcpy( &xstep, header+24, 8 );
  590. memcpy( &rows, header+32, 4 );
  591. memcpy( &columns, header+36, 4 );
  592. // gtx file header has invalid extents, corrupt?
  593. if( xorigin < -360 || xorigin > 360
  594. || yorigin < -90 || yorigin > 90 )
  595. {
  596. return false;
  597. }
  598. // Fill in CTABLE structure.
  599. pj_ctable ct;
  600. ct.id = "GTX Vertical Grid Shift File";
  601. ct.ll.lam = xorigin;
  602. ct.ll.phi = yorigin;
  603. ct.del.lam = xstep;
  604. ct.del.phi = ystep;
  605. ct.lim.lam = columns;
  606. ct.lim.phi = rows;
  607. // some GTX files come in 0-360 and we shift them back into the
  608. // expected -180 to 180 range if possible. This does not solve
  609. // problems with grids spanning the dateline.
  610. if( ct.ll.lam >= 180.0 )
  611. ct.ll.lam -= 360.0;
  612. if( ct.ll.lam >= 0.0 && ct.ll.lam + ct.del.lam * ct.lim.lam > 180.0 )
  613. {
  614. //"This GTX spans the dateline! This will cause problems." );
  615. }
  616. ct.ll.lam *= d2r;
  617. ct.ll.phi *= d2r;
  618. ct.del.lam *= d2r;
  619. ct.del.phi *= d2r;
  620. //ct.cvs.clear();
  621. // is offset needed?
  622. gridinfo.push_back(pj_gi(gridname, pj_gi::gtx, 40));
  623. gridinfo.back().ct = ct;
  624. return true;
  625. }
  626. /************************************************************************/
  627. /* pj_gridinfo_init_ctable2() */
  628. /* */
  629. /* Read the header portion of a "ctable2" format grid. */
  630. /************************************************************************/
  631. // Originally nad_ctable2_init() defined in nad_init.c
  632. template <typename IStream>
  633. inline bool pj_gridinfo_init_ctable2(std::string const& gridname,
  634. IStream & is,
  635. pj_gridinfo & gridinfo)
  636. {
  637. BOOST_STATIC_ASSERT( sizeof(boost::int32_t) == 4 );
  638. BOOST_STATIC_ASSERT( sizeof(double) == 8 );
  639. char header[160];
  640. is.read(header, sizeof(header));
  641. if( is.fail() )
  642. {
  643. return false;
  644. }
  645. if( !is_lsb() )
  646. {
  647. swap_words( header + 96, 8, 4 );
  648. swap_words( header + 128, 4, 2 );
  649. }
  650. // ctable2 - wrong header!
  651. if (! cstr_equal(header, "CTABLE V2", 9))
  652. {
  653. return false;
  654. }
  655. // read the table header
  656. pj_ctable ct;
  657. ct.id = std::string(header + 16, std::find(header + 16, header + 16 + 80, '\0'));
  658. //memcpy( &ct.ll.lam, header + 96, 8 );
  659. //memcpy( &ct.ll.phi, header + 104, 8 );
  660. //memcpy( &ct.del.lam, header + 112, 8 );
  661. //memcpy( &ct.del.phi, header + 120, 8 );
  662. //memcpy( &ct.lim.lam, header + 128, 4 );
  663. //memcpy( &ct.lim.phi, header + 132, 4 );
  664. memcpy( &ct.ll, header + 96, 40 );
  665. // do some minimal validation to ensure the structure isn't corrupt
  666. if ( ct.lim.lam < 1 || ct.lim.lam > 100000
  667. || ct.lim.phi < 1 || ct.lim.phi > 100000 )
  668. {
  669. return false;
  670. }
  671. // trim white space and newlines off id
  672. boost::trim_right_if(ct.id, is_trimmable_char());
  673. //ct.cvs.clear();
  674. gridinfo.push_back(pj_gi(gridname, pj_gi::ctable2));
  675. gridinfo.back().ct = ct;
  676. return true;
  677. }
  678. /************************************************************************/
  679. /* pj_gridinfo_init_ctable() */
  680. /* */
  681. /* Read the header portion of a "ctable" format grid. */
  682. /************************************************************************/
  683. // Originally nad_ctable_init() defined in nad_init.c
  684. template <typename IStream>
  685. inline bool pj_gridinfo_init_ctable(std::string const& gridname,
  686. IStream & is,
  687. pj_gridinfo & gridinfo)
  688. {
  689. BOOST_STATIC_ASSERT( sizeof(boost::int32_t) == 4 );
  690. BOOST_STATIC_ASSERT( sizeof(double) == 8 );
  691. // 80 + 2*8 + 2*8 + 2*4
  692. char header[120];
  693. // NOTE: in proj4 data is loaded directly into CTABLE
  694. is.read(header, sizeof(header));
  695. if( is.fail() )
  696. {
  697. return false;
  698. }
  699. // NOTE: in proj4 LSB is not checked here
  700. // read the table header
  701. pj_ctable ct;
  702. ct.id = std::string(header, std::find(header, header + 80, '\0'));
  703. memcpy( &ct.ll, header + 80, 40 );
  704. // do some minimal validation to ensure the structure isn't corrupt
  705. if ( ct.lim.lam < 1 || ct.lim.lam > 100000
  706. || ct.lim.phi < 1 || ct.lim.phi > 100000 )
  707. {
  708. return false;
  709. }
  710. // trim white space and newlines off id
  711. boost::trim_right_if(ct.id, is_trimmable_char());
  712. //ct.cvs.clear();
  713. gridinfo.push_back(pj_gi(gridname, pj_gi::ctable));
  714. gridinfo.back().ct = ct;
  715. return true;
  716. }
  717. /************************************************************************/
  718. /* pj_gridinfo_init() */
  719. /* */
  720. /* Open and parse header details from a datum gridshift file */
  721. /* returning a list of PJ_GRIDINFOs for the grids in that */
  722. /* file. This superceeds use of nad_init() for modern */
  723. /* applications. */
  724. /************************************************************************/
  725. template <typename IStream>
  726. inline bool pj_gridinfo_init(std::string const& gridname,
  727. IStream & is,
  728. pj_gridinfo & gridinfo)
  729. {
  730. char header[160];
  731. // Check if the stream is opened.
  732. if (! is.is_open()) {
  733. return false;
  734. }
  735. // Load a header, to determine the file type.
  736. is.read(header, sizeof(header));
  737. if ( is.fail() ) {
  738. return false;
  739. }
  740. is.seekg(0);
  741. // Determine file type.
  742. if ( cstr_equal(header + 0, "HEADER", 6)
  743. && cstr_equal(header + 96, "W GRID", 6)
  744. && cstr_equal(header + 144, "TO NAD83 ", 16) )
  745. {
  746. return pj_gridinfo_init_ntv1(gridname, is, gridinfo);
  747. }
  748. else if( cstr_equal(header + 0, "NUM_OREC", 8)
  749. && cstr_equal(header + 48, "GS_TYPE", 7) )
  750. {
  751. return pj_gridinfo_init_ntv2(gridname, is, gridinfo);
  752. }
  753. else if( boost::algorithm::ends_with(gridname, "gtx")
  754. || boost::algorithm::ends_with(gridname, "GTX") )
  755. {
  756. return pj_gridinfo_init_gtx(gridname, is, gridinfo);
  757. }
  758. else if( cstr_equal(header + 0, "CTABLE V2", 9) )
  759. {
  760. return pj_gridinfo_init_ctable2(gridname, is, gridinfo);
  761. }
  762. else
  763. {
  764. return pj_gridinfo_init_ctable(gridname, is, gridinfo);
  765. }
  766. }
  767. } // namespace detail
  768. }}} // namespace boost::geometry::projections
  769. #endif // BOOST_GEOMETRY_SRS_PROJECTIONS_IMPL_PJ_GRIDINFO_HPP