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