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