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