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