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