1 : /******************************************************************************
2 : *
3 : * Purpose: Implementation of the CPCIDSKGeoref class.
4 : *
5 : ******************************************************************************
6 : * Copyright (c) 2009
7 : * PCI Geomatics, 50 West Wilmot Street, Richmond Hill, Ont, Canada
8 : *
9 : * Permission is hereby granted, free of charge, to any person obtaining a
10 : * copy of this software and associated documentation files (the "Software"),
11 : * to deal in the Software without restriction, including without limitation
12 : * the rights to use, copy, modify, merge, publish, distribute, sublicense,
13 : * and/or sell copies of the Software, and to permit persons to whom the
14 : * Software is furnished to do so, subject to the following conditions:
15 : *
16 : * The above copyright notice and this permission notice shall be included
17 : * in all copies or substantial portions of the Software.
18 : *
19 : * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
20 : * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
21 : * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
22 : * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
23 : * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
24 : * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
25 : * DEALINGS IN THE SOFTWARE.
26 : ****************************************************************************/
27 :
28 : #include "pcidsk_exception.h"
29 : #include "segment/cpcidskgeoref.h"
30 : #include "core/pcidsk_utils.h"
31 : #include <cassert>
32 : #include <cstring>
33 : #include <cstdlib>
34 : #include <cstdio>
35 : #include <cctype>
36 :
37 : using namespace PCIDSK;
38 :
39 : static double PAK2PCI( double deg, int function );
40 :
41 : #ifndef ABS
42 : # define ABS(x) ((x<0) ? (-1*(x)) : x)
43 : #endif
44 :
45 : /************************************************************************/
46 : /* CPCIDSKGeoref() */
47 : /************************************************************************/
48 :
49 : CPCIDSKGeoref::CPCIDSKGeoref( PCIDSKFile *file, int segment,
50 88 : const char *segment_pointer )
51 88 : : CPCIDSKSegment( file, segment, segment_pointer )
52 :
53 : {
54 88 : loaded = false;
55 88 : }
56 :
57 : /************************************************************************/
58 : /* ~CPCIDSKGeoref() */
59 : /************************************************************************/
60 :
61 88 : CPCIDSKGeoref::~CPCIDSKGeoref()
62 :
63 : {
64 88 : }
65 :
66 : /************************************************************************/
67 : /* Initialize() */
68 : /************************************************************************/
69 :
70 44 : void CPCIDSKGeoref::Initialize()
71 :
72 : {
73 : // Note: we depend on Load() reacting gracefully to an uninitialized
74 : // georeferencing segment.
75 :
76 44 : WriteSimple( "PIXEL", 0.0, 1.0, 0.0, 0.0, 0.0, 1.0 );
77 44 : }
78 :
79 : /************************************************************************/
80 : /* Load() */
81 : /************************************************************************/
82 :
83 261 : void CPCIDSKGeoref::Load()
84 :
85 : {
86 261 : if( loaded )
87 0 : return;
88 :
89 : // TODO: this should likely be protected by a mutex.
90 :
91 : /* -------------------------------------------------------------------- */
92 : /* Load the segment contents into a buffer. */
93 : /* -------------------------------------------------------------------- */
94 261 : seg_data.SetSize( (int) (data_size - 1024) );
95 :
96 261 : ReadFromFile( seg_data.buffer, 0, data_size - 1024 );
97 :
98 : /* -------------------------------------------------------------------- */
99 : /* Handle simple case of a POLYNOMIAL. */
100 : /* -------------------------------------------------------------------- */
101 261 : if( strncmp(seg_data.buffer,"POLYNOMIAL",10) == 0 )
102 : {
103 4 : seg_data.Get(32,16,geosys);
104 :
105 4 : if( seg_data.GetInt(48,8) != 3 || seg_data.GetInt(56,8) != 3 )
106 0 : ThrowPCIDSKException( "Unexpected number of coefficients in POLYNOMIAL GEO segment." );
107 :
108 4 : a1 = seg_data.GetDouble(212+26*0,26);
109 4 : a2 = seg_data.GetDouble(212+26*1,26);
110 4 : xrot = seg_data.GetDouble(212+26*2,26);
111 :
112 4 : b1 = seg_data.GetDouble(1642+26*0,26);
113 4 : yrot = seg_data.GetDouble(1642+26*1,26);
114 4 : b3 = seg_data.GetDouble(1642+26*2,26);
115 : }
116 :
117 : /* -------------------------------------------------------------------- */
118 : /* Handle the case of a PROJECTION segment - for now we ignore */
119 : /* the actual projection parameters. */
120 : /* -------------------------------------------------------------------- */
121 257 : else if( strncmp(seg_data.buffer,"PROJECTION",10) == 0 )
122 : {
123 213 : seg_data.Get(32,16,geosys);
124 :
125 213 : if( seg_data.GetInt(48,8) != 3 || seg_data.GetInt(56,8) != 3 )
126 0 : ThrowPCIDSKException( "Unexpected number of coefficients in POLYNOMIAL GEO segment." );
127 :
128 213 : a1 = seg_data.GetDouble(1980+26*0,26);
129 213 : a2 = seg_data.GetDouble(1980+26*1,26);
130 213 : xrot = seg_data.GetDouble(1980+26*2,26);
131 :
132 213 : b1 = seg_data.GetDouble(2526+26*0,26);
133 213 : yrot = seg_data.GetDouble(2526+26*1,26);
134 213 : b3 = seg_data.GetDouble(2526+26*2,26);
135 : }
136 :
137 : /* -------------------------------------------------------------------- */
138 : /* Blank segment, just created and we just initialize things a bit.*/
139 : /* -------------------------------------------------------------------- */
140 44 : else if( memcmp(seg_data.buffer,
141 : "\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0\0",16) == 0 )
142 : {
143 44 : geosys = "";
144 :
145 44 : a1 = 0.0;
146 44 : a2 = 1.0;
147 44 : xrot = 0.0;
148 44 : b1 = 0.0;
149 44 : yrot = 0.0;
150 44 : b3 = 1.0;
151 : }
152 :
153 : else
154 : {
155 : ThrowPCIDSKException( "Unexpected GEO segment type: %s",
156 0 : seg_data.Get(0,16) );
157 : }
158 : }
159 :
160 : /************************************************************************/
161 : /* GetGeosys() */
162 : /************************************************************************/
163 :
164 44 : std::string CPCIDSKGeoref::GetGeosys()
165 :
166 : {
167 44 : Load();
168 44 : return geosys;
169 : }
170 :
171 : /************************************************************************/
172 : /* GetTransform() */
173 : /************************************************************************/
174 :
175 : void CPCIDSKGeoref::GetTransform( double &a1, double &a2, double &xrot,
176 45 : double &b1, double &yrot, double &b3 )
177 :
178 : {
179 45 : Load();
180 :
181 45 : a1 = this->a1;
182 45 : a2 = this->a2;
183 45 : xrot = this->xrot;
184 45 : b1 = this->b1;
185 45 : yrot = this->yrot;
186 45 : b3 = this->b3;
187 45 : }
188 :
189 : /************************************************************************/
190 : /* GetParameters() */
191 : /************************************************************************/
192 :
193 2 : std::vector<double> CPCIDSKGeoref::GetParameters()
194 :
195 : {
196 : unsigned int i;
197 2 : std::vector<double> parms;
198 :
199 2 : Load();
200 :
201 2 : parms.resize(18);
202 :
203 2 : if( strncmp(seg_data.buffer,"PROJECTION",10) != 0 )
204 : {
205 18 : for( i = 0; i < 17; i++ )
206 17 : parms[i] = 0.0;
207 1 : parms[17] = -1.0;
208 : }
209 : else
210 : {
211 18 : for( i = 0; i < 17; i++ )
212 17 : parms[i] = seg_data.GetDouble(80+26*i,26);
213 :
214 1 : std::string grid_units;
215 1 : seg_data.Get(64,16,grid_units);
216 :
217 1 : if( EQUALN(grid_units.c_str(),"DEGREE",3) )
218 0 : parms[17] = (double) (int) UNIT_DEGREE;
219 1 : else if( EQUALN(grid_units.c_str(),"MET",3) )
220 1 : parms[17] = (double) (int) UNIT_METER;
221 0 : else if( EQUALN(grid_units.c_str(),"FOOT",4) )
222 0 : parms[17] = (double) (int) UNIT_US_FOOT;
223 0 : else if( EQUALN(grid_units.c_str(),"FEET",4) )
224 0 : parms[17] = (double) (int) UNIT_US_FOOT;
225 0 : else if( EQUALN(grid_units.c_str(),"INTL FOOT",5) )
226 0 : parms[17] = (double) (int) UNIT_INTL_FOOT;
227 : else
228 0 : parms[17] = -1.0; /* unknown */
229 : }
230 :
231 0 : return parms;
232 : }
233 :
234 : /************************************************************************/
235 : /* WriteSimple() */
236 : /************************************************************************/
237 :
238 : void CPCIDSKGeoref::WriteSimple( std::string const& geosys,
239 : double a1, double a2, double xrot,
240 128 : double b1, double yrot, double b3 )
241 :
242 : {
243 128 : Load();
244 :
245 128 : std::string geosys_clean(ReformatGeosys( geosys ));
246 :
247 : /* -------------------------------------------------------------------- */
248 : /* Establish the appropriate units code when possible. */
249 : /* -------------------------------------------------------------------- */
250 128 : std::string units_code = "METER";
251 :
252 256 : if( EQUALN(geosys_clean.c_str(),"FOOT",4) )
253 0 : units_code = "FOOT";
254 128 : else if( EQUALN(geosys_clean.c_str(),"SPAF",4) )
255 0 : units_code = "FOOT";
256 128 : else if( EQUALN(geosys_clean.c_str(),"SPIF",4) )
257 0 : units_code = "INTL FOOT";
258 128 : else if( EQUALN(geosys_clean.c_str(),"LONG",4) )
259 41 : units_code = "DEEGREE";
260 :
261 : /* -------------------------------------------------------------------- */
262 : /* Write a fairly simple PROJECTION segment. */
263 : /* -------------------------------------------------------------------- */
264 128 : seg_data.SetSize( 6 * 512 );
265 :
266 128 : seg_data.Put( " ", 0, seg_data.buffer_size );
267 :
268 : // SD.PRO.P1
269 128 : seg_data.Put( "PROJECTION", 0, 16 );
270 :
271 : // SD.PRO.P2
272 128 : seg_data.Put( "PIXEL", 16, 16 );
273 :
274 : // SD.PRO.P3
275 128 : seg_data.Put( geosys_clean.c_str(), 32, 16 );
276 :
277 : // SD.PRO.P4
278 128 : seg_data.Put( 3, 48, 8 );
279 :
280 : // SD.PRO.P5
281 128 : seg_data.Put( 3, 56, 8 );
282 :
283 : // SD.PRO.P6
284 128 : seg_data.Put( units_code.c_str(), 64, 16 );
285 :
286 : // SD.PRO.P7 - P22
287 2304 : for( int i = 0; i < 17; i++ )
288 2176 : seg_data.Put( 0.0, 80 + i*26, 26, "%26.18E" );
289 :
290 : // SD.PRO.P24
291 128 : PrepareGCTPFields();
292 :
293 : // SD.PRO.P26
294 128 : seg_data.Put( a1, 1980 + 0*26, 26, "%26.18E" );
295 128 : seg_data.Put( a2, 1980 + 1*26, 26, "%26.18E" );
296 128 : seg_data.Put( xrot,1980 + 2*26, 26, "%26.18E" );
297 :
298 : // SD.PRO.P27
299 128 : seg_data.Put( b1, 2526 + 0*26, 26, "%26.18E" );
300 128 : seg_data.Put( yrot, 2526 + 1*26, 26, "%26.18E" );
301 128 : seg_data.Put( b3, 2526 + 2*26, 26, "%26.18E" );
302 :
303 128 : WriteToFile( seg_data.buffer, 0, seg_data.buffer_size );
304 :
305 128 : loaded = false;
306 128 : }
307 :
308 : /************************************************************************/
309 : /* WriteParameters() */
310 : /************************************************************************/
311 :
312 42 : void CPCIDSKGeoref::WriteParameters( std::vector<double> const& parms )
313 :
314 : {
315 42 : Load();
316 :
317 42 : if( parms.size() < 17 )
318 0 : ThrowPCIDSKException( "Did not get expected number of paramters in WriteParameters()" );
319 :
320 : unsigned int i;
321 :
322 756 : for( i = 0; i < 17; i++ )
323 714 : seg_data.Put(parms[i],80+26*i,26,"%26.16f");
324 :
325 42 : if( parms.size() >= 18 )
326 : {
327 42 : switch( (UnitCode) (int) parms[17] )
328 : {
329 : case UNIT_DEGREE:
330 41 : seg_data.Put( "DEGREE", 64, 16 );
331 41 : break;
332 :
333 : case UNIT_METER:
334 1 : seg_data.Put( "METER", 64, 16 );
335 1 : break;
336 :
337 : case UNIT_US_FOOT:
338 0 : seg_data.Put( "FOOT", 64, 16 );
339 0 : break;
340 :
341 : case UNIT_INTL_FOOT:
342 0 : seg_data.Put( "INTL FOOT", 64, 16 );
343 : break;
344 : }
345 : }
346 :
347 42 : PrepareGCTPFields();
348 :
349 42 : WriteToFile( seg_data.buffer, 0, seg_data.buffer_size );
350 :
351 : // no need to mark loaded false, since we don't cache these parameters.
352 42 : }
353 :
354 : /************************************************************************/
355 : /* GetUSGSParameters() */
356 : /************************************************************************/
357 :
358 0 : std::vector<double> CPCIDSKGeoref::GetUSGSParameters()
359 :
360 : {
361 : unsigned int i;
362 0 : std::vector<double> parms;
363 :
364 0 : Load();
365 :
366 0 : parms.resize(19);
367 0 : if( strncmp(seg_data.buffer,"PROJECTION",10) != 0 )
368 : {
369 0 : for( i = 0; i < 19; i++ )
370 0 : parms[i] = 0.0;
371 : }
372 : else
373 : {
374 0 : for( i = 0; i < 19; i++ )
375 0 : parms[i] = seg_data.GetDouble(1458+26*i,26);
376 : }
377 :
378 0 : return parms;
379 : }
380 :
381 : /************************************************************************/
382 : /* ReformatGeosys() */
383 : /* */
384 : /* Put a geosys string into standard form. Similar to what the */
385 : /* DecodeGeosys() function in the PCI SDK does. */
386 : /************************************************************************/
387 :
388 298 : std::string CPCIDSKGeoref::ReformatGeosys( std::string const& geosys )
389 :
390 : {
391 : /* -------------------------------------------------------------------- */
392 : /* Put into a local buffer and pad out to 16 characters with */
393 : /* spaces. */
394 : /* -------------------------------------------------------------------- */
395 : char local_buf[33];
396 :
397 298 : strncpy( local_buf, geosys.c_str(), 16 );
398 298 : local_buf[16] = '\0';
399 298 : strcat( local_buf, " " );
400 298 : local_buf[16] = '\0';
401 :
402 : /* -------------------------------------------------------------------- */
403 : /* Extract the earth model from the geosys string. */
404 : /* -------------------------------------------------------------------- */
405 : char earthmodel[5];
406 : const char *cp;
407 : int i;
408 : char last;
409 :
410 298 : cp = local_buf;
411 5066 : while( cp < local_buf + 16 && cp[1] != '\0' )
412 4470 : cp++;
413 :
414 2488 : while( cp > local_buf && isspace(*cp) )
415 1892 : cp--;
416 :
417 298 : last = '\0';
418 974 : while( cp > local_buf
419 : && (isdigit((unsigned char)*cp)
420 : || *cp == '-' || *cp == '+' ) )
421 : {
422 378 : if( last == '\0' ) last = *cp;
423 378 : cp--;
424 : }
425 :
426 424 : if( isdigit( (unsigned char)last ) &&
427 : ( *cp == 'D' || *cp == 'd' ||
428 : *cp == 'E' || *cp == 'e' ) )
429 : {
430 126 : i = atoi( cp+1 );
431 252 : if( i > -100 && i < 1000
432 : && (cp == local_buf
433 : || ( cp > local_buf && isspace( *(cp-1) ) )
434 : )
435 : )
436 : {
437 249 : if( *cp == 'D' || *cp == 'd' )
438 123 : sprintf( earthmodel, "D%03d", i );
439 : else
440 3 : sprintf( earthmodel, "E%03d", i );
441 : }
442 : else
443 : {
444 0 : sprintf( earthmodel, " " );
445 : }
446 : }
447 : else
448 : {
449 172 : sprintf( earthmodel, " " );
450 : }
451 :
452 : /* -------------------------------------------------------------------- */
453 : /* Identify by geosys string. */
454 : /* -------------------------------------------------------------------- */
455 : const char *ptr;
456 : int zone, ups_zone;
457 298 : char zone_code = ' ';
458 :
459 298 : if( EQUALN(local_buf,"PIX",3) )
460 : {
461 172 : strcpy( local_buf, "PIXEL " );
462 : }
463 126 : else if( EQUALN(local_buf,"UTM",3) )
464 : {
465 : /* Attempt to find a zone and ellipsoid */
466 3 : for( ptr=local_buf+3; isspace(*ptr); ptr++ ) {}
467 6 : if( isdigit( (unsigned char)*ptr ) || *ptr == '-' )
468 : {
469 3 : zone = atoi(ptr);
470 3 : for( ; isdigit((unsigned char)*ptr) || *ptr == '-'; ptr++ ) {}
471 3 : for( ; isspace(*ptr); ptr++ ) {}
472 3 : if( isalpha(*ptr) && !isdigit((unsigned char)*(ptr+1)) )
473 0 : zone_code = *(ptr++);
474 : }
475 : else
476 0 : zone = -100;
477 :
478 6 : if( zone >= -60 && zone <= 60 && zone != 0 )
479 : {
480 3 : if( zone_code >= 'a' && zone_code <= 'z' )
481 0 : zone_code = zone_code - 'a' + 'A';
482 :
483 3 : if( zone_code == ' ' && zone < 0 )
484 0 : zone_code = 'C';
485 :
486 3 : zone = ABS(zone);
487 :
488 : sprintf( local_buf,
489 : "UTM %3d %c %4s",
490 3 : zone, zone_code, earthmodel );
491 : }
492 : else
493 : {
494 : sprintf( local_buf,
495 : "UTM %4s",
496 0 : earthmodel );
497 : }
498 3 : if( local_buf[14] == ' ' )
499 0 : local_buf[14] = '0';
500 3 : if( local_buf[13] == ' ' )
501 0 : local_buf[13] = '0';
502 : }
503 123 : else if( EQUALN(local_buf,"MET",3) )
504 : {
505 0 : sprintf( local_buf, "METRE %4s", earthmodel );
506 : }
507 123 : else if( EQUALN(local_buf,"FEET",4) || EQUALN(local_buf,"FOOT",4))
508 : {
509 0 : sprintf( local_buf, "FOOT %4s", earthmodel );
510 : }
511 123 : else if( EQUALN(local_buf,"LAT",3) ||
512 : EQUALN(local_buf,"LON",3) )
513 : {
514 : sprintf( local_buf,
515 : "LONG/LAT %4s",
516 123 : earthmodel );
517 : }
518 0 : else if( EQUALN(local_buf,"SPCS ",5) ||
519 : EQUALN(local_buf,"SPAF ",5) ||
520 : EQUALN(local_buf,"SPIF ",5) )
521 : {
522 0 : int nSPZone = 0;
523 :
524 0 : for( ptr=local_buf+4; isspace(*ptr); ptr++ ) {}
525 0 : nSPZone = atoi(ptr);
526 :
527 0 : if ( EQUALN(local_buf,"SPCS ",5) )
528 0 : strcpy( local_buf, "SPCS " );
529 0 : else if ( EQUALN(local_buf,"SPAF ",5) )
530 0 : strcpy( local_buf, "SPAF " );
531 : else
532 0 : strcpy( local_buf, "SPIF " );
533 :
534 0 : if( nSPZone != 0 )
535 0 : sprintf( local_buf + 5, "%4d %4s",nSPZone,earthmodel);
536 : else
537 0 : sprintf( local_buf + 5, " %4s",earthmodel);
538 :
539 : }
540 0 : else if( EQUALN(local_buf,"ACEA ",5) )
541 : {
542 : sprintf( local_buf,
543 : "ACEA %4s",
544 0 : earthmodel );
545 : }
546 0 : else if( EQUALN(local_buf,"AE ",3) )
547 : {
548 : sprintf( local_buf,
549 : "AE %4s",
550 0 : earthmodel );
551 : }
552 0 : else if( EQUALN(local_buf,"EC ",3) )
553 : {
554 : sprintf( local_buf,
555 : "EC %4s",
556 0 : earthmodel );
557 : }
558 0 : else if( EQUALN(local_buf,"ER ",3) )
559 : {
560 : sprintf( local_buf,
561 : "ER %4s",
562 0 : earthmodel );
563 : }
564 0 : else if( EQUALN(local_buf,"GNO ",4) )
565 : {
566 : sprintf( local_buf,
567 : "GNO %4s",
568 0 : earthmodel );
569 : }
570 0 : else if( EQUALN(local_buf,"GVNP",4) )
571 : {
572 : sprintf( local_buf,
573 : "GVNP %4s",
574 0 : earthmodel );
575 : }
576 0 : else if( EQUALN(local_buf,"LAEA_ELL",8) )
577 : {
578 : sprintf( local_buf,
579 : "LAEA_ELL %4s",
580 0 : earthmodel );
581 : }
582 0 : else if( EQUALN(local_buf,"LAEA",4) )
583 : {
584 : sprintf( local_buf,
585 : "LAEA %4s",
586 0 : earthmodel );
587 : }
588 0 : else if( EQUALN(local_buf,"LCC_1SP",7) )
589 : {
590 : sprintf( local_buf,
591 : "LCC_1SP %4s",
592 0 : earthmodel );
593 : }
594 0 : else if( EQUALN(local_buf,"LCC ",4) )
595 : {
596 : sprintf( local_buf,
597 : "LCC %4s",
598 0 : earthmodel );
599 : }
600 0 : else if( EQUALN(local_buf,"MC ",3) )
601 : {
602 : sprintf( local_buf,
603 : "MC %4s",
604 0 : earthmodel );
605 : }
606 0 : else if( EQUALN(local_buf,"MER ",4) )
607 : {
608 : sprintf( local_buf,
609 : "MER %4s",
610 0 : earthmodel );
611 : }
612 0 : else if( EQUALN(local_buf,"MSC ",4) )
613 : {
614 : sprintf( local_buf,
615 : "MSC %4s",
616 0 : earthmodel );
617 : }
618 0 : else if( EQUALN(local_buf,"OG ",3) )
619 : {
620 : sprintf( local_buf,
621 : "OG %4s",
622 0 : earthmodel );
623 : }
624 0 : else if( EQUALN(local_buf,"OM ",3) )
625 : {
626 : sprintf( local_buf,
627 : "OM %4s",
628 0 : earthmodel );
629 : }
630 0 : else if( EQUALN(local_buf,"PC ",3) )
631 : {
632 : sprintf( local_buf,
633 : "PC %4s",
634 0 : earthmodel );
635 : }
636 0 : else if( EQUALN(local_buf,"PS ",3) )
637 : {
638 : sprintf( local_buf,
639 : "PS %4s",
640 0 : earthmodel );
641 : }
642 0 : else if( EQUALN(local_buf,"ROB ",4) )
643 : {
644 : sprintf( local_buf,
645 : "ROB %4s",
646 0 : earthmodel );
647 : }
648 0 : else if( EQUALN(local_buf,"SG ",3) )
649 : {
650 : sprintf( local_buf,
651 : "SG %4s",
652 0 : earthmodel );
653 : }
654 0 : else if( EQUALN(local_buf,"SIN ",4) )
655 : {
656 : sprintf( local_buf,
657 : "SIN %4s",
658 0 : earthmodel );
659 : }
660 0 : else if( EQUALN(local_buf,"SOM ",4) )
661 : {
662 : sprintf( local_buf,
663 : "SOM %4s",
664 0 : earthmodel );
665 : }
666 0 : else if( EQUALN(local_buf,"TM ",3) )
667 : {
668 : sprintf( local_buf,
669 : "TM %4s",
670 0 : earthmodel );
671 : }
672 0 : else if( EQUALN(local_buf,"VDG ",4) )
673 : {
674 : sprintf( local_buf,
675 : "VDG %4s",
676 0 : earthmodel );
677 : }
678 0 : else if( EQUALN(local_buf,"UPSA",4) )
679 : {
680 : sprintf( local_buf,
681 : "UPSA %4s",
682 0 : earthmodel );
683 : }
684 0 : else if( EQUALN(local_buf,"UPS ",4) )
685 : {
686 : /* Attempt to find UPS zone */
687 0 : for( ptr=local_buf+3; isspace(*ptr); ptr++ ) {}
688 0 : if( *ptr == 'A' || *ptr == 'B' || *ptr == 'Y' || *ptr == 'Z' )
689 0 : ups_zone = *ptr;
690 0 : else if( *ptr == 'a' || *ptr == 'b' || *ptr == 'y' || *ptr == 'z' )
691 0 : ups_zone = toupper( *ptr );
692 : else
693 0 : ups_zone = ' ';
694 :
695 : sprintf( local_buf,
696 : "UPS %c %4s",
697 0 : ups_zone, earthmodel );
698 : }
699 0 : else if( EQUALN(local_buf,"GOOD",4) )
700 : {
701 : sprintf( local_buf,
702 : "GOOD %4s",
703 0 : earthmodel );
704 : }
705 0 : else if( EQUALN(local_buf,"NZMG",4) )
706 : {
707 : sprintf( local_buf,
708 : "NZMG %4s",
709 0 : earthmodel );
710 : }
711 0 : else if( EQUALN(local_buf,"CASS",4) )
712 : {
713 0 : if( EQUALN( earthmodel, "D000", 4 ) )
714 0 : sprintf( local_buf, "CASS %4s", "E010" );
715 : else
716 0 : sprintf( local_buf, "CASS %4s", earthmodel );
717 : }
718 0 : else if( EQUALN(local_buf,"RSO ",4) )
719 : {
720 0 : if( EQUALN( earthmodel, "D000", 4 ) )
721 0 : sprintf( local_buf, "RSO %4s", "E010" );
722 : else
723 0 : sprintf( local_buf, "RSO %4s", earthmodel );
724 : }
725 0 : else if( EQUALN(local_buf,"KROV",4) )
726 : {
727 0 : if( EQUALN( earthmodel, "D000", 4 ) )
728 0 : sprintf( local_buf, "KROV %4s", "E002" );
729 : else
730 0 : sprintf( local_buf, "KROV %4s", earthmodel );
731 : }
732 0 : else if( EQUALN(local_buf,"KRON",4) )
733 : {
734 0 : if( EQUALN( earthmodel, "D000", 4 ) )
735 0 : sprintf( local_buf, "KRON %4s", "E002" );
736 : else
737 0 : sprintf( local_buf, "KRON %4s", earthmodel );
738 : }
739 0 : else if( EQUALN(local_buf,"SGDO",4) )
740 : {
741 0 : if( EQUALN( earthmodel, "D000", 4 ) )
742 0 : sprintf( local_buf, "SGDO %4s", "E910" );
743 : else
744 0 : sprintf( local_buf, "SGDO %4s", earthmodel );
745 : }
746 0 : else if( EQUALN(local_buf,"LBSG",4) )
747 : {
748 0 : if( EQUALN( earthmodel, "D000", 4 ) )
749 0 : sprintf( local_buf, "LBSG %4s", "E202" );
750 : else
751 0 : sprintf( local_buf, "LBSG %4s", earthmodel );
752 : }
753 0 : else if( EQUALN(local_buf,"ISIN",4) )
754 : {
755 0 : if( EQUALN( earthmodel, "D000", 4 ) )
756 0 : sprintf( local_buf, "ISIN %4s", "E700" );
757 : else
758 0 : sprintf( local_buf, "ISIN %4s", earthmodel );
759 : }
760 : /* -------------------------------------------------------------------- */
761 : /* This may be a user projection. Just reformat the earth model */
762 : /* portion. */
763 : /* -------------------------------------------------------------------- */
764 : else
765 : {
766 0 : sprintf( local_buf, "%-11.11s %4s", geosys.c_str(), earthmodel );
767 : }
768 :
769 298 : return local_buf;
770 : }
771 :
772 : /*
773 : C PAK2PCI converts a Latitude or Longitude value held in decimal
774 : C form to or from the standard packed DMS format DDDMMMSSS.SSS.
775 : C The standard packed DMS format is the required format for any
776 : C Latitude or Longitude value in the projection parameter array
777 : C (TPARIN and/or TPARIO) in calling the U.S.G.S. GCTP package,
778 : C but is not required for the actual coordinates to be converted.
779 : C This routine has been coverted from the PAKPCI fortran routine.
780 : C
781 : C When function is 1, the value returned is made up as follows:
782 : C
783 : C PACK2PCI = (DDD * 1000000) + (MMM * 1000) + SSS.SSS
784 : C
785 : C When function is 0, the value returned is made up as follows:
786 : C
787 : C PACK2PCI = DDD + MMM/60 + SSS/3600
788 : C
789 : C where: DDD are the degrees
790 : C MMM are the minutes
791 : C SSS.SSS are the seconds
792 : C
793 : C The sign of the input value is retained and will denote the
794 : C hemisphere (For longitude, (-) is West and (+) is East of
795 : C Greenwich; For latitude, (-) is South and (+) is North of
796 : C the equator).
797 : C
798 : C
799 : C CALL SEQUENCE
800 : C
801 : C double = PACK2PCI (degrees, function)
802 : C
803 : C degrees - (double) Latitude or Longitude value in decimal
804 : C degrees.
805 : C
806 : C function - (Int) Function to perform
807 : C 1, convert decimal degrees to DDDMMMSSS.SSS
808 : C 0, convert DDDMMMSSS.SSS to decimal degrees
809 : C
810 : C
811 : C EXAMPLE
812 : C
813 : C double degrees, packed, unpack
814 : C
815 : C degrees = -125.425 ! Same as 125d 25' 30" W
816 : C packed = PACK2PCI (degrees, 1) ! PACKED will equal -125025030.000
817 : C unpack = PACK2PCI (degrees, 0) ! UNPACK will equal -125.425
818 : */
819 :
820 : /************************************************************************/
821 : /* PAK2PCI() */
822 : /************************************************************************/
823 :
824 4 : static double PAK2PCI( double deg, int function )
825 : {
826 : double new_deg;
827 : int sign;
828 : double result;
829 :
830 : double degrees;
831 : double temp1, temp2, temp3;
832 : int dd, mm;
833 : double ss;
834 :
835 4 : sign = (int)(1.0);
836 4 : degrees = deg;
837 :
838 4 : if ( degrees < 0 )
839 : {
840 2 : sign = (int)(-1.0);
841 2 : degrees = degrees * sign;
842 : }
843 :
844 : /* -------------------------------------------------------------------- */
845 : /* Unpack the value. */
846 : /* -------------------------------------------------------------------- */
847 4 : if ( function == 0 )
848 : {
849 0 : new_deg = (double) ABS( degrees );
850 :
851 0 : dd = (int)( new_deg / 1000000.0);
852 :
853 0 : new_deg = ( new_deg - (dd * 1000000) );
854 0 : mm = (int)(new_deg/(1000));
855 :
856 0 : new_deg = ( new_deg - (mm * 1000) );
857 :
858 0 : ss = new_deg;
859 :
860 0 : result = (double) sign * ( dd + mm/60.0 + ss/3600.0 );
861 : }
862 : else
863 : {
864 : /* -------------------------------------------------------------------- */
865 : /* Procduce DDDMMMSSS.SSS from decimal degrees. */
866 : /* -------------------------------------------------------------------- */
867 4 : new_deg = (double) ((int)degrees % 360);
868 4 : temp1 = degrees - new_deg;
869 :
870 4 : temp2 = temp1 * 60;
871 :
872 4 : mm = (int)((temp2 * 60) / 60);
873 :
874 4 : temp3 = temp2 - mm;
875 4 : ss = temp3 * 60;
876 :
877 : result = (double) sign *
878 4 : ( (new_deg * 1000000) + (mm * 1000) + ss);
879 : }
880 4 : return result;
881 : }
882 :
883 : /************************************************************************/
884 : /* PrepareGCTPFields() */
885 : /* */
886 : /* Fill the GCTP fields in the seg_data image based on the */
887 : /* non-GCTP values. */
888 : /************************************************************************/
889 :
890 170 : void CPCIDSKGeoref::PrepareGCTPFields()
891 :
892 : {
893 : enum GCTP_UNIT_CODES {
894 : GCTP_UNIT_UNKNOWN = -1, /* Default, NOT a valid code */
895 : GCTP_UNIT_RADIAN = 0, /* 0, NOT used at present */
896 : GCTP_UNIT_US_FOOT, /* 1, Used for GEO_SPAF */
897 : GCTP_UNIT_METRE, /* 2, Used for most map projections */
898 : GCTP_UNIT_SECOND, /* 3, NOT used at present */
899 : GCTP_UNIT_DEGREE, /* 4, Used for GEO_LONG */
900 : GCTP_UNIT_INTL_FOOT, /* 5, Used for GEO_SPIF */
901 : GCTP_UNIT_TABLE /* 6, NOT used at present */
902 : };
903 :
904 170 : seg_data.Get(32,16,geosys);
905 170 : std::string geosys_clean(ReformatGeosys( geosys ));
906 :
907 : /* -------------------------------------------------------------------- */
908 : /* Establish the GCTP units code. */
909 : /* -------------------------------------------------------------------- */
910 170 : double IOmultiply = 1.0;
911 170 : int UnitsCode = GCTP_UNIT_METRE;
912 :
913 170 : std::string grid_units;
914 170 : seg_data.Get(64,16,grid_units);
915 :
916 170 : if( EQUALN(grid_units.c_str(),"MET",3) )
917 88 : UnitsCode = GCTP_UNIT_METRE;
918 82 : else if( EQUALN( grid_units.c_str(), "FOOT", 4 ) )
919 : {
920 0 : UnitsCode = GCTP_UNIT_US_FOOT;
921 0 : IOmultiply = 1.0 / 0.3048006096012192;
922 : }
923 82 : else if( EQUALN( grid_units.c_str(), "INTL FOOT", 9 ) )
924 : {
925 0 : UnitsCode = GCTP_UNIT_INTL_FOOT;
926 0 : IOmultiply = 1.0 / 0.3048;
927 : }
928 82 : else if( EQUALN( grid_units.c_str(), "DEGREE", 6 ) )
929 41 : UnitsCode = GCTP_UNIT_DEGREE;
930 :
931 : /* -------------------------------------------------------------------- */
932 : /* Extract the non-GCTP style parameters. */
933 : /* -------------------------------------------------------------------- */
934 : double pci_parms[17];
935 : int i;
936 :
937 3060 : for( i = 0; i < 17; i++ )
938 2890 : pci_parms[i] = seg_data.GetDouble(80+26*i,26);
939 :
940 : #define Dearth0 pci_parms[0]
941 : #define Dearth1 pci_parms[1]
942 : #define RefLong pci_parms[2]
943 : #define RefLat pci_parms[3]
944 : #define StdParallel1 pci_parms[4]
945 : #define StdParallel2 pci_parms[5]
946 : #define FalseEasting pci_parms[6]
947 : #define FalseNorthing pci_parms[7]
948 : #define Scale pci_parms[8]
949 : #define Height pci_parms[9]
950 : #define Long1 pci_parms[10]
951 : #define Lat1 pci_parms[11]
952 : #define Long2 pci_parms[12]
953 : #define Lat2 pci_parms[13]
954 : #define Azimuth pci_parms[14]
955 : #define LandsatNum pci_parms[15]
956 : #define LandsatPath pci_parms[16]
957 :
958 : /* -------------------------------------------------------------------- */
959 : /* Get the zone code. */
960 : /* -------------------------------------------------------------------- */
961 170 : int ProjectionZone = 0;
962 :
963 170 : if( strncmp(geosys_clean.c_str(),"UTM ",4) == 0
964 : || strncmp(geosys_clean.c_str(),"SPCS ",5) == 0
965 : || strncmp(geosys_clean.c_str(),"SPAF ",5) == 0
966 : || strncmp(geosys_clean.c_str(),"SPIF ",5) == 0 )
967 : {
968 2 : ProjectionZone = atoi(geosys_clean.c_str() + 5);
969 : }
970 :
971 : /* -------------------------------------------------------------------- */
972 : /* Handle the ellipsoid. We depend on applications properly */
973 : /* setting proj_parms[0], and proj_parms[1] with the semi-major */
974 : /* and semi-minor axes in all other cases. */
975 : /* */
976 : /* I wish we could lookup datum codes to find their GCTP */
977 : /* ellipsoid values here! */
978 : /* -------------------------------------------------------------------- */
979 170 : int Spheroid = -1;
980 170 : if( geosys_clean[12] == 'E' )
981 2 : Spheroid = atoi(geosys_clean.c_str() + 13);
982 :
983 170 : if( Spheroid < 0 || Spheroid > 19 )
984 168 : Spheroid = -1;
985 :
986 : /* -------------------------------------------------------------------- */
987 : /* Initialize the USGS Parameters. */
988 : /* -------------------------------------------------------------------- */
989 : double USGSParms[15];
990 : int gsys;
991 :
992 2720 : for ( i = 0; i < 15; i++ )
993 2550 : USGSParms[i] = 0;
994 :
995 : /* -------------------------------------------------------------------- */
996 : /* Projection 0: Geographic (no projection) */
997 : /* -------------------------------------------------------------------- */
998 170 : if( strncmp(geosys_clean.c_str(),"LONG ",5) == 0 )
999 : {
1000 0 : gsys = 0;
1001 0 : UnitsCode = GCTP_UNIT_DEGREE;
1002 : }
1003 :
1004 : /* -------------------------------------------------------------------- */
1005 : /* Projection 1: Universal Transverse Mercator */
1006 : /* -------------------------------------------------------------------- */
1007 170 : else if( strncmp(geosys_clean.c_str(),"UTM ",4) == 0 )
1008 : {
1009 2 : char row_char = geosys_clean[10];
1010 2 : gsys = 1;
1011 :
1012 : // Southern hemisphere?
1013 2 : if( (row_char >= 'C') && (row_char <= 'M') && ProjectionZone > 0 )
1014 : {
1015 0 : ProjectionZone *= -1;
1016 : }
1017 :
1018 : /* -------------------------------------------------------------------- */
1019 : /* Process UTM as TM. The reason for this is the GCTP software */
1020 : /* does not provide for input of an Earth Model for UTM, but does */
1021 : /* for TM. */
1022 : /* -------------------------------------------------------------------- */
1023 2 : gsys = 9;
1024 2 : USGSParms[0] = Dearth0;
1025 2 : USGSParms[1] = Dearth1;
1026 2 : USGSParms[2] = 0.9996;
1027 :
1028 : USGSParms[4] = PAK2PCI(
1029 2 : ( ABS(ProjectionZone) * 6.0 - 183.0 ), 1 );
1030 2 : USGSParms[5] = PAK2PCI( 0.0, 1 );
1031 2 : USGSParms[6] = 500000.0;
1032 2 : USGSParms[7] = ( ProjectionZone < 0 ) ? 10000000.0 : 0.0;
1033 : }
1034 :
1035 : /* -------------------------------------------------------------------- */
1036 : /* Projection 2: State Plane Coordinate System */
1037 : /* -------------------------------------------------------------------- */
1038 168 : else if( strncmp(geosys_clean.c_str(),"SPCS ",5) == 0 )
1039 : {
1040 0 : gsys = 2;
1041 0 : if( UnitsCode != GCTP_UNIT_METRE
1042 : && UnitsCode != GCTP_UNIT_US_FOOT
1043 : && UnitsCode != GCTP_UNIT_INTL_FOOT )
1044 0 : UnitsCode = GCTP_UNIT_METRE;
1045 : }
1046 :
1047 168 : else if( strncmp(geosys_clean.c_str(),"SPAF ",5) == 0 )
1048 : {
1049 0 : gsys = 2;
1050 0 : if( UnitsCode != GCTP_UNIT_METRE
1051 : && UnitsCode != GCTP_UNIT_US_FOOT
1052 : && UnitsCode != GCTP_UNIT_INTL_FOOT )
1053 0 : UnitsCode = GCTP_UNIT_US_FOOT;
1054 : }
1055 :
1056 168 : else if( strncmp(geosys_clean.c_str(),"SPIF ",5) == 0 )
1057 : {
1058 0 : gsys = 2;
1059 0 : if( UnitsCode != GCTP_UNIT_METRE
1060 : && UnitsCode != GCTP_UNIT_US_FOOT
1061 : && UnitsCode != GCTP_UNIT_INTL_FOOT )
1062 0 : UnitsCode = GCTP_UNIT_INTL_FOOT;
1063 : }
1064 :
1065 : /* -------------------------------------------------------------------- */
1066 : /* Projection 3: Albers Conical Equal-Area */
1067 : /* -------------------------------------------------------------------- */
1068 168 : else if( strncmp(geosys_clean.c_str(),"ACEA ",5) == 0 )
1069 : {
1070 0 : gsys = 3;
1071 0 : USGSParms[0] = Dearth0;
1072 0 : USGSParms[1] = Dearth1;
1073 0 : USGSParms[2] = PAK2PCI(StdParallel1, 1);
1074 0 : USGSParms[3] = PAK2PCI(StdParallel2, 1);
1075 0 : USGSParms[4] = PAK2PCI(RefLong, 1);
1076 0 : USGSParms[5] = PAK2PCI(RefLat, 1);
1077 0 : USGSParms[6] = FalseEasting * IOmultiply;
1078 0 : USGSParms[7] = FalseNorthing * IOmultiply;
1079 : }
1080 :
1081 : /* -------------------------------------------------------------------- */
1082 : /* Projection 4: Lambert Conformal Conic */
1083 : /* -------------------------------------------------------------------- */
1084 168 : else if( strncmp(geosys_clean.c_str(),"LCC ",5) == 0 )
1085 : {
1086 0 : gsys = 4;
1087 0 : USGSParms[0] = Dearth0;
1088 0 : USGSParms[1] = Dearth1;
1089 0 : USGSParms[2] = PAK2PCI(StdParallel1, 1);
1090 0 : USGSParms[3] = PAK2PCI(StdParallel2, 1);
1091 0 : USGSParms[4] = PAK2PCI(RefLong, 1);
1092 0 : USGSParms[5] = PAK2PCI(RefLat, 1);
1093 0 : USGSParms[6] = FalseEasting * IOmultiply;
1094 0 : USGSParms[7] = FalseNorthing * IOmultiply;
1095 : }
1096 :
1097 : /* -------------------------------------------------------------------- */
1098 : /* Projection 5: Mercator */
1099 : /* -------------------------------------------------------------------- */
1100 168 : else if( strncmp(geosys_clean.c_str(),"MER ",5) == 0 )
1101 : {
1102 0 : gsys = 5;
1103 0 : USGSParms[0] = Dearth0;
1104 0 : USGSParms[1] = Dearth1;
1105 :
1106 0 : USGSParms[4] = PAK2PCI(RefLong, 1);
1107 0 : USGSParms[5] = PAK2PCI(RefLat, 1);
1108 0 : USGSParms[6] = FalseEasting * IOmultiply;
1109 0 : USGSParms[7] = FalseNorthing * IOmultiply;
1110 : }
1111 :
1112 : /* -------------------------------------------------------------------- */
1113 : /* Projection 6: Polar Stereographic */
1114 : /* -------------------------------------------------------------------- */
1115 168 : else if( strncmp(geosys_clean.c_str(),"PS ",5) == 0 )
1116 : {
1117 0 : gsys = 6;
1118 0 : USGSParms[0] = Dearth0;
1119 0 : USGSParms[1] = Dearth1;
1120 :
1121 0 : USGSParms[4] = PAK2PCI(RefLong, 1);
1122 0 : USGSParms[5] = PAK2PCI(RefLat, 1);
1123 0 : USGSParms[6] = FalseEasting * IOmultiply;
1124 0 : USGSParms[7] = FalseNorthing * IOmultiply;
1125 : }
1126 :
1127 : /* -------------------------------------------------------------------- */
1128 : /* Projection 7: Polyconic */
1129 : /* -------------------------------------------------------------------- */
1130 168 : else if( strncmp(geosys_clean.c_str(),"PC ",5) == 0 )
1131 : {
1132 0 : gsys = 7;
1133 0 : USGSParms[0] = Dearth0;
1134 0 : USGSParms[1] = Dearth1;
1135 :
1136 0 : USGSParms[4] = PAK2PCI(RefLong, 1);
1137 0 : USGSParms[5] = PAK2PCI(RefLat, 1);
1138 0 : USGSParms[6] = FalseEasting * IOmultiply;
1139 0 : USGSParms[7] = FalseNorthing * IOmultiply;
1140 : }
1141 :
1142 : /* -------------------------------------------------------------------- */
1143 : /* Projection 8: Equidistant Conic */
1144 : /* Format A, one standard parallel, usgs_params[8] = 0 */
1145 : /* Format B, two standard parallels, usgs_params[8] = not 0 */
1146 : /* -------------------------------------------------------------------- */
1147 168 : else if( strncmp(geosys_clean.c_str(),"EC ",5) == 0 )
1148 : {
1149 0 : gsys = 8;
1150 0 : USGSParms[0] = Dearth0;
1151 0 : USGSParms[1] = Dearth1;
1152 0 : USGSParms[2] = PAK2PCI(StdParallel1, 1);
1153 0 : USGSParms[3] = PAK2PCI(StdParallel2, 1);
1154 0 : USGSParms[4] = PAK2PCI(RefLong, 1);
1155 0 : USGSParms[5] = PAK2PCI(RefLat, 1);
1156 0 : USGSParms[6] = FalseEasting * IOmultiply;
1157 0 : USGSParms[7] = FalseNorthing * IOmultiply;
1158 :
1159 0 : if ( StdParallel2 != 0 )
1160 : {
1161 0 : USGSParms[8] = 1;
1162 : }
1163 : }
1164 :
1165 : /* -------------------------------------------------------------------- */
1166 : /* Projection 9: Transverse Mercator */
1167 : /* -------------------------------------------------------------------- */
1168 168 : else if( strncmp(geosys_clean.c_str(),"TM ",5) == 0 )
1169 : {
1170 0 : gsys = 9;
1171 0 : USGSParms[0] = Dearth0;
1172 0 : USGSParms[1] = Dearth1;
1173 0 : USGSParms[2] = Scale;
1174 :
1175 0 : USGSParms[4] = PAK2PCI(RefLong, 1);
1176 0 : USGSParms[5] = PAK2PCI(RefLat, 1);
1177 0 : USGSParms[6] = FalseEasting * IOmultiply;
1178 0 : USGSParms[7] = FalseNorthing * IOmultiply;
1179 : }
1180 :
1181 : /* -------------------------------------------------------------------- */
1182 : /* Projection 10: Stereographic */
1183 : /* -------------------------------------------------------------------- */
1184 168 : else if( strncmp(geosys_clean.c_str(),"SG ",5) == 0 )
1185 : {
1186 0 : gsys = 10;
1187 0 : USGSParms[0] = Dearth0;
1188 :
1189 0 : USGSParms[4] = PAK2PCI(RefLong, 1);
1190 0 : USGSParms[5] = PAK2PCI(RefLat, 1);
1191 0 : USGSParms[6] = FalseEasting * IOmultiply;
1192 0 : USGSParms[7] = FalseNorthing * IOmultiply;
1193 : }
1194 :
1195 : /* -------------------------------------------------------------------- */
1196 : /* Projection 11: Lambert Azimuthal Equal-Area */
1197 : /* -------------------------------------------------------------------- */
1198 168 : else if( strncmp(geosys_clean.c_str(),"LAEA ",5) == 0 )
1199 : {
1200 0 : gsys = 11;
1201 :
1202 0 : USGSParms[0] = Dearth0;
1203 :
1204 0 : USGSParms[4] = PAK2PCI(RefLong, 1);
1205 0 : USGSParms[5] = PAK2PCI(RefLat, 1);
1206 0 : USGSParms[6] = FalseEasting * IOmultiply;
1207 0 : USGSParms[7] = FalseNorthing * IOmultiply;
1208 : }
1209 :
1210 : /* -------------------------------------------------------------------- */
1211 : /* Projection 12: Azimuthal Equidistant */
1212 : /* -------------------------------------------------------------------- */
1213 168 : else if( strncmp(geosys_clean.c_str(),"AE ",5) == 0 )
1214 : {
1215 0 : gsys = 12;
1216 0 : USGSParms[0] = Dearth0;
1217 :
1218 0 : USGSParms[4] = PAK2PCI(RefLong, 1);
1219 0 : USGSParms[5] = PAK2PCI(RefLat, 1);
1220 0 : USGSParms[6] = FalseEasting * IOmultiply;
1221 0 : USGSParms[7] = FalseNorthing * IOmultiply;
1222 : }
1223 :
1224 : /* -------------------------------------------------------------------- */
1225 : /* Projection 13: Gnomonic */
1226 : /* -------------------------------------------------------------------- */
1227 168 : else if( strncmp(geosys_clean.c_str(),"GNO ",5) == 0 )
1228 : {
1229 0 : gsys = 13;
1230 0 : USGSParms[0] = Dearth0;
1231 :
1232 0 : USGSParms[4] = PAK2PCI(RefLong, 1);
1233 0 : USGSParms[5] = PAK2PCI(RefLat, 1);
1234 0 : USGSParms[6] = FalseEasting * IOmultiply;
1235 0 : USGSParms[7] = FalseNorthing * IOmultiply;
1236 : }
1237 :
1238 : /* -------------------------------------------------------------------- */
1239 : /* Projection 14: Orthographic */
1240 : /* -------------------------------------------------------------------- */
1241 168 : else if( strncmp(geosys_clean.c_str(),"OG ",5) == 0 )
1242 : {
1243 0 : gsys = 14;
1244 0 : USGSParms[0] = Dearth0;
1245 :
1246 0 : USGSParms[4] = PAK2PCI(RefLong, 1);
1247 0 : USGSParms[5] = PAK2PCI(RefLat, 1);
1248 0 : USGSParms[6] = FalseEasting * IOmultiply;
1249 0 : USGSParms[7] = FalseNorthing * IOmultiply;
1250 : }
1251 :
1252 : /* -------------------------------------------------------------------- */
1253 : /* Projection 15: General Vertical Near-Side Perspective */
1254 : /* -------------------------------------------------------------------- */
1255 168 : else if( strncmp(geosys_clean.c_str(),"GVNP ",5) == 0 )
1256 : {
1257 0 : gsys = 15;
1258 0 : USGSParms[0] = Dearth0;
1259 :
1260 0 : USGSParms[2] = Height;
1261 :
1262 0 : USGSParms[4] = PAK2PCI(RefLong, 1);
1263 0 : USGSParms[5] = PAK2PCI(RefLat, 1);
1264 0 : USGSParms[6] = FalseEasting * IOmultiply;
1265 0 : USGSParms[7] = FalseNorthing * IOmultiply;
1266 : }
1267 :
1268 : /* -------------------------------------------------------------------- */
1269 : /* Projection 16: Sinusoidal */
1270 : /* -------------------------------------------------------------------- */
1271 168 : else if( strncmp(geosys_clean.c_str(),"SIN ",5) == 0 )
1272 : {
1273 0 : gsys = 16;
1274 0 : USGSParms[0] = Dearth0;
1275 0 : USGSParms[4] = PAK2PCI(RefLong, 1);
1276 0 : USGSParms[6] = FalseEasting * IOmultiply;
1277 0 : USGSParms[7] = FalseNorthing * IOmultiply;
1278 : }
1279 :
1280 : /* -------------------------------------------------------------------- */
1281 : /* Projection 17: Equirectangular */
1282 : /* -------------------------------------------------------------------- */
1283 168 : else if( strncmp(geosys_clean.c_str(),"ER ",5) == 0 )
1284 : {
1285 0 : gsys = 17;
1286 0 : USGSParms[0] = Dearth0;
1287 0 : USGSParms[4] = PAK2PCI(RefLong, 1);
1288 0 : USGSParms[5] = PAK2PCI(RefLat, 1);
1289 0 : USGSParms[6] = FalseEasting * IOmultiply;
1290 0 : USGSParms[7] = FalseNorthing * IOmultiply;
1291 : }
1292 : /* -------------------------------------------------------------------- */
1293 : /* Projection 18: Miller Cylindrical */
1294 : /* -------------------------------------------------------------------- */
1295 168 : else if( strncmp(geosys_clean.c_str(),"MC ",5) == 0 )
1296 : {
1297 0 : gsys = 18;
1298 0 : USGSParms[0] = Dearth0;
1299 :
1300 0 : USGSParms[4] = PAK2PCI(RefLong, 1);
1301 :
1302 0 : USGSParms[6] = FalseEasting * IOmultiply;
1303 0 : USGSParms[7] = FalseNorthing * IOmultiply;
1304 : }
1305 :
1306 : /* -------------------------------------------------------------------- */
1307 : /* Projection 19: Van der Grinten */
1308 : /* -------------------------------------------------------------------- */
1309 168 : else if( strncmp(geosys_clean.c_str(),"VDG ",5) == 0 )
1310 : {
1311 0 : gsys = 19;
1312 0 : USGSParms[0] = Dearth0;
1313 :
1314 0 : USGSParms[4] = PAK2PCI(RefLong, 1);
1315 :
1316 0 : USGSParms[6] = FalseEasting * IOmultiply;
1317 0 : USGSParms[7] = FalseNorthing * IOmultiply;
1318 : }
1319 :
1320 : /* -------------------------------------------------------------------- */
1321 : /* Projection 20: Oblique Mercator (Hotine) */
1322 : /* Format A, Azimuth and RefLong defined (Long1, Lat1, */
1323 : /* Long2, Lat2 not defined), usgs_params[12] = 0 */
1324 : /* Format B, Long1, Lat1, Long2, Lat2 defined (Azimuth */
1325 : /* and RefLong not defined), usgs_params[12] = not 0 */
1326 : /* -------------------------------------------------------------------- */
1327 168 : else if( strncmp(geosys_clean.c_str(),"OM ",5) == 0 )
1328 : {
1329 0 : gsys = 20;
1330 0 : USGSParms[0] = Dearth0;
1331 0 : USGSParms[1] = Dearth1;
1332 0 : USGSParms[2] = Scale;
1333 0 : USGSParms[3] = PAK2PCI(Azimuth ,1);
1334 :
1335 0 : USGSParms[4] = PAK2PCI(RefLong, 1);
1336 0 : USGSParms[5] = PAK2PCI(RefLat, 1);
1337 0 : USGSParms[6] = FalseEasting * IOmultiply;
1338 0 : USGSParms[7] = FalseNorthing * IOmultiply;
1339 :
1340 0 : USGSParms[8] = PAK2PCI(Long1, 1);
1341 0 : USGSParms[9] = PAK2PCI(Lat1, 1);
1342 0 : USGSParms[10] = PAK2PCI(Long2, 1);
1343 0 : USGSParms[11] = PAK2PCI(Lat2, 1);
1344 0 : if ( (Long1 != 0) || (Lat1 != 0) ||
1345 : (Long2 != 0) || (Lat2 != 0) )
1346 0 : USGSParms[12] = 0.0;
1347 : else
1348 0 : USGSParms[12] = 1.0;
1349 : }
1350 : /* -------------------------------------------------------------------- */
1351 : /* Projection 21: Robinson */
1352 : /* -------------------------------------------------------------------- */
1353 168 : else if( strncmp(geosys_clean.c_str(),"ROB ",5) == 0 )
1354 : {
1355 0 : gsys = 21;
1356 0 : USGSParms[0] = Dearth0;
1357 :
1358 0 : USGSParms[4] = PAK2PCI(RefLong, 1);
1359 : USGSParms[6] = FalseEasting
1360 0 : * IOmultiply;
1361 : USGSParms[7] = FalseNorthing
1362 0 : * IOmultiply;
1363 :
1364 : }
1365 : /* -------------------------------------------------------------------- */
1366 : /* Projection 22: Space Oblique Mercator */
1367 : /* -------------------------------------------------------------------- */
1368 168 : else if( strncmp(geosys_clean.c_str(),"SOM ",5) == 0 )
1369 : {
1370 0 : gsys = 22;
1371 0 : USGSParms[0] = Dearth0;
1372 0 : USGSParms[1] = Dearth1;
1373 0 : USGSParms[2] = LandsatNum;
1374 0 : USGSParms[3] = LandsatPath;
1375 : USGSParms[6] = FalseEasting
1376 0 : * IOmultiply;
1377 : USGSParms[7] = FalseNorthing
1378 0 : * IOmultiply;
1379 : }
1380 : /* -------------------------------------------------------------------- */
1381 : /* Projection 23: Modified Stereographic Conformal (Alaska) */
1382 : /* -------------------------------------------------------------------- */
1383 168 : else if( strncmp(geosys_clean.c_str(),"MSC ",5) == 0 )
1384 : {
1385 0 : gsys = 23;
1386 0 : USGSParms[0] = Dearth0;
1387 0 : USGSParms[1] = Dearth1;
1388 : USGSParms[6] = FalseEasting
1389 0 : * IOmultiply;
1390 : USGSParms[7] = FalseNorthing
1391 0 : * IOmultiply;
1392 : }
1393 :
1394 : /* -------------------------------------------------------------------- */
1395 : /* Projection 6: Universal Polar Stereographic is just Polar */
1396 : /* Stereographic with certain assumptions. */
1397 : /* -------------------------------------------------------------------- */
1398 168 : else if( strncmp(geosys_clean.c_str(),"UPS ",5) == 0 )
1399 : {
1400 0 : gsys = 6;
1401 :
1402 0 : USGSParms[0] = Dearth0;
1403 0 : USGSParms[1] = Dearth1;
1404 :
1405 0 : USGSParms[4] = PAK2PCI(0.0, 1);
1406 :
1407 0 : USGSParms[6] = 2000000.0;
1408 0 : USGSParms[7] = 2000000.0;
1409 :
1410 0 : double dwork = 81.0 + 6.0/60.0 + 52.3/3600.0;
1411 :
1412 0 : if( geosys_clean[10] == 'A' || geosys_clean[10] == 'B' )
1413 : {
1414 0 : USGSParms[5] = PAK2PCI(-dwork,1);
1415 : }
1416 0 : else if( geosys_clean[10] == 'Y' || geosys_clean[10]=='Z')
1417 : {
1418 0 : USGSParms[5] = PAK2PCI(dwork,1);
1419 : }
1420 : else
1421 : {
1422 0 : USGSParms[4] = PAK2PCI(RefLong, 1);
1423 0 : USGSParms[5] = PAK2PCI(RefLat, 1);
1424 : USGSParms[6] = FalseEasting
1425 0 : * IOmultiply;
1426 : USGSParms[7] = FalseNorthing
1427 0 : * IOmultiply;
1428 : }
1429 : }
1430 :
1431 : /* -------------------------------------------------------------------- */
1432 : /* Unknown code. */
1433 : /* -------------------------------------------------------------------- */
1434 : else
1435 : {
1436 168 : gsys = -1;
1437 : }
1438 :
1439 170 : if( ProjectionZone == 0 )
1440 168 : ProjectionZone = 10000 + gsys;
1441 :
1442 : /* -------------------------------------------------------------------- */
1443 : /* Place USGS values in the formatted segment. */
1444 : /* -------------------------------------------------------------------- */
1445 170 : seg_data.Put( (double) gsys, 1458 , 26, "%26.18lE" );
1446 170 : seg_data.Put( (double) ProjectionZone, 1458+26, 26, "%26.18lE" );
1447 :
1448 2720 : for( i = 0; i < 15; i++ )
1449 2550 : seg_data.Put( USGSParms[i], 1458+26*(2+i), 26, "%26.18lE" );
1450 :
1451 170 : seg_data.Put( (double) UnitsCode, 1458+26*17, 26, "%26.18lE" );
1452 170 : seg_data.Put( (double) Spheroid, 1458+26*18, 26, "%26.18lE" );
1453 170 : }
1454 :
1455 :
|
