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