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_GRIDINFO_HPP
41 #define BOOST_GEOMETRY_SRS_PROJECTIONS_IMPL_PJ_GRIDINFO_HPP
42
43
44 #include <boost/algorithm/string.hpp>
45
46 #include <boost/geometry/core/assert.hpp>
47 #include <boost/geometry/util/math.hpp>
48
49 #include <boost/cstdint.hpp>
50
51 #include <algorithm>
52 #include <string>
53 #include <vector>
54
55
56 namespace boost { namespace geometry { namespace projections
57 {
58
59 namespace detail
60 {
61
62 /************************************************************************/
63 /* swap_words() */
64 /* */
65 /* Convert the byte order of the given word(s) in place. */
66 /************************************************************************/
67
is_lsb()68 inline bool is_lsb()
69 {
70 static const int byte_order_test = 1;
71 static bool result = (1 == ((const unsigned char *) (&byte_order_test))[0]);
72 return result;
73 }
74
swap_words(char * data,int word_size,int word_count)75 inline void swap_words( char *data, int word_size, int word_count )
76 {
77 for (int word = 0; word < word_count; word++)
78 {
79 for (int i = 0; i < word_size/2; i++)
80 {
81 std::swap(data[i], data[word_size-i-1]);
82 }
83
84 data += word_size;
85 }
86 }
87
cstr_equal(const char * s1,const char * s2,std::size_t n)88 inline bool cstr_equal(const char * s1, const char * s2, std::size_t n)
89 {
90 return std::equal(s1, s1 + n, s2);
91 }
92
93 struct is_trimmable_char
94 {
operator ()boost::geometry::projections::detail::is_trimmable_char95 inline bool operator()(char ch)
96 {
97 return ch == '\n' || ch == ' ';
98 }
99 };
100
101 // structs originally defined in projects.h
102
103 struct pj_ctable
104 {
105 struct lp_t { double lam, phi; };
106 struct flp_t { float lam, phi; };
107 struct ilp_t { boost::int32_t lam, phi; };
108
109 std::string id; // ascii info
110 lp_t ll; // lower left corner coordinates
111 lp_t del; // size of cells
112 ilp_t lim; // limits of conversion matrix
113 std::vector<flp_t> cvs; // conversion matrix
114
swapboost::geometry::projections::detail::pj_ctable115 inline void swap(pj_ctable & r)
116 {
117 id.swap(r.id);
118 std::swap(ll, r.ll);
119 std::swap(del, r.del);
120 std::swap(lim, r.lim);
121 cvs.swap(r.cvs);
122 }
123 };
124
125 struct pj_gi_load
126 {
127 enum format_t { missing = 0, ntv1, ntv2, gtx, ctable, ctable2 };
128 typedef boost::long_long_type offset_t;
129
pj_gi_loadboost::geometry::projections::detail::pj_gi_load130 explicit pj_gi_load(std::string const& gname = "",
131 format_t f = missing,
132 offset_t off = 0,
133 bool swap = false)
134 : gridname(gname)
135 , format(f)
136 , grid_offset(off)
137 , must_swap(swap)
138 {}
139
140 std::string gridname; // identifying name of grid, eg "conus" or ntv2_0.gsb
141
142 format_t format; // format of this grid, ie "ctable", "ntv1", "ntv2" or "missing".
143
144 offset_t grid_offset; // offset in file, for delayed loading
145 bool must_swap; // only for NTv2
146
147 pj_ctable ct;
148
swapboost::geometry::projections::detail::pj_gi_load149 inline void swap(pj_gi_load & r)
150 {
151 gridname.swap(r.gridname);
152 std::swap(format, r.format);
153 std::swap(grid_offset, r.grid_offset);
154 std::swap(must_swap, r.must_swap);
155 ct.swap(r.ct);
156 }
157
158 };
159
160 struct pj_gi
161 : pj_gi_load
162 {
pj_giboost::geometry::projections::detail::pj_gi163 explicit pj_gi(std::string const& gname = "",
164 pj_gi_load::format_t f = missing,
165 pj_gi_load::offset_t off = 0,
166 bool swap = false)
167 : pj_gi_load(gname, f, off, swap)
168 {}
169
170 std::vector<pj_gi> children;
171
swapboost::geometry::projections::detail::pj_gi172 inline void swap(pj_gi & r)
173 {
174 pj_gi_load::swap(r);
175 children.swap(r.children);
176 }
177 };
178
179 typedef std::vector<pj_gi> pj_gridinfo;
180
181
182 /************************************************************************/
183 /* pj_gridinfo_load_ctable() */
184 /* */
185 /* Load the data portion of a ctable formatted grid. */
186 /************************************************************************/
187
188 // Originally nad_ctable_load() defined in nad_init.c
189 template <typename IStream>
pj_gridinfo_load_ctable(IStream & is,pj_gi_load & gi)190 bool pj_gridinfo_load_ctable(IStream & is, pj_gi_load & gi)
191 {
192 pj_ctable & ct = gi.ct;
193
194 // Move the input stream by the size of the proj4 original CTABLE
195 std::size_t header_size = 80
196 + 2 * sizeof(pj_ctable::lp_t)
197 + sizeof(pj_ctable::ilp_t)
198 + sizeof(pj_ctable::flp_t*);
199 is.seekg(header_size);
200
201 // read all the actual shift values
202 std::size_t a_size = ct.lim.lam * ct.lim.phi;
203 ct.cvs.resize(a_size);
204
205 std::size_t ch_size = sizeof(pj_ctable::flp_t) * a_size;
206 is.read(reinterpret_cast<char*>(&ct.cvs[0]), ch_size);
207
208 if (is.fail() || std::size_t(is.gcount()) != ch_size)
209 {
210 ct.cvs.clear();
211 //ctable loading failed on fread() - binary incompatible?
212 return false;
213 }
214
215 return true;
216 }
217
218 /************************************************************************/
219 /* pj_gridinfo_load_ctable2() */
220 /* */
221 /* Load the data portion of a ctable2 formatted grid. */
222 /************************************************************************/
223
224 // Originally nad_ctable2_load() defined in nad_init.c
225 template <typename IStream>
pj_gridinfo_load_ctable2(IStream & is,pj_gi_load & gi)226 bool pj_gridinfo_load_ctable2(IStream & is, pj_gi_load & gi)
227 {
228 pj_ctable & ct = gi.ct;
229
230 is.seekg(160);
231
232 // read all the actual shift values
233 std::size_t a_size = ct.lim.lam * ct.lim.phi;
234 ct.cvs.resize(a_size);
235
236 std::size_t ch_size = sizeof(pj_ctable::flp_t) * a_size;
237 is.read(reinterpret_cast<char*>(&ct.cvs[0]), ch_size);
238
239 if (is.fail() || std::size_t(is.gcount()) != ch_size)
240 {
241 //ctable2 loading failed on fread() - binary incompatible?
242 ct.cvs.clear();
243 return false;
244 }
245
246 if (! is_lsb())
247 {
248 swap_words(reinterpret_cast<char*>(&ct.cvs[0]), 4, (int)a_size * 2);
249 }
250
251 return true;
252 }
253
254 /************************************************************************/
255 /* pj_gridinfo_load_ntv1() */
256 /* */
257 /* NTv1 format. */
258 /* We process one line at a time. Note that the array storage */
259 /* direction (e-w) is different in the NTv1 file and what */
260 /* the CTABLE is supposed to have. The phi/lam are also */
261 /* reversed, and we have to be aware of byte swapping. */
262 /************************************************************************/
263
264 // originally in pj_gridinfo_load() function
265 template <typename IStream>
pj_gridinfo_load_ntv1(IStream & is,pj_gi_load & gi)266 inline bool pj_gridinfo_load_ntv1(IStream & is, pj_gi_load & gi)
267 {
268 static const double s2r = math::d2r<double>() / 3600.0;
269
270 std::size_t const r_size = gi.ct.lim.lam * 2;
271 std::size_t const ch_size = sizeof(double) * r_size;
272
273 is.seekg(gi.grid_offset);
274
275 std::vector<double> row_buf(r_size);
276 gi.ct.cvs.resize(gi.ct.lim.lam * gi.ct.lim.phi);
277
278 for (boost::int32_t row = 0; row < gi.ct.lim.phi; row++ )
279 {
280 is.read(reinterpret_cast<char*>(&row_buf[0]), ch_size);
281
282 if (is.fail() || std::size_t(is.gcount()) != ch_size)
283 {
284 gi.ct.cvs.clear();
285 return false;
286 }
287
288 if (is_lsb())
289 swap_words(reinterpret_cast<char*>(&row_buf[0]), 8, (int)r_size);
290
291 // convert seconds to radians
292 for (boost::int32_t i = 0; i < gi.ct.lim.lam; i++ )
293 {
294 pj_ctable::flp_t & cvs = gi.ct.cvs[row * gi.ct.lim.lam + (gi.ct.lim.lam - i - 1)];
295
296 cvs.phi = (float) (row_buf[i*2] * s2r);
297 cvs.lam = (float) (row_buf[i*2+1] * s2r);
298 }
299 }
300
301 return true;
302 }
303
304 /* -------------------------------------------------------------------- */
305 /* pj_gridinfo_load_ntv2() */
306 /* */
307 /* NTv2 format. */
308 /* We process one line at a time. Note that the array storage */
309 /* direction (e-w) is different in the NTv2 file and what */
310 /* the CTABLE is supposed to have. The phi/lam are also */
311 /* reversed, and we have to be aware of byte swapping. */
312 /* -------------------------------------------------------------------- */
313
314 // originally in pj_gridinfo_load() function
315 template <typename IStream>
pj_gridinfo_load_ntv2(IStream & is,pj_gi_load & gi)316 inline bool pj_gridinfo_load_ntv2(IStream & is, pj_gi_load & gi)
317 {
318 static const double s2r = math::d2r<double>() / 3600.0;
319
320 std::size_t const r_size = gi.ct.lim.lam * 4;
321 std::size_t const ch_size = sizeof(float) * r_size;
322
323 is.seekg(gi.grid_offset);
324
325 std::vector<float> row_buf(r_size);
326 gi.ct.cvs.resize(gi.ct.lim.lam * gi.ct.lim.phi);
327
328 for (boost::int32_t row = 0; row < gi.ct.lim.phi; row++ )
329 {
330 is.read(reinterpret_cast<char*>(&row_buf[0]), ch_size);
331
332 if (is.fail() || std::size_t(is.gcount()) != ch_size)
333 {
334 gi.ct.cvs.clear();
335 return false;
336 }
337
338 if (gi.must_swap)
339 {
340 swap_words(reinterpret_cast<char*>(&row_buf[0]), 4, (int)r_size);
341 }
342
343 // convert seconds to radians
344 for (boost::int32_t i = 0; i < gi.ct.lim.lam; i++ )
345 {
346 pj_ctable::flp_t & cvs = gi.ct.cvs[row * gi.ct.lim.lam + (gi.ct.lim.lam - i - 1)];
347
348 // skip accuracy values
349 cvs.phi = (float) (row_buf[i*4] * s2r);
350 cvs.lam = (float) (row_buf[i*4+1] * s2r);
351 }
352 }
353
354 return true;
355 }
356
357 /************************************************************************/
358 /* pj_gridinfo_load_gtx() */
359 /* */
360 /* GTX format. */
361 /************************************************************************/
362
363 // originally in pj_gridinfo_load() function
364 template <typename IStream>
pj_gridinfo_load_gtx(IStream & is,pj_gi_load & gi)365 inline bool pj_gridinfo_load_gtx(IStream & is, pj_gi_load & gi)
366 {
367 boost::int32_t words = gi.ct.lim.lam * gi.ct.lim.phi;
368 std::size_t const ch_size = sizeof(float) * words;
369
370 is.seekg(gi.grid_offset);
371
372 // TODO: Consider changing this unintuitive code
373 // NOTE: Vertical shift data (one float per point) is stored in a container
374 // holding horizontal shift data (two floats per point).
375 gi.ct.cvs.resize((words + 1) / 2);
376
377 is.read(reinterpret_cast<char*>(&gi.ct.cvs[0]), ch_size);
378
379 if (is.fail() || std::size_t(is.gcount()) != ch_size)
380 {
381 gi.ct.cvs.clear();
382 return false;
383 }
384
385 if (is_lsb())
386 {
387 swap_words(reinterpret_cast<char*>(&gi.ct.cvs[0]), 4, words);
388 }
389
390 return true;
391 }
392
393 /************************************************************************/
394 /* pj_gridinfo_load() */
395 /* */
396 /* This function is intended to implement delayed loading of */
397 /* the data contents of a grid file. The header and related */
398 /* stuff are loaded by pj_gridinfo_init(). */
399 /************************************************************************/
400
401 template <typename IStream>
pj_gridinfo_load(IStream & is,pj_gi_load & gi)402 inline bool pj_gridinfo_load(IStream & is, pj_gi_load & gi)
403 {
404 if (! gi.ct.cvs.empty())
405 {
406 return true;
407 }
408
409 if (! is.is_open())
410 {
411 return false;
412 }
413
414 // Original platform specific CTable format.
415 if (gi.format == pj_gi::ctable)
416 {
417 return pj_gridinfo_load_ctable(is, gi);
418 }
419 // CTable2 format.
420 else if (gi.format == pj_gi::ctable2)
421 {
422 return pj_gridinfo_load_ctable2(is, gi);
423 }
424 // NTv1 format.
425 else if (gi.format == pj_gi::ntv1)
426 {
427 return pj_gridinfo_load_ntv1(is, gi);
428 }
429 // NTv2 format.
430 else if (gi.format == pj_gi::ntv2)
431 {
432 return pj_gridinfo_load_ntv2(is, gi);
433 }
434 // GTX format.
435 else if (gi.format == pj_gi::gtx)
436 {
437 return pj_gridinfo_load_gtx(is, gi);
438 }
439 else
440 {
441 return false;
442 }
443 }
444
445 /************************************************************************/
446 /* pj_gridinfo_parent() */
447 /* */
448 /* Seek a parent grid file by name from a grid list */
449 /************************************************************************/
450
451 template <typename It>
pj_gridinfo_parent(It first,It last,std::string const & name)452 inline It pj_gridinfo_parent(It first, It last, std::string const& name)
453 {
454 for ( ; first != last ; ++first)
455 {
456 if (first->ct.id == name)
457 return first;
458
459 It parent = pj_gridinfo_parent(first->children.begin(), first->children.end(), name);
460 if( parent != first->children.end() )
461 return parent;
462 }
463
464 return last;
465 }
466
467 /************************************************************************/
468 /* pj_gridinfo_init_ntv2() */
469 /* */
470 /* Load a ntv2 (.gsb) file. */
471 /************************************************************************/
472
473 template <typename IStream>
pj_gridinfo_init_ntv2(std::string const & gridname,IStream & is,pj_gridinfo & gridinfo)474 inline bool pj_gridinfo_init_ntv2(std::string const& gridname,
475 IStream & is,
476 pj_gridinfo & gridinfo)
477 {
478 BOOST_STATIC_ASSERT( sizeof(boost::int32_t) == 4 );
479 BOOST_STATIC_ASSERT( sizeof(double) == 8 );
480
481 static const double s2r = math::d2r<double>() / 3600.0;
482
483 std::size_t gridinfo_orig_size = gridinfo.size();
484
485 // Read the overview header.
486 char header[11*16];
487
488 is.read(header, sizeof(header));
489 if( is.fail() )
490 {
491 return false;
492 }
493
494 bool must_swap = (header[8] == 11)
495 ? !is_lsb()
496 : is_lsb();
497
498 // NOTE: This check is not implemented in proj4
499 if (! cstr_equal(header + 56, "SECONDS", 7))
500 {
501 return false;
502 }
503
504 // Byte swap interesting fields if needed.
505 if( must_swap )
506 {
507 swap_words( header+8, 4, 1 );
508 swap_words( header+8+16, 4, 1 );
509 swap_words( header+8+32, 4, 1 );
510 swap_words( header+8+7*16, 8, 1 );
511 swap_words( header+8+8*16, 8, 1 );
512 swap_words( header+8+9*16, 8, 1 );
513 swap_words( header+8+10*16, 8, 1 );
514 }
515
516 // Get the subfile count out ... all we really use for now.
517 boost::int32_t num_subfiles;
518 memcpy( &num_subfiles, header+8+32, 4 );
519
520 // Step through the subfiles, creating a PJ_GRIDINFO for each.
521 for( boost::int32_t subfile = 0; subfile < num_subfiles; subfile++ )
522 {
523 // Read header.
524 is.read(header, sizeof(header));
525 if( is.fail() )
526 {
527 return false;
528 }
529
530 if(! cstr_equal(header, "SUB_NAME", 8))
531 {
532 return false;
533 }
534
535 // Byte swap interesting fields if needed.
536 if( must_swap )
537 {
538 swap_words( header+8+16*4, 8, 1 );
539 swap_words( header+8+16*5, 8, 1 );
540 swap_words( header+8+16*6, 8, 1 );
541 swap_words( header+8+16*7, 8, 1 );
542 swap_words( header+8+16*8, 8, 1 );
543 swap_words( header+8+16*9, 8, 1 );
544 swap_words( header+8+16*10, 4, 1 );
545 }
546
547 // Initialize a corresponding "ct" structure.
548 pj_ctable ct;
549 pj_ctable::lp_t ur;
550
551 ct.id = std::string(header + 8, 8);
552
553 ct.ll.lam = - *((double *) (header+7*16+8)); /* W_LONG */
554 ct.ll.phi = *((double *) (header+4*16+8)); /* S_LAT */
555
556 ur.lam = - *((double *) (header+6*16+8)); /* E_LONG */
557 ur.phi = *((double *) (header+5*16+8)); /* N_LAT */
558
559 ct.del.lam = *((double *) (header+9*16+8));
560 ct.del.phi = *((double *) (header+8*16+8));
561
562 ct.lim.lam = (boost::int32_t) (fabs(ur.lam-ct.ll.lam)/ct.del.lam + 0.5) + 1;
563 ct.lim.phi = (boost::int32_t) (fabs(ur.phi-ct.ll.phi)/ct.del.phi + 0.5) + 1;
564
565 ct.ll.lam *= s2r;
566 ct.ll.phi *= s2r;
567 ct.del.lam *= s2r;
568 ct.del.phi *= s2r;
569
570 boost::int32_t gs_count;
571 memcpy( &gs_count, header + 8 + 16*10, 4 );
572 if( gs_count != ct.lim.lam * ct.lim.phi )
573 {
574 return false;
575 }
576
577 //ct.cvs.clear();
578
579 // Create a new gridinfo for this if we aren't processing the
580 // 1st subfile, and initialize our grid info.
581
582 // Attach to the correct list or sublist.
583
584 // TODO is offset needed?
585 pj_gi gi(gridname, pj_gi::ntv2, is.tellg(), must_swap);
586 gi.ct = ct;
587
588 if( subfile == 0 )
589 {
590 gridinfo.push_back(gi);
591 }
592 else if( cstr_equal(header+24, "NONE", 4) )
593 {
594 gridinfo.push_back(gi);
595 }
596 else
597 {
598 pj_gridinfo::iterator git = pj_gridinfo_parent(gridinfo.begin() + gridinfo_orig_size,
599 gridinfo.end(),
600 std::string((const char*)header+24, 8));
601
602 if( git == gridinfo.end() )
603 {
604 gridinfo.push_back(gi);
605 }
606 else
607 {
608 git->children.push_back(gi);
609 }
610 }
611
612 // Seek past the data.
613 is.seekg(gs_count * 16, std::ios::cur);
614 }
615
616 return true;
617 }
618
619 /************************************************************************/
620 /* pj_gridinfo_init_ntv1() */
621 /* */
622 /* Load an NTv1 style Canadian grid shift file. */
623 /************************************************************************/
624
625 template <typename IStream>
pj_gridinfo_init_ntv1(std::string const & gridname,IStream & is,pj_gridinfo & gridinfo)626 inline bool pj_gridinfo_init_ntv1(std::string const& gridname,
627 IStream & is,
628 pj_gridinfo & gridinfo)
629 {
630 BOOST_STATIC_ASSERT( sizeof(boost::int32_t) == 4 );
631 BOOST_STATIC_ASSERT( sizeof(double) == 8 );
632
633 static const double d2r = math::d2r<double>();
634
635 // Read the header.
636 char header[176];
637
638 is.read(header, sizeof(header));
639 if( is.fail() )
640 {
641 return false;
642 }
643
644 // Regularize fields of interest.
645 if( is_lsb() )
646 {
647 swap_words( header+8, 4, 1 );
648 swap_words( header+24, 8, 1 );
649 swap_words( header+40, 8, 1 );
650 swap_words( header+56, 8, 1 );
651 swap_words( header+72, 8, 1 );
652 swap_words( header+88, 8, 1 );
653 swap_words( header+104, 8, 1 );
654 }
655
656 // NTv1 grid shift file has wrong record count, corrupt?
657 if( *((boost::int32_t *) (header+8)) != 12 )
658 {
659 return false;
660 }
661
662 // NOTE: This check is not implemented in proj4
663 if (! cstr_equal(header + 120, "SECONDS", 7))
664 {
665 return false;
666 }
667
668 // Fill in CTABLE structure.
669 pj_ctable ct;
670 pj_ctable::lp_t ur;
671
672 ct.id = "NTv1 Grid Shift File";
673
674 ct.ll.lam = - *((double *) (header+72));
675 ct.ll.phi = *((double *) (header+24));
676 ur.lam = - *((double *) (header+56));
677 ur.phi = *((double *) (header+40));
678 ct.del.lam = *((double *) (header+104));
679 ct.del.phi = *((double *) (header+88));
680 ct.lim.lam = (boost::int32_t) (fabs(ur.lam-ct.ll.lam)/ct.del.lam + 0.5) + 1;
681 ct.lim.phi = (boost::int32_t) (fabs(ur.phi-ct.ll.phi)/ct.del.phi + 0.5) + 1;
682
683 ct.ll.lam *= d2r;
684 ct.ll.phi *= d2r;
685 ct.del.lam *= d2r;
686 ct.del.phi *= d2r;
687 //ct.cvs.clear();
688
689 // is offset needed?
690 gridinfo.push_back(pj_gi(gridname, pj_gi::ntv1, is.tellg()));
691 gridinfo.back().ct = ct;
692
693 return true;
694 }
695
696 /************************************************************************/
697 /* pj_gridinfo_init_gtx() */
698 /* */
699 /* Load a NOAA .gtx vertical datum shift file. */
700 /************************************************************************/
701
702 template <typename IStream>
pj_gridinfo_init_gtx(std::string const & gridname,IStream & is,pj_gridinfo & gridinfo)703 inline bool pj_gridinfo_init_gtx(std::string const& gridname,
704 IStream & is,
705 pj_gridinfo & gridinfo)
706 {
707 BOOST_STATIC_ASSERT( sizeof(boost::int32_t) == 4 );
708 BOOST_STATIC_ASSERT( sizeof(double) == 8 );
709
710 static const double d2r = math::d2r<double>();
711
712 // Read the header.
713 char header[40];
714
715 is.read(header, sizeof(header));
716 if( is.fail() )
717 {
718 return false;
719 }
720
721 // Regularize fields of interest and extract.
722 double xorigin, yorigin, xstep, ystep;
723 boost::int32_t rows, columns;
724
725 if( is_lsb() )
726 {
727 swap_words( header+0, 8, 4 );
728 swap_words( header+32, 4, 2 );
729 }
730
731 memcpy( &yorigin, header+0, 8 );
732 memcpy( &xorigin, header+8, 8 );
733 memcpy( &ystep, header+16, 8 );
734 memcpy( &xstep, header+24, 8 );
735
736 memcpy( &rows, header+32, 4 );
737 memcpy( &columns, header+36, 4 );
738
739 // gtx file header has invalid extents, corrupt?
740 if( xorigin < -360 || xorigin > 360
741 || yorigin < -90 || yorigin > 90 )
742 {
743 return false;
744 }
745
746 // Fill in CTABLE structure.
747 pj_ctable ct;
748
749 ct.id = "GTX Vertical Grid Shift File";
750
751 ct.ll.lam = xorigin;
752 ct.ll.phi = yorigin;
753 ct.del.lam = xstep;
754 ct.del.phi = ystep;
755 ct.lim.lam = columns;
756 ct.lim.phi = rows;
757
758 // some GTX files come in 0-360 and we shift them back into the
759 // expected -180 to 180 range if possible. This does not solve
760 // problems with grids spanning the dateline.
761 if( ct.ll.lam >= 180.0 )
762 ct.ll.lam -= 360.0;
763
764 if( ct.ll.lam >= 0.0 && ct.ll.lam + ct.del.lam * ct.lim.lam > 180.0 )
765 {
766 //"This GTX spans the dateline! This will cause problems." );
767 }
768
769 ct.ll.lam *= d2r;
770 ct.ll.phi *= d2r;
771 ct.del.lam *= d2r;
772 ct.del.phi *= d2r;
773 //ct.cvs.clear();
774
775 // is offset needed?
776 gridinfo.push_back(pj_gi(gridname, pj_gi::gtx, 40));
777 gridinfo.back().ct = ct;
778
779 return true;
780 }
781
782 /************************************************************************/
783 /* pj_gridinfo_init_ctable2() */
784 /* */
785 /* Read the header portion of a "ctable2" format grid. */
786 /************************************************************************/
787
788 // Originally nad_ctable2_init() defined in nad_init.c
789 template <typename IStream>
pj_gridinfo_init_ctable2(std::string const & gridname,IStream & is,pj_gridinfo & gridinfo)790 inline bool pj_gridinfo_init_ctable2(std::string const& gridname,
791 IStream & is,
792 pj_gridinfo & gridinfo)
793 {
794 BOOST_STATIC_ASSERT( sizeof(boost::int32_t) == 4 );
795 BOOST_STATIC_ASSERT( sizeof(double) == 8 );
796
797 char header[160];
798
799 is.read(header, sizeof(header));
800 if( is.fail() )
801 {
802 return false;
803 }
804
805 if( !is_lsb() )
806 {
807 swap_words( header + 96, 8, 4 );
808 swap_words( header + 128, 4, 2 );
809 }
810
811 // ctable2 - wrong header!
812 if (! cstr_equal(header, "CTABLE V2", 9))
813 {
814 return false;
815 }
816
817 // read the table header
818 pj_ctable ct;
819
820 ct.id = std::string(header + 16, std::find(header + 16, header + 16 + 80, '\0'));
821 //memcpy( &ct.ll.lam, header + 96, 8 );
822 //memcpy( &ct.ll.phi, header + 104, 8 );
823 //memcpy( &ct.del.lam, header + 112, 8 );
824 //memcpy( &ct.del.phi, header + 120, 8 );
825 //memcpy( &ct.lim.lam, header + 128, 4 );
826 //memcpy( &ct.lim.phi, header + 132, 4 );
827 memcpy( &ct.ll, header + 96, 40 );
828
829 // do some minimal validation to ensure the structure isn't corrupt
830 if ( (ct.lim.lam < 1) || (ct.lim.lam > 100000)
831 || (ct.lim.phi < 1) || (ct.lim.phi > 100000) )
832 {
833 return false;
834 }
835
836 // trim white space and newlines off id
837 boost::trim_right_if(ct.id, is_trimmable_char());
838
839 //ct.cvs.clear();
840
841 gridinfo.push_back(pj_gi(gridname, pj_gi::ctable2));
842 gridinfo.back().ct = ct;
843
844 return true;
845 }
846
847 /************************************************************************/
848 /* pj_gridinfo_init_ctable() */
849 /* */
850 /* Read the header portion of a "ctable" format grid. */
851 /************************************************************************/
852
853 // Originally nad_ctable_init() defined in nad_init.c
854 template <typename IStream>
pj_gridinfo_init_ctable(std::string const & gridname,IStream & is,pj_gridinfo & gridinfo)855 inline bool pj_gridinfo_init_ctable(std::string const& gridname,
856 IStream & is,
857 pj_gridinfo & gridinfo)
858 {
859 BOOST_STATIC_ASSERT( sizeof(boost::int32_t) == 4 );
860 BOOST_STATIC_ASSERT( sizeof(double) == 8 );
861
862 // 80 + 2*8 + 2*8 + 2*4
863 char header[120];
864
865 // NOTE: in proj4 data is loaded directly into CTABLE
866
867 is.read(header, sizeof(header));
868 if( is.fail() )
869 {
870 return false;
871 }
872
873 // NOTE: in proj4 LSB is not checked here
874
875 // read the table header
876 pj_ctable ct;
877
878 ct.id = std::string(header, std::find(header, header + 80, '\0'));
879 memcpy( &ct.ll, header + 80, 40 );
880
881 // do some minimal validation to ensure the structure isn't corrupt
882 if ( (ct.lim.lam < 1) || (ct.lim.lam > 100000)
883 || (ct.lim.phi < 1) || (ct.lim.phi > 100000) )
884 {
885 return false;
886 }
887
888 // trim white space and newlines off id
889 boost::trim_right_if(ct.id, is_trimmable_char());
890
891 //ct.cvs.clear();
892
893 gridinfo.push_back(pj_gi(gridname, pj_gi::ctable));
894 gridinfo.back().ct = ct;
895
896 return true;
897 }
898
899 /************************************************************************/
900 /* pj_gridinfo_init() */
901 /* */
902 /* Open and parse header details from a datum gridshift file */
903 /* returning a list of PJ_GRIDINFOs for the grids in that */
904 /* file. This superceeds use of nad_init() for modern */
905 /* applications. */
906 /************************************************************************/
907
908 template <typename IStream>
pj_gridinfo_init(std::string const & gridname,IStream & is,pj_gridinfo & gridinfo)909 inline bool pj_gridinfo_init(std::string const& gridname,
910 IStream & is,
911 pj_gridinfo & gridinfo)
912 {
913 char header[160];
914
915 // Check if the stream is opened.
916 if (! is.is_open()) {
917 return false;
918 }
919
920 // Load a header, to determine the file type.
921 is.read(header, sizeof(header));
922
923 if ( is.fail() ) {
924 return false;
925 }
926
927 is.seekg(0);
928
929 // Determine file type.
930 if ( cstr_equal(header + 0, "HEADER", 6)
931 && cstr_equal(header + 96, "W GRID", 6)
932 && cstr_equal(header + 144, "TO NAD83 ", 16) )
933 {
934 return pj_gridinfo_init_ntv1(gridname, is, gridinfo);
935 }
936 else if( cstr_equal(header + 0, "NUM_OREC", 8)
937 && cstr_equal(header + 48, "GS_TYPE", 7) )
938 {
939 return pj_gridinfo_init_ntv2(gridname, is, gridinfo);
940 }
941 else if( boost::algorithm::ends_with(gridname, "gtx")
942 || boost::algorithm::ends_with(gridname, "GTX") )
943 {
944 return pj_gridinfo_init_gtx(gridname, is, gridinfo);
945 }
946 else if( cstr_equal(header + 0, "CTABLE V2", 9) )
947 {
948 return pj_gridinfo_init_ctable2(gridname, is, gridinfo);
949 }
950 else
951 {
952 return pj_gridinfo_init_ctable(gridname, is, gridinfo);
953 }
954 }
955
956
957 } // namespace detail
958
959 }}} // namespace boost::geometry::projections
960
961 #endif // BOOST_GEOMETRY_SRS_PROJECTIONS_IMPL_PJ_GRIDINFO_HPP
962