Coverage Report

Created: 2025-11-15 08:43

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