Coverage Report

Created: 2026-08-11 08:26

next uncovered line (L), next uncovered region (R), next uncovered branch (B)
/src/gdal/frmts/gtiff/libgeotiff/geo_normalize.c
Line
Count
Source
1
/******************************************************************************
2
 *
3
 * Project:  libgeotiff
4
 * Purpose:  Code to normalize PCS and other composite codes in a GeoTIFF file.
5
 * Author:   Frank Warmerdam, warmerdam@pobox.com
6
 *
7
 ******************************************************************************
8
 * Copyright (c) 1999, Frank Warmerdam
9
 * Copyright (c) 2018, Even Rouault <even.rouault at spatialys.com>
10
 *
11
 * SPDX-License-Identifier: MIT
12
 *****************************************************************************/
13
14
#include <assert.h>
15
#include <stdio.h>
16
#include <stdlib.h>
17
#include <string.h>
18
19
#include "cpl_serv.h"
20
#include "geo_tiffp.h"
21
#include "geovalues.h"
22
#include "geo_normalize.h"
23
#include "geo_keyp.h"
24
25
#include "proj.h"
26
27
#ifndef KvUserDefined
28
#  define KvUserDefined 32767
29
#endif
30
31
#ifndef M_PI
32
#  define M_PI 3.14159265358979323846
33
#endif
34
35
/* EPSG Codes for projection parameters.  Unfortunately, these bear no
36
   relationship to the GeoTIFF codes even though the names are so similar. */
37
38
21.1k
#define EPSGNatOriginLat         8801
39
21.0k
#define EPSGNatOriginLong        8802
40
141k
#define EPSGNatOriginScaleFactor 8805
41
21.0k
#define EPSGFalseEasting         8806
42
21.0k
#define EPSGFalseNorthing        8807
43
266
#define EPSGProjCenterLat        8811
44
266
#define EPSGProjCenterLong       8812
45
266
#define EPSGAzimuth              8813
46
61.0k
#define EPSGAngleRectifiedToSkewedGrid 8814
47
115k
#define EPSGInitialLineScaleFactor 8815
48
324
#define EPSGProjCenterEasting    8816
49
291
#define EPSGProjCenterNorthing   8817
50
#define EPSGPseudoStdParallelLat 8818
51
54.7k
#define EPSGPseudoStdParallelScaleFactor 8819
52
4.40k
#define EPSGFalseOriginLat       8821
53
4.40k
#define EPSGFalseOriginLong      8822
54
4.55k
#define EPSGStdParallel1Lat      8823
55
4.40k
#define EPSGStdParallel2Lat      8824
56
4.40k
#define EPSGFalseOriginEasting   8826
57
4.40k
#define EPSGFalseOriginNorthing  8827
58
#define EPSGSphericalOriginLat   8828
59
#define EPSGSphericalOriginLong  8829
60
#define EPSGInitialLongitude     8830
61
#define EPSGZoneWidth            8831
62
54
#define EPSGLatOfStdParallel     8832
63
54
#define EPSGOriginLong           8833
64
#define EPSGTopocentricOriginLat 8834
65
#define EPSGTopocentricOriginLong 8835
66
#define EPSGTopocentricOriginHeight 8836
67
68
0
#define CT_Ext_Mercator_2SP     -CT_Mercator
69
70
#ifndef CPL_INLINE
71
#  if (defined(__GNUC__) && !defined(__NO_INLINE__)) || defined(_MSC_VER)
72
#    define HAS_CPL_INLINE  1
73
#    define CPL_INLINE __inline
74
#  elif defined(__SUNPRO_CC)
75
#    define HAS_CPL_INLINE  1
76
#    define CPL_INLINE inline
77
#  else
78
#    define CPL_INLINE
79
#  endif
80
#endif
81
82
#ifndef CPL_UNUSED
83
#if defined(__GNUC__) && __GNUC__ >= 4
84
#  define CPL_UNUSED __attribute((__unused__))
85
#else
86
/* TODO: add cases for other compilers */
87
#  define CPL_UNUSED
88
#endif
89
#endif
90
91
166k
CPL_INLINE static void CPL_IGNORE_RET_VAL_INT(CPL_UNUSED int unused) {}
92
93
/************************************************************************/
94
/*                          GTIFKeyGetSSHORT()                          */
95
/************************************************************************/
96
97
// Geotiff SHORT keys are supposed to be unsigned, but geo_normalize interface
98
// uses signed short...
99
static int GTIFKeyGetSSHORT( GTIF *gtif, geokey_t key, short* pnVal )
100
592k
{
101
592k
    unsigned short sVal;
102
592k
    if( GTIFKeyGetSHORT(gtif, key, &sVal, 0, 1) == 1 )
103
168k
    {
104
168k
        memcpy(pnVal, &sVal, 2);
105
168k
        return 1;
106
168k
    }
107
424k
    return 0;
108
592k
}
109
110
/************************************************************************/
111
/*                           GTIFGetPCSInfo()                           */
112
/************************************************************************/
113
114
int GTIFGetPCSInfoEx( void* ctxIn,
115
                      int nPCSCode, char **ppszEPSGName,
116
                      short *pnProjOp, short *pnUOMLengthCode,
117
                      short *pnGeogCS )
118
119
68.7k
{
120
68.7k
    int         nDatum;
121
68.7k
    int         nZone;
122
123
    /* Deal with a few well known CRS */
124
68.7k
    int Proj = GTIFPCSToMapSys( nPCSCode, &nDatum, &nZone );
125
68.7k
    if ((Proj == MapSys_UTM_North || Proj == MapSys_UTM_South) &&
126
16.0k
        nDatum != KvUserDefined)
127
16.0k
    {
128
16.0k
        const char* pszDatumName = NULL;
129
16.0k
        switch (nDatum)
130
16.0k
        {
131
318
            case GCS_NAD27: pszDatumName = "NAD27"; break;
132
1.28k
            case GCS_NAD83: pszDatumName = "NAD83"; break;
133
1.19k
            case GCS_WGS_72: pszDatumName = "WGS 72"; break;
134
12.3k
            case GCS_WGS_72BE: pszDatumName = "WGS 72BE"; break;
135
842
            case GCS_WGS_84: pszDatumName = "WGS 84"; break;
136
0
            default: break;
137
16.0k
        }
138
139
16.0k
        if (pszDatumName)
140
16.0k
        {
141
16.0k
            if (ppszEPSGName)
142
7.94k
            {
143
7.94k
                char szEPSGName[64];
144
7.94k
                snprintf(
145
7.94k
                    szEPSGName, sizeof(szEPSGName), "%s / UTM zone %d%c",
146
7.94k
                    pszDatumName, nZone,
147
7.94k
                    (Proj == MapSys_UTM_North) ? 'N' : 'S');
148
7.94k
                *ppszEPSGName = CPLStrdup(szEPSGName);
149
7.94k
            }
150
151
16.0k
            if (pnProjOp)
152
8.08k
                *pnProjOp = (short) (((Proj == MapSys_UTM_North) ? Proj_UTM_zone_1N - 1 : Proj_UTM_zone_1S - 1) + nZone);
153
154
16.0k
            if (pnUOMLengthCode)
155
8.08k
                *pnUOMLengthCode = 9001; /* Linear_Meter */
156
157
16.0k
            if (pnGeogCS)
158
8.08k
                *pnGeogCS = (short) nDatum;
159
160
16.0k
            return TRUE;
161
16.0k
        }
162
16.0k
    }
163
164
52.7k
    if( nPCSCode == KvUserDefined )
165
0
        return FALSE;
166
167
52.7k
    PJ_CONTEXT* ctx = (PJ_CONTEXT*)ctxIn;
168
52.7k
    {
169
52.7k
        char szCode[12];
170
171
52.7k
        snprintf(szCode, sizeof(szCode), "%d", nPCSCode);
172
52.7k
        PJ* proj_crs = proj_create_from_database(
173
52.7k
            ctx, "EPSG", szCode, PJ_CATEGORY_CRS, 0, NULL);
174
52.7k
        if( !proj_crs )
175
33.5k
        {
176
33.5k
            return FALSE;
177
33.5k
        }
178
179
19.1k
        if( proj_get_type(proj_crs) != PJ_TYPE_PROJECTED_CRS )
180
2.65k
        {
181
2.65k
            proj_destroy(proj_crs);
182
2.65k
            return FALSE;
183
2.65k
        }
184
185
16.4k
        if( ppszEPSGName )
186
8.16k
        {
187
8.16k
            const char* pszName = proj_get_name(proj_crs);
188
8.16k
            if( !pszName )
189
0
            {
190
                // shouldn't happen
191
0
                proj_destroy(proj_crs);
192
0
                return FALSE;
193
0
            }
194
8.16k
            *ppszEPSGName = CPLStrdup(pszName);
195
8.16k
        }
196
197
16.4k
        if( pnProjOp )
198
8.33k
        {
199
8.33k
            PJ* conversion = proj_crs_get_coordoperation(
200
8.33k
                ctx, proj_crs);
201
8.33k
            if( !conversion )
202
0
            {
203
                // shouldn't happen except out of memory
204
0
                proj_destroy(proj_crs);
205
0
                return FALSE;
206
0
            }
207
208
8.33k
            {
209
8.33k
                const char* pszConvCode = proj_get_id_code(conversion, 0);
210
8.33k
                assert( pszConvCode );
211
8.33k
                *pnProjOp = (short) atoi(pszConvCode);
212
8.33k
            }
213
214
8.33k
            proj_destroy(conversion);
215
8.33k
        }
216
217
16.4k
        if( pnUOMLengthCode )
218
8.33k
        {
219
8.33k
            PJ* coordSys = proj_crs_get_coordinate_system(
220
8.33k
                ctx, proj_crs);
221
8.33k
            if( !coordSys )
222
0
            {
223
                // shouldn't happen except out of memory
224
0
                proj_destroy(proj_crs);
225
0
                return FALSE;
226
0
            }
227
228
8.33k
            {
229
8.33k
                const char* pszUnitCode = NULL;
230
8.33k
                if( !proj_cs_get_axis_info(
231
8.33k
                    ctx, coordSys, 0,
232
8.33k
                    NULL, /* name */
233
8.33k
                    NULL, /* abbreviation*/
234
8.33k
                    NULL, /* direction */
235
8.33k
                    NULL, /* conversion factor */
236
8.33k
                    NULL, /* unit name */
237
8.33k
                    NULL, /* unit auth name (should be EPSG) */
238
8.33k
                    &pszUnitCode) || pszUnitCode == NULL )
239
0
                {
240
0
                    proj_destroy(coordSys);
241
0
                    return FALSE;
242
0
                }
243
8.33k
                *pnUOMLengthCode = (short) atoi(pszUnitCode);
244
8.33k
                proj_destroy(coordSys);
245
8.33k
            }
246
8.33k
        }
247
248
16.4k
        if( pnGeogCS )
249
8.33k
        {
250
8.33k
            PJ* geod_crs = proj_crs_get_geodetic_crs(ctx, proj_crs);
251
8.33k
            if( !geod_crs )
252
0
            {
253
                // shouldn't happen except out of memory
254
0
                proj_destroy(proj_crs);
255
0
                return FALSE;
256
0
            }
257
258
8.33k
            {
259
8.33k
                const char* pszGeodCode = proj_get_id_code(geod_crs, 0);
260
8.33k
                assert( pszGeodCode );
261
8.33k
                *pnGeogCS = (short) atoi(pszGeodCode);
262
8.33k
            }
263
264
8.33k
            proj_destroy(geod_crs);
265
8.33k
        }
266
267
268
16.4k
        proj_destroy(proj_crs);
269
16.4k
        return TRUE;
270
16.4k
    }
271
16.4k
}
272
273
274
int GTIFGetPCSInfo( int nPCSCode, char **ppszEPSGName,
275
                    short *pnProjOp, short *pnUOMLengthCode,
276
                    short *pnGeogCS )
277
278
0
{
279
0
    PJ_CONTEXT* ctx = proj_context_create();
280
0
    int ret = GTIFGetPCSInfoEx(ctx, nPCSCode, ppszEPSGName, pnProjOp,
281
0
                               pnUOMLengthCode, pnGeogCS);
282
0
    proj_context_destroy(ctx);
283
0
    return ret;
284
0
}
285
286
/************************************************************************/
287
/*                           GTIFAngleToDD()                            */
288
/*                                                                      */
289
/*      Convert a numeric angle to decimal degrees.                     */
290
/************************************************************************/
291
292
double GTIFAngleToDD( double dfAngle, int nUOMAngle )
293
294
35.8k
{
295
35.8k
    if( nUOMAngle == 9110 )    /* DDD.MMSSsss */
296
56
    {
297
56
        if( dfAngle > -999.9 && dfAngle < 999.9 )
298
51
        {
299
51
            char  szAngleString[32];
300
301
51
            snprintf(szAngleString, sizeof(szAngleString), "%12.7f", dfAngle);
302
51
            dfAngle = GTIFAngleStringToDD( szAngleString, nUOMAngle );
303
51
        }
304
56
    }
305
35.7k
    else if ( nUOMAngle != KvUserDefined )
306
26.9k
    {
307
26.9k
        double    dfInDegrees = 1.0;
308
309
26.9k
        GTIFGetUOMAngleInfo( nUOMAngle, NULL, &dfInDegrees );
310
26.9k
        dfAngle = dfAngle * dfInDegrees;
311
26.9k
    }
312
313
35.8k
    return dfAngle;
314
35.8k
}
315
316
/************************************************************************/
317
/*                        GTIFAngleStringToDD()                         */
318
/*                                                                      */
319
/*      Convert an angle in the specified units to decimal degrees.     */
320
/************************************************************************/
321
322
double GTIFAngleStringToDD( const char * pszAngle, int nUOMAngle )
323
324
51
{
325
51
    double  dfAngle;
326
327
51
    if( nUOMAngle == 9110 )    /* DDD.MMSSsss */
328
51
    {
329
51
        char  *pszDecimal;
330
331
51
        dfAngle = ABS(atoi(pszAngle));
332
51
        pszDecimal = strchr(pszAngle,'.');
333
51
        if( pszDecimal != NULL && strlen(pszDecimal) > 1 )
334
51
        {
335
51
            char  szMinutes[3];
336
51
            char  szSeconds[64];
337
338
51
            szMinutes[0] = pszDecimal[1];
339
51
            if( pszDecimal[2] >= '0' && pszDecimal[2] <= '9' )
340
51
                szMinutes[1] = pszDecimal[2];
341
0
            else
342
0
                szMinutes[1] = '0';
343
344
51
            szMinutes[2] = '\0';
345
51
            dfAngle += atoi(szMinutes) / 60.0;
346
347
51
            if( strlen(pszDecimal) > 3 )
348
51
            {
349
51
                szSeconds[0] = pszDecimal[3];
350
51
                if( pszDecimal[4] >= '0' && pszDecimal[4] <= '9' )
351
51
                {
352
51
                    szSeconds[1] = pszDecimal[4];
353
51
                    szSeconds[2] = '.';
354
51
                    strncpy( szSeconds+3, pszDecimal + 5, sizeof(szSeconds) - 3 );
355
51
                    szSeconds[sizeof(szSeconds) - 1] = 0;
356
51
                }
357
0
                else
358
0
                {
359
0
                    szSeconds[1] = '0';
360
0
                    szSeconds[2] = '\0';
361
0
                }
362
51
                dfAngle += GTIFAtof(szSeconds) / 3600.0;
363
51
            }
364
51
        }
365
366
51
        if( pszAngle[0] == '-' )
367
0
            dfAngle *= -1;
368
51
    }
369
0
    else if( nUOMAngle == 9105 || nUOMAngle == 9106 )  /* grad */
370
0
    {
371
0
        dfAngle = 180 * (GTIFAtof(pszAngle ) / 200);
372
0
    }
373
0
    else if( nUOMAngle == 9101 )      /* radians */
374
0
    {
375
0
        dfAngle = 180 * (GTIFAtof(pszAngle ) / M_PI);
376
0
    }
377
0
    else if( nUOMAngle == 9103 )      /* arc-minute */
378
0
    {
379
0
        dfAngle = GTIFAtof(pszAngle) / 60;
380
0
    }
381
0
    else if( nUOMAngle == 9104 )      /* arc-second */
382
0
    {
383
0
        dfAngle = GTIFAtof(pszAngle) / 3600;
384
0
    }
385
0
    else /* decimal degrees ... some cases missing but seemingly never used */
386
0
    {
387
0
        CPLAssert( nUOMAngle == 9102 || nUOMAngle == KvUserDefined
388
0
                   || nUOMAngle == 0 );
389
390
0
        dfAngle = GTIFAtof(pszAngle );
391
0
    }
392
393
51
    return dfAngle;
394
51
}
395
396
/************************************************************************/
397
/*                           GTIFGetGCSInfo()                           */
398
/*                                                                      */
399
/*      Fetch the datum, and prime meridian related to a particular     */
400
/*      GCS.                                                            */
401
/************************************************************************/
402
403
int GTIFGetGCSInfoEx( void* ctxIn,
404
                      int nGCSCode, char ** ppszName,
405
                      short * pnDatum, short * pnPM, short *pnUOMAngle )
406
407
77.3k
{
408
77.3k
    int   nDatum=0;
409
410
/* -------------------------------------------------------------------- */
411
/*      Handle some "well known" GCS codes directly                     */
412
/* -------------------------------------------------------------------- */
413
77.3k
    const char * pszName = NULL;
414
77.3k
    const int nPM = PM_Greenwich;
415
77.3k
    const int nUOMAngle = Angular_DMS_Hemisphere;
416
77.3k
    if( nGCSCode == GCS_NAD27 )
417
4.82k
    {
418
4.82k
        nDatum = Datum_North_American_Datum_1927;
419
4.82k
        pszName = "NAD27";
420
4.82k
    }
421
72.5k
    else if( nGCSCode == GCS_NAD83 )
422
8.53k
    {
423
8.53k
        nDatum = Datum_North_American_Datum_1983;
424
8.53k
        pszName = "NAD83";
425
8.53k
    }
426
63.9k
    else if( nGCSCode == GCS_WGS_84 )
427
1.61k
    {
428
1.61k
        nDatum = Datum_WGS84;
429
1.61k
        pszName = "WGS 84";
430
1.61k
    }
431
62.3k
    else if( nGCSCode == GCS_WGS_72 )
432
899
    {
433
899
        nDatum = Datum_WGS72;
434
899
        pszName = "WGS 72";
435
899
    }
436
61.4k
    else if ( nGCSCode == KvUserDefined )
437
31.3k
    {
438
31.3k
        return FALSE;
439
31.3k
    }
440
441
46.0k
    if (pszName != NULL)
442
15.8k
    {
443
15.8k
        if( ppszName != NULL )
444
7.80k
            *ppszName = CPLStrdup( pszName );
445
15.8k
        if( pnDatum != NULL )
446
8.08k
            *pnDatum = (short) nDatum;
447
15.8k
        if( pnPM != NULL )
448
8.08k
            *pnPM = (short) nPM;
449
15.8k
        if( pnUOMAngle != NULL )
450
8.08k
            *pnUOMAngle = (short) nUOMAngle;
451
452
15.8k
        return TRUE;
453
15.8k
    }
454
455
/* -------------------------------------------------------------------- */
456
/*      Search the database.                                            */
457
/* -------------------------------------------------------------------- */
458
459
30.1k
    PJ_CONTEXT* ctx = (PJ_CONTEXT*)ctxIn;
460
30.1k
    {
461
30.1k
        char szCode[12];
462
463
30.1k
        snprintf(szCode, sizeof(szCode), "%d", nGCSCode);
464
30.1k
        PJ* geod_crs = proj_create_from_database(
465
30.1k
            ctx, "EPSG", szCode, PJ_CATEGORY_CRS, 0, NULL);
466
30.1k
        if( !geod_crs )
467
3.29k
        {
468
3.29k
            return FALSE;
469
3.29k
        }
470
471
26.8k
        {
472
26.8k
            const int objType = proj_get_type(geod_crs);
473
26.8k
            if( objType != PJ_TYPE_GEODETIC_CRS &&
474
26.8k
                objType != PJ_TYPE_GEOCENTRIC_CRS &&
475
25.9k
                objType != PJ_TYPE_GEOGRAPHIC_2D_CRS &&
476
4.51k
                objType != PJ_TYPE_GEOGRAPHIC_3D_CRS )
477
3.64k
            {
478
3.64k
                proj_destroy(geod_crs);
479
3.64k
                return FALSE;
480
3.64k
            }
481
26.8k
        }
482
483
23.1k
        if( ppszName )
484
11.5k
        {
485
11.5k
            pszName = proj_get_name(geod_crs);
486
11.5k
            if( !pszName )
487
0
            {
488
                // shouldn't happen
489
0
                proj_destroy(geod_crs);
490
0
                return FALSE;
491
0
            }
492
11.5k
            *ppszName = CPLStrdup(pszName);
493
11.5k
        }
494
495
23.1k
        if( pnDatum )
496
11.6k
        {
497
11.6k
#if PROJ_VERSION_MAJOR >= 8
498
11.6k
            PJ* datum = proj_crs_get_datum_forced(ctx, geod_crs);
499
#else
500
            PJ* datum = proj_crs_get_datum(ctx, geod_crs);
501
#endif
502
11.6k
            if( !datum )
503
0
            {
504
0
                proj_destroy(geod_crs);
505
0
                return FALSE;
506
0
            }
507
508
11.6k
            {
509
11.6k
                const char* pszDatumCode = proj_get_id_code(datum, 0);
510
11.6k
                assert( pszDatumCode );
511
11.6k
                *pnDatum = (short) atoi(pszDatumCode);
512
11.6k
            }
513
514
11.6k
            proj_destroy(datum);
515
11.6k
        }
516
517
23.1k
        if( pnPM )
518
11.6k
        {
519
11.6k
            PJ* pm = proj_get_prime_meridian(ctx, geod_crs);
520
11.6k
            if( !pm )
521
0
            {
522
0
                proj_destroy(geod_crs);
523
0
                return FALSE;
524
0
            }
525
526
11.6k
            {
527
11.6k
                const char* pszPMCode = proj_get_id_code(pm, 0);
528
11.6k
                assert( pszPMCode );
529
11.6k
                *pnPM = (short) atoi(pszPMCode);
530
11.6k
            }
531
532
11.6k
            proj_destroy(pm);
533
11.6k
        }
534
535
23.1k
        if( pnUOMAngle )
536
11.6k
        {
537
11.6k
            PJ* coordSys = proj_crs_get_coordinate_system(
538
11.6k
                ctx, geod_crs);
539
11.6k
            if( !coordSys )
540
0
            {
541
                // shouldn't happen except out of memory
542
0
                proj_destroy(geod_crs);
543
0
                return FALSE;
544
0
            }
545
546
11.6k
            {
547
11.6k
                const char* pszUnitCode = NULL;
548
11.6k
                if( !proj_cs_get_axis_info(
549
11.6k
                    ctx, coordSys, 0,
550
11.6k
                    NULL, /* name */
551
11.6k
                    NULL, /* abbreviation*/
552
11.6k
                    NULL, /* direction */
553
11.6k
                    NULL, /* conversion factor */
554
11.6k
                    NULL, /* unit name */
555
11.6k
                    NULL, /* unit auth name (should be EPSG) */
556
11.6k
                    &pszUnitCode) || pszUnitCode == NULL )
557
0
                {
558
0
                    proj_destroy(coordSys);
559
0
                    return FALSE;
560
0
                }
561
11.6k
                *pnUOMAngle = (short) atoi(pszUnitCode);
562
11.6k
                proj_destroy(coordSys);
563
11.6k
            }
564
11.6k
        }
565
566
23.1k
        proj_destroy(geod_crs);
567
23.1k
        return TRUE;
568
23.1k
    }
569
23.1k
}
570
571
int GTIFGetGCSInfo( int nGCSCode, char ** ppszName,
572
                    short * pnDatum, short * pnPM, short *pnUOMAngle )
573
574
0
{
575
0
    PJ_CONTEXT* ctx = proj_context_create();
576
0
    const int ret = GTIFGetGCSInfoEx(ctx, nGCSCode, ppszName, pnDatum,
577
0
                                     pnPM, pnUOMAngle);
578
0
    proj_context_destroy(ctx);
579
0
    return ret;
580
0
}
581
582
/************************************************************************/
583
/*                        GTIFGetEllipsoidInfo()                        */
584
/*                                                                      */
585
/*      Fetch info about an ellipsoid.  Axes are always returned in     */
586
/*      meters.  SemiMajor computed based on inverse flattening         */
587
/*      where that is provided.                                         */
588
/************************************************************************/
589
590
int GTIFGetEllipsoidInfoEx( void* ctxIn,
591
                            int nEllipseCode, char ** ppszName,
592
                            double * pdfSemiMajor, double * pdfSemiMinor )
593
594
73.4k
{
595
73.4k
    PJ_CONTEXT* ctx = (PJ_CONTEXT*)ctxIn;
596
/* -------------------------------------------------------------------- */
597
/*      Try some well known ellipsoids.                                 */
598
/* -------------------------------------------------------------------- */
599
73.4k
    double  dfSemiMajor=0.0;
600
73.4k
    double     dfInvFlattening=0.0, dfSemiMinor=0.0;
601
73.4k
    const char *pszName = NULL;
602
603
73.4k
    if( nEllipseCode == Ellipse_Clarke_1866 )
604
5.11k
    {
605
5.11k
        pszName = "Clarke 1866";
606
5.11k
        dfSemiMajor = 6378206.4;
607
5.11k
        dfSemiMinor = 6356583.8;
608
5.11k
        dfInvFlattening = 0.0;
609
5.11k
    }
610
68.2k
    else if( nEllipseCode == Ellipse_GRS_1980 )
611
13.4k
    {
612
13.4k
        pszName = "GRS 1980";
613
13.4k
        dfSemiMajor = 6378137.0;
614
13.4k
        dfSemiMinor = 0.0;
615
13.4k
        dfInvFlattening = 298.257222101;
616
13.4k
    }
617
54.8k
    else if( nEllipseCode == Ellipse_WGS_84 )
618
2.04k
    {
619
2.04k
        pszName = "WGS 84";
620
2.04k
        dfSemiMajor = 6378137.0;
621
2.04k
        dfSemiMinor = 0.0;
622
2.04k
        dfInvFlattening = 298.257223563;
623
2.04k
    }
624
52.7k
    else if( nEllipseCode == 7043 )
625
8.88k
    {
626
8.88k
        pszName = "WGS 72";
627
8.88k
        dfSemiMajor = 6378135.0;
628
8.88k
        dfSemiMinor = 0.0;
629
8.88k
        dfInvFlattening = 298.26;
630
8.88k
    }
631
632
73.4k
    if (pszName != NULL)
633
29.5k
    {
634
29.5k
        if( dfSemiMinor == 0.0 )
635
24.4k
            dfSemiMinor = dfSemiMajor * (1 - 1.0/dfInvFlattening);
636
637
29.5k
        if( pdfSemiMinor != NULL )
638
14.9k
            *pdfSemiMinor = dfSemiMinor;
639
29.5k
        if( pdfSemiMajor != NULL )
640
14.9k
            *pdfSemiMajor = dfSemiMajor;
641
29.5k
        if( ppszName != NULL )
642
14.5k
            *ppszName = CPLStrdup( pszName );
643
644
29.5k
        return TRUE;
645
29.5k
    }
646
647
43.8k
    if( nEllipseCode == KvUserDefined )
648
34.8k
        return FALSE;
649
650
/* -------------------------------------------------------------------- */
651
/*      Search the database.                                            */
652
/* -------------------------------------------------------------------- */
653
9.04k
    {
654
9.04k
        char szCode[12];
655
656
9.04k
        snprintf(szCode, sizeof(szCode), "%d", nEllipseCode);
657
9.04k
        PJ* ellipsoid = proj_create_from_database(
658
9.04k
            ctx, "EPSG", szCode, PJ_CATEGORY_ELLIPSOID, 0, NULL);
659
9.04k
        if( !ellipsoid )
660
3.09k
        {
661
3.09k
            return FALSE;
662
3.09k
        }
663
664
5.94k
        if( ppszName )
665
2.95k
        {
666
2.95k
            pszName = proj_get_name(ellipsoid);
667
2.95k
            if( !pszName )
668
0
            {
669
                // shouldn't happen
670
0
                proj_destroy(ellipsoid);
671
0
                return FALSE;
672
0
            }
673
2.95k
            *ppszName = CPLStrdup(pszName);
674
2.95k
        }
675
676
5.94k
        proj_ellipsoid_get_parameters(
677
5.94k
            ctx, ellipsoid, pdfSemiMajor, pdfSemiMinor, NULL, NULL);
678
679
5.94k
        proj_destroy(ellipsoid);
680
681
5.94k
        return TRUE;
682
5.94k
    }
683
5.94k
}
684
685
int GTIFGetEllipsoidInfo( int nEllipseCode, char ** ppszName,
686
                          double * pdfSemiMajor, double * pdfSemiMinor )
687
688
0
{
689
0
    PJ_CONTEXT* ctx = proj_context_create();
690
0
    const int ret = GTIFGetEllipsoidInfoEx(ctx, nEllipseCode, ppszName, pdfSemiMajor,
691
0
                                           pdfSemiMinor);
692
0
    proj_context_destroy(ctx);
693
0
    return ret;
694
0
}
695
696
/************************************************************************/
697
/*                           GTIFGetPMInfo()                            */
698
/*                                                                      */
699
/*      Get the offset between a given prime meridian and Greenwich     */
700
/*      in degrees.                                                     */
701
/************************************************************************/
702
703
int GTIFGetPMInfoEx( void* ctxIn,
704
                     int nPMCode, char ** ppszName, double *pdfOffset )
705
706
73.8k
{
707
/* -------------------------------------------------------------------- */
708
/*      Use a special short cut for Greenwich, since it is so common.   */
709
/* -------------------------------------------------------------------- */
710
73.8k
    if( nPMCode == PM_Greenwich )
711
38.6k
    {
712
38.6k
        if( pdfOffset != NULL )
713
19.5k
            *pdfOffset = 0.0;
714
38.6k
        if( ppszName != NULL )
715
19.1k
            *ppszName = CPLStrdup( "Greenwich" );
716
38.6k
        return TRUE;
717
38.6k
    }
718
719
720
35.1k
    if( nPMCode == KvUserDefined )
721
34.2k
        return FALSE;
722
723
/* -------------------------------------------------------------------- */
724
/*      Search the database.                                            */
725
/* -------------------------------------------------------------------- */
726
948
    PJ_CONTEXT* ctx = (PJ_CONTEXT*)ctxIn;
727
948
    {
728
948
        char szCode[12];
729
730
948
        snprintf(szCode, sizeof(szCode), "%d", nPMCode);
731
948
        PJ* pm = proj_create_from_database(
732
948
            ctx, "EPSG", szCode, PJ_CATEGORY_PRIME_MERIDIAN, 0, NULL);
733
948
        if( !pm )
734
633
        {
735
633
            return FALSE;
736
633
        }
737
738
315
        if( ppszName )
739
155
        {
740
155
            const char* pszName = proj_get_name(pm);
741
155
            if( !pszName )
742
0
            {
743
                // shouldn't happen
744
0
                proj_destroy(pm);
745
0
                return FALSE;
746
0
            }
747
155
            *ppszName = CPLStrdup(pszName);
748
155
        }
749
750
315
        if( pdfOffset )
751
160
        {
752
160
            double conv_factor = 0;
753
160
            proj_prime_meridian_get_parameters(
754
160
                ctx, pm, pdfOffset, &conv_factor, NULL);
755
160
            *pdfOffset *= conv_factor * 180.0 / M_PI;
756
160
        }
757
758
315
        proj_destroy(pm);
759
760
315
        return TRUE;
761
315
    }
762
315
}
763
764
int GTIFGetPMInfo( int nPMCode, char ** ppszName, double *pdfOffset )
765
766
0
{
767
0
    PJ_CONTEXT* ctx = proj_context_create();
768
0
    const int ret = GTIFGetPMInfoEx(ctx, nPMCode, ppszName, pdfOffset);
769
0
    proj_context_destroy(ctx);
770
0
    return ret;
771
0
}
772
773
/************************************************************************/
774
/*                          GTIFGetDatumInfo()                          */
775
/*                                                                      */
776
/*      Fetch the ellipsoid, and name for a datum.                      */
777
/************************************************************************/
778
779
int GTIFGetDatumInfoEx( void* ctxIn,
780
                        int nDatumCode, char ** ppszName, short * pnEllipsoid )
781
782
73.2k
{
783
73.2k
    const char* pszName = NULL;
784
73.2k
    int   nEllipsoid = 0;
785
786
/* -------------------------------------------------------------------- */
787
/*      Handle a few built-in datums.                                   */
788
/* -------------------------------------------------------------------- */
789
73.2k
    if( nDatumCode == Datum_North_American_Datum_1927 )
790
4.76k
    {
791
4.76k
        nEllipsoid = Ellipse_Clarke_1866;
792
4.76k
        pszName = "North American Datum 1927";
793
4.76k
    }
794
68.4k
    else if( nDatumCode == Datum_North_American_Datum_1983 )
795
8.23k
    {
796
8.23k
        nEllipsoid = Ellipse_GRS_1980;
797
8.23k
        pszName = "North American Datum 1983";
798
8.23k
    }
799
60.2k
    else if( nDatumCode == Datum_WGS84 )
800
1.62k
    {
801
1.62k
        nEllipsoid = Ellipse_WGS_84;
802
1.62k
        pszName = "World Geodetic System 1984";
803
1.62k
    }
804
58.6k
    else if( nDatumCode == Datum_WGS72 )
805
811
    {
806
811
        nEllipsoid = 7043; /* WGS72 */
807
811
        pszName = "World Geodetic System 1972";
808
811
    }
809
810
73.2k
    if (pszName != NULL)
811
15.4k
    {
812
15.4k
        if( pnEllipsoid != NULL )
813
7.85k
            *pnEllipsoid = (short) nEllipsoid;
814
815
15.4k
        if( ppszName != NULL )
816
7.57k
            *ppszName = CPLStrdup( pszName );
817
818
15.4k
        return TRUE;
819
15.4k
    }
820
821
57.8k
    if( nDatumCode == KvUserDefined )
822
28.8k
        return FALSE;
823
824
/* -------------------------------------------------------------------- */
825
/*      Search the database.                                            */
826
/* -------------------------------------------------------------------- */
827
28.9k
    PJ_CONTEXT* ctx = (PJ_CONTEXT*)ctxIn;
828
28.9k
    {
829
28.9k
        char szCode[12];
830
831
28.9k
        snprintf(szCode, sizeof(szCode), "%d", nDatumCode);
832
28.9k
        PJ* datum = proj_create_from_database(
833
28.9k
            ctx, "EPSG", szCode, PJ_CATEGORY_DATUM, 0, NULL);
834
28.9k
        if( !datum )
835
8.73k
        {
836
8.73k
            return FALSE;
837
8.73k
        }
838
839
20.2k
        const PJ_TYPE pjType = proj_get_type(datum);
840
20.2k
        if( pjType != PJ_TYPE_GEODETIC_REFERENCE_FRAME &&
841
8.91k
            pjType != PJ_TYPE_DYNAMIC_GEODETIC_REFERENCE_FRAME )
842
90
        {
843
90
            proj_destroy(datum);
844
90
            return FALSE;
845
90
        }
846
847
20.1k
        if( ppszName )
848
10.0k
        {
849
10.0k
            pszName = proj_get_name(datum);
850
10.0k
            if( !pszName )
851
0
            {
852
                // shouldn't happen
853
0
                proj_destroy(datum);
854
0
                return FALSE;
855
0
            }
856
10.0k
            *ppszName = CPLStrdup(pszName);
857
10.0k
        }
858
859
20.1k
        if( pnEllipsoid )
860
10.1k
        {
861
10.1k
            PJ* ellipsoid = proj_get_ellipsoid(ctx, datum);
862
10.1k
            if( !ellipsoid )
863
0
            {
864
0
                proj_destroy(datum);
865
0
                return FALSE;
866
0
            }
867
868
10.1k
            {
869
10.1k
                const char* pszEllipsoidCode = proj_get_id_code(
870
10.1k
                    ellipsoid, 0);
871
10.1k
                assert( pszEllipsoidCode );
872
10.1k
                *pnEllipsoid = (short) atoi(pszEllipsoidCode);
873
10.1k
            }
874
875
10.1k
            proj_destroy(ellipsoid);
876
10.1k
        }
877
878
20.1k
        proj_destroy(datum);
879
880
20.1k
        return TRUE;
881
20.1k
    }
882
20.1k
}
883
884
int GTIFGetDatumInfo( int nDatumCode, char ** ppszName, short * pnEllipsoid )
885
886
0
{
887
0
    PJ_CONTEXT* ctx = proj_context_create();
888
0
    const int ret = GTIFGetDatumInfoEx(ctx, nDatumCode, ppszName, pnEllipsoid);
889
0
    proj_context_destroy(ctx);
890
0
    return ret;
891
0
}
892
893
/************************************************************************/
894
/*                        GTIFGetUOMLengthInfo()                        */
895
/*                                                                      */
896
/*      Note: This function should eventually also know how to          */
897
/*      lookup length aliases in the UOM_LE_ALIAS table.                */
898
/************************************************************************/
899
900
int GTIFGetUOMLengthInfoEx( void* ctxIn,
901
                            int nUOMLengthCode,
902
                            char **ppszUOMName,
903
                            double * pdfInMeters )
904
905
77.2k
{
906
/* -------------------------------------------------------------------- */
907
/*      We short cut meter to save work and avoid failure for missing   */
908
/*      in the most common cases.               */
909
/* -------------------------------------------------------------------- */
910
77.2k
    if( nUOMLengthCode == 9001 )
911
36.6k
    {
912
36.6k
        if( ppszUOMName != NULL )
913
18.8k
            *ppszUOMName = CPLStrdup( "metre" );
914
36.6k
        if( pdfInMeters != NULL )
915
17.8k
            *pdfInMeters = 1.0;
916
917
36.6k
        return TRUE;
918
36.6k
    }
919
920
40.5k
    if( nUOMLengthCode == 9002 )
921
16.7k
    {
922
16.7k
        if( ppszUOMName != NULL )
923
8.46k
            *ppszUOMName = CPLStrdup( "foot" );
924
16.7k
        if( pdfInMeters != NULL )
925
8.24k
            *pdfInMeters = 0.3048;
926
927
16.7k
        return TRUE;
928
16.7k
    }
929
930
23.8k
    if( nUOMLengthCode == 9003 )
931
2.60k
    {
932
2.60k
        if( ppszUOMName != NULL )
933
1.29k
            *ppszUOMName = CPLStrdup( "US survey foot" );
934
2.60k
        if( pdfInMeters != NULL )
935
1.30k
            *pdfInMeters = 12.0 / 39.37;
936
937
2.60k
        return TRUE;
938
2.60k
    }
939
940
21.2k
    if( nUOMLengthCode == KvUserDefined )
941
5.66k
        return FALSE;
942
943
/* -------------------------------------------------------------------- */
944
/*      Search the units database for this unit.  If we don't find      */
945
/*      it return failure.                                              */
946
/* -------------------------------------------------------------------- */
947
15.6k
    char szCode[12];
948
15.6k
    const char* pszName = NULL;
949
950
15.6k
    snprintf(szCode, sizeof(szCode), "%d", nUOMLengthCode);
951
15.6k
    PJ_CONTEXT* ctx = (PJ_CONTEXT*)ctxIn;
952
15.6k
    if( !proj_uom_get_info_from_database(
953
15.6k
        ctx, "EPSG", szCode, &pszName, pdfInMeters,  NULL) )
954
14.6k
    {
955
14.6k
        return FALSE;
956
14.6k
    }
957
931
    if( ppszUOMName )
958
524
    {
959
524
        *ppszUOMName = CPLStrdup(pszName);
960
524
    }
961
931
    return TRUE;
962
15.6k
}
963
964
int GTIFGetUOMLengthInfo( int nUOMLengthCode,
965
                          char **ppszUOMName,
966
                          double * pdfInMeters )
967
968
10.2k
{
969
10.2k
    PJ_CONTEXT* ctx = proj_context_create();
970
10.2k
    const int ret = GTIFGetUOMLengthInfoEx(
971
10.2k
        ctx, nUOMLengthCode, ppszUOMName, pdfInMeters);
972
10.2k
    proj_context_destroy(ctx);
973
10.2k
    return ret;
974
10.2k
}
975
976
/************************************************************************/
977
/*                        GTIFGetUOMAngleInfo()                         */
978
/************************************************************************/
979
980
int GTIFGetUOMAngleInfoEx( void* ctxIn,
981
                           int nUOMAngleCode,
982
                           char **ppszUOMName,
983
                           double * pdfInDegrees )
984
985
127k
{
986
127k
    const char  *pszUOMName = NULL;
987
127k
    double  dfInDegrees = 1.0;
988
989
127k
    switch( nUOMAngleCode )
990
127k
    {
991
4.66k
      case 9101:
992
4.66k
        pszUOMName = "radian";
993
4.66k
        dfInDegrees = 180.0 / M_PI;
994
4.66k
        break;
995
996
2.47k
      case 9102:
997
2.51k
      case 9107:
998
7.17k
      case 9108:
999
7.28k
      case 9110:
1000
14.5k
      case 9122:
1001
14.5k
        pszUOMName = "degree";
1002
14.5k
        dfInDegrees = 1.0;
1003
14.5k
        break;
1004
1005
38.7k
      case 9103:
1006
38.7k
        pszUOMName = "arc-minute";
1007
38.7k
        dfInDegrees = 1 / 60.0;
1008
38.7k
        break;
1009
1010
1.12k
      case 9104:
1011
1.12k
        pszUOMName = "arc-second";
1012
1.12k
        dfInDegrees = 1 / 3600.0;
1013
1.12k
        break;
1014
1015
2.89k
      case 9105:
1016
2.89k
        pszUOMName = "grad";
1017
2.89k
        dfInDegrees = 180.0 / 200.0;
1018
2.89k
        break;
1019
1020
14.9k
      case 9106:
1021
14.9k
        pszUOMName = "gon";
1022
14.9k
        dfInDegrees = 180.0 / 200.0;
1023
14.9k
        break;
1024
1025
66
      case 9109:
1026
66
        pszUOMName = "microradian";
1027
66
        dfInDegrees = 180.0 / (M_PI * 1000000.0);
1028
66
        break;
1029
1030
50.6k
      default:
1031
50.6k
        break;
1032
127k
    }
1033
1034
127k
    if (pszUOMName)
1035
77.0k
    {
1036
77.0k
        if( ppszUOMName != NULL )
1037
30.3k
        {
1038
30.3k
            *ppszUOMName = CPLStrdup( pszUOMName );
1039
30.3k
        }
1040
1041
77.0k
        if( pdfInDegrees != NULL )
1042
77.0k
            *pdfInDegrees = dfInDegrees;
1043
1044
77.0k
        return TRUE;
1045
77.0k
    }
1046
1047
50.6k
    if( nUOMAngleCode == KvUserDefined )
1048
8.35k
        return FALSE;
1049
1050
/* -------------------------------------------------------------------- */
1051
/*      Search the units database for this unit.  If we don't find      */
1052
/*      it return failure.                                              */
1053
/* -------------------------------------------------------------------- */
1054
42.2k
    char szCode[12];
1055
42.2k
    const char* pszName = NULL;
1056
42.2k
    double dfConvFactorToRadians = 0;
1057
1058
42.2k
    snprintf(szCode, sizeof(szCode), "%d", nUOMAngleCode);
1059
42.2k
    PJ_CONTEXT* ctx = (PJ_CONTEXT*)ctxIn;
1060
42.2k
    if( !proj_uom_get_info_from_database(
1061
42.2k
        ctx, "EPSG", szCode, &pszName, &dfConvFactorToRadians, NULL) )
1062
36.1k
    {
1063
36.1k
        return FALSE;
1064
36.1k
    }
1065
6.12k
    if( ppszUOMName )
1066
2.31k
    {
1067
2.31k
        *ppszUOMName = CPLStrdup(pszName);
1068
2.31k
    }
1069
6.12k
    if( pdfInDegrees )
1070
6.12k
    {
1071
6.12k
        *pdfInDegrees = dfConvFactorToRadians * 180.0 / M_PI;
1072
6.12k
    }
1073
6.12k
    return TRUE;
1074
42.2k
}
1075
1076
int GTIFGetUOMAngleInfo( int nUOMAngleCode,
1077
                         char **ppszUOMName,
1078
                         double * pdfInDegrees )
1079
1080
26.9k
{
1081
26.9k
    PJ_CONTEXT* ctx = proj_context_create();
1082
26.9k
    const int ret = GTIFGetUOMAngleInfoEx(
1083
26.9k
        ctx, nUOMAngleCode, ppszUOMName, pdfInDegrees);
1084
26.9k
    proj_context_destroy(ctx);
1085
26.9k
    return ret;
1086
26.9k
}
1087
1088
/************************************************************************/
1089
/*                    EPSGProjMethodToCTProjMethod()                    */
1090
/*                                                                      */
1091
/*      Convert between the EPSG enumeration for projection methods,    */
1092
/*      and the GeoTIFF CT codes.                                       */
1093
/************************************************************************/
1094
1095
static int EPSGProjMethodToCTProjMethod( int nEPSG, int bReturnExtendedCTCode )
1096
1097
56.2k
{
1098
56.2k
    switch( nEPSG )
1099
56.2k
    {
1100
753
      case 9801:
1101
753
        return CT_LambertConfConic_1SP;
1102
1103
6.19k
      case 9802:
1104
6.19k
        return CT_LambertConfConic_2SP;
1105
1106
0
      case 9803:
1107
0
        return CT_LambertConfConic_2SP; /* Belgian variant not supported */
1108
1109
108
      case 9804:
1110
108
        return CT_Mercator;  /* 1SP and 2SP not differentiated */
1111
1112
0
      case 9805:
1113
0
        if( bReturnExtendedCTCode )
1114
0
            return CT_Ext_Mercator_2SP;
1115
0
        else
1116
0
            return CT_Mercator;  /* 1SP and 2SP not differentiated */
1117
1118
      /* Mercator 1SP (Spherical) For EPSG:3785 */
1119
51
      case 9841:
1120
51
        return CT_Mercator;  /* 1SP and 2SP not differentiated */
1121
1122
      /* Google Mercator For EPSG:3857 */
1123
141
      case 1024:
1124
141
        return CT_Mercator;  /* 1SP and 2SP not differentiated */
1125
1126
312
      case 9806:
1127
312
        return CT_CassiniSoldner;
1128
1129
32.6k
      case 9807:
1130
32.6k
        return CT_TransverseMercator;
1131
1132
360
      case 9808:
1133
360
        return CT_TransvMercator_SouthOriented;
1134
1135
54
      case 9809:
1136
54
        return CT_ObliqueStereographic;
1137
1138
744
      case 9810:
1139
906
      case 9829: /* variant B not quite the same - not sure how to handle */
1140
906
        return CT_PolarStereographic;
1141
1142
3
      case 9811:
1143
3
        return CT_NewZealandMapGrid;
1144
1145
99
      case 9812:
1146
99
        return CT_ObliqueMercator; /* is hotine actually different? */
1147
1148
12
      case 9813:
1149
12
        return CT_ObliqueMercator_Laborde;
1150
1151
0
      case 9814:
1152
0
        return CT_ObliqueMercator_Rosenmund; /* swiss  */
1153
1154
288
      case 9815:
1155
288
        return CT_HotineObliqueMercatorAzimuthCenter;
1156
1157
0
      case 9816: /* tunesia mining grid has no counterpart */
1158
0
        return KvUserDefined;
1159
1160
6
      case 9818:
1161
6
        return CT_Polyconic;
1162
1163
234
      case 9820:
1164
315
      case 1027:
1165
315
        return CT_LambertAzimEqualArea;
1166
1167
414
      case 9822:
1168
414
        return CT_AlbersEqualArea;
1169
1170
3
      case 9834:
1171
3
        return CT_CylindricalEqualArea;
1172
1173
6
      case 1028:
1174
54
      case 1029:
1175
144
      case 9823: /* spherical */
1176
225
      case 9842: /* elliptical */
1177
225
        return CT_Equirectangular;
1178
1179
13.3k
      default: /* use the EPSG code for other methods */
1180
13.3k
        return nEPSG;
1181
56.2k
    }
1182
56.2k
}
1183
1184
/************************************************************************/
1185
/*                           SetGTParamIds()                            */
1186
/*                                                                      */
1187
/*      This is hardcoded logic to set the GeoTIFF parameter            */
1188
/*      identifiers for all the EPSG supported projections.  As new     */
1189
/*      projection methods are added, this code will need to be updated */
1190
/************************************************************************/
1191
1192
static int SetGTParamIds( int nCTProjection,
1193
                          int nEPSGProjMethod,
1194
                          int *panProjParamId,
1195
                          int *panEPSGCodes )
1196
1197
32.4k
{
1198
32.4k
    int anWorkingDummy[7];
1199
1200
32.4k
    if( panEPSGCodes == NULL )
1201
23.7k
        panEPSGCodes = anWorkingDummy;
1202
32.4k
    if( panProjParamId == NULL )
1203
8.68k
        panProjParamId = anWorkingDummy;
1204
1205
32.4k
    memset( panEPSGCodes, 0, sizeof(int) * 7 );
1206
1207
    /* psDefn->nParms = 7; */
1208
1209
32.4k
    switch( nCTProjection )
1210
32.4k
    {
1211
208
      case CT_CassiniSoldner:
1212
210
      case CT_NewZealandMapGrid:
1213
214
      case CT_Polyconic:
1214
214
        panProjParamId[0] = ProjNatOriginLatGeoKey;
1215
214
        panProjParamId[1] = ProjNatOriginLongGeoKey;
1216
214
        panProjParamId[5] = ProjFalseEastingGeoKey;
1217
214
        panProjParamId[6] = ProjFalseNorthingGeoKey;
1218
1219
214
        panEPSGCodes[0] = EPSGNatOriginLat;
1220
214
        panEPSGCodes[1] = EPSGNatOriginLong;
1221
214
        panEPSGCodes[5] = EPSGFalseEasting;
1222
214
        panEPSGCodes[6] = EPSGFalseNorthing;
1223
214
        return TRUE;
1224
1225
66
      case CT_ObliqueMercator:
1226
258
      case CT_HotineObliqueMercatorAzimuthCenter:
1227
258
        panProjParamId[0] = ProjCenterLatGeoKey;
1228
258
        panProjParamId[1] = ProjCenterLongGeoKey;
1229
258
        panProjParamId[2] = ProjAzimuthAngleGeoKey;
1230
258
        panProjParamId[3] = ProjRectifiedGridAngleGeoKey;
1231
258
        panProjParamId[4] = ProjScaleAtCenterGeoKey;
1232
258
        panProjParamId[5] = ProjFalseEastingGeoKey;
1233
258
        panProjParamId[6] = ProjFalseNorthingGeoKey;
1234
1235
258
        panEPSGCodes[0] = EPSGProjCenterLat;
1236
258
        panEPSGCodes[1] = EPSGProjCenterLong;
1237
258
        panEPSGCodes[2] = EPSGAzimuth;
1238
258
        panEPSGCodes[3] = EPSGAngleRectifiedToSkewedGrid;
1239
258
        panEPSGCodes[4] = EPSGInitialLineScaleFactor;
1240
258
        panEPSGCodes[5] = EPSGProjCenterEasting; /* EPSG proj method 9812 uses EPSGFalseEasting, but 9815 uses EPSGProjCenterEasting */
1241
258
        panEPSGCodes[6] = EPSGProjCenterNorthing; /* EPSG proj method 9812 uses EPSGFalseNorthing, but 9815 uses EPSGProjCenterNorthing */
1242
258
        return TRUE;
1243
1244
8
      case CT_ObliqueMercator_Laborde:
1245
8
        panProjParamId[0] = ProjCenterLatGeoKey;
1246
8
        panProjParamId[1] = ProjCenterLongGeoKey;
1247
8
        panProjParamId[2] = ProjAzimuthAngleGeoKey;
1248
8
        panProjParamId[4] = ProjScaleAtCenterGeoKey;
1249
8
        panProjParamId[5] = ProjFalseEastingGeoKey;
1250
8
        panProjParamId[6] = ProjFalseNorthingGeoKey;
1251
1252
8
        panEPSGCodes[0] = EPSGProjCenterLat;
1253
8
        panEPSGCodes[1] = EPSGProjCenterLong;
1254
8
        panEPSGCodes[2] = EPSGAzimuth;
1255
8
        panEPSGCodes[4] = EPSGInitialLineScaleFactor;
1256
8
        panEPSGCodes[5] = EPSGFalseEasting;
1257
8
        panEPSGCodes[6] = EPSGFalseNorthing;
1258
8
        return TRUE;
1259
1260
502
      case CT_LambertConfConic_1SP:
1261
702
      case CT_Mercator:
1262
738
      case CT_ObliqueStereographic:
1263
1.34k
      case CT_PolarStereographic:
1264
20.1k
      case CT_TransverseMercator:
1265
20.3k
      case CT_TransvMercator_SouthOriented:
1266
20.3k
        panProjParamId[0] = ProjNatOriginLatGeoKey;
1267
20.3k
        if( nCTProjection == CT_PolarStereographic )
1268
604
        {
1269
604
            panProjParamId[1] = ProjStraightVertPoleLongGeoKey;
1270
604
        }
1271
19.7k
        else
1272
19.7k
        {
1273
19.7k
            panProjParamId[1] = ProjNatOriginLongGeoKey;
1274
19.7k
        }
1275
20.3k
        if( nEPSGProjMethod == 9805 ) /* Mercator_2SP */
1276
0
        {
1277
0
            panProjParamId[2] = ProjStdParallel1GeoKey;
1278
0
        }
1279
20.3k
        panProjParamId[4] = ProjScaleAtNatOriginGeoKey;
1280
20.3k
        panProjParamId[5] = ProjFalseEastingGeoKey;
1281
20.3k
        panProjParamId[6] = ProjFalseNorthingGeoKey;
1282
1283
20.3k
        panEPSGCodes[0] = EPSGNatOriginLat;
1284
20.3k
        panEPSGCodes[1] = EPSGNatOriginLong;
1285
20.3k
        if( nEPSGProjMethod == 9805 ) /* Mercator_2SP */
1286
0
        {
1287
0
            panEPSGCodes[2] = EPSGStdParallel1Lat;
1288
0
        }
1289
20.3k
        panEPSGCodes[4] = EPSGNatOriginScaleFactor;
1290
20.3k
        panEPSGCodes[5] = EPSGFalseEasting;
1291
20.3k
        panEPSGCodes[6] = EPSGFalseNorthing;
1292
20.3k
        return TRUE;
1293
1294
4.13k
      case CT_LambertConfConic_2SP:
1295
4.13k
        panProjParamId[0] = ProjFalseOriginLatGeoKey;
1296
4.13k
        panProjParamId[1] = ProjFalseOriginLongGeoKey;
1297
4.13k
        panProjParamId[2] = ProjStdParallel1GeoKey;
1298
4.13k
        panProjParamId[3] = ProjStdParallel2GeoKey;
1299
4.13k
        panProjParamId[5] = ProjFalseEastingGeoKey;
1300
4.13k
        panProjParamId[6] = ProjFalseNorthingGeoKey;
1301
1302
4.13k
        panEPSGCodes[0] = EPSGFalseOriginLat;
1303
4.13k
        panEPSGCodes[1] = EPSGFalseOriginLong;
1304
4.13k
        panEPSGCodes[2] = EPSGStdParallel1Lat;
1305
4.13k
        panEPSGCodes[3] = EPSGStdParallel2Lat;
1306
4.13k
        panEPSGCodes[5] = EPSGFalseOriginEasting;
1307
4.13k
        panEPSGCodes[6] = EPSGFalseOriginNorthing;
1308
4.13k
        return TRUE;
1309
1310
276
      case CT_AlbersEqualArea:
1311
276
        panProjParamId[0] = ProjStdParallel1GeoKey;
1312
276
        panProjParamId[1] = ProjStdParallel2GeoKey;
1313
276
        panProjParamId[2] = ProjNatOriginLatGeoKey;
1314
276
        panProjParamId[3] = ProjNatOriginLongGeoKey;
1315
276
        panProjParamId[5] = ProjFalseEastingGeoKey;
1316
276
        panProjParamId[6] = ProjFalseNorthingGeoKey;
1317
1318
276
        panEPSGCodes[0] = EPSGStdParallel1Lat;
1319
276
        panEPSGCodes[1] = EPSGStdParallel2Lat;
1320
276
        panEPSGCodes[2] = EPSGFalseOriginLat;
1321
276
        panEPSGCodes[3] = EPSGFalseOriginLong;
1322
276
        panEPSGCodes[5] = EPSGFalseOriginEasting;
1323
276
        panEPSGCodes[6] = EPSGFalseOriginNorthing;
1324
276
        return TRUE;
1325
1326
0
      case CT_SwissObliqueCylindrical:
1327
0
        panProjParamId[0] = ProjCenterLatGeoKey;
1328
0
        panProjParamId[1] = ProjCenterLongGeoKey;
1329
0
        panProjParamId[5] = ProjFalseEastingGeoKey;
1330
0
        panProjParamId[6] = ProjFalseNorthingGeoKey;
1331
1332
        /* EPSG codes? */
1333
0
        return TRUE;
1334
1335
210
      case CT_LambertAzimEqualArea:
1336
210
        panProjParamId[0] = ProjCenterLatGeoKey;
1337
210
        panProjParamId[1] = ProjCenterLongGeoKey;
1338
210
        panProjParamId[5] = ProjFalseEastingGeoKey;
1339
210
        panProjParamId[6] = ProjFalseNorthingGeoKey;
1340
1341
210
        panEPSGCodes[0] = EPSGNatOriginLat;
1342
210
        panEPSGCodes[1] = EPSGNatOriginLong;
1343
210
        panEPSGCodes[5] = EPSGFalseEasting;
1344
210
        panEPSGCodes[6] = EPSGFalseNorthing;
1345
210
        return TRUE;
1346
1347
2
      case CT_CylindricalEqualArea:
1348
2
        panProjParamId[0] = ProjStdParallel1GeoKey;
1349
2
        panProjParamId[1] = ProjNatOriginLongGeoKey;
1350
2
        panProjParamId[5] = ProjFalseEastingGeoKey;
1351
2
        panProjParamId[6] = ProjFalseNorthingGeoKey;
1352
1353
2
        panEPSGCodes[0] = EPSGStdParallel1Lat;
1354
2
        panEPSGCodes[1] = EPSGFalseOriginLong;
1355
2
        panEPSGCodes[5] = EPSGFalseOriginEasting;
1356
2
        panEPSGCodes[6] = EPSGFalseOriginNorthing;
1357
2
        return TRUE;
1358
1359
150
      case CT_Equirectangular:
1360
150
        panProjParamId[0] = ProjCenterLatGeoKey;
1361
150
        panProjParamId[1] = ProjCenterLongGeoKey;
1362
150
        panProjParamId[2] = ProjStdParallel1GeoKey;
1363
150
        panProjParamId[5] = ProjFalseEastingGeoKey;
1364
150
        panProjParamId[6] = ProjFalseNorthingGeoKey;
1365
1366
150
        panEPSGCodes[0] = EPSGNatOriginLat;
1367
150
        panEPSGCodes[1] = EPSGNatOriginLong;
1368
150
        panEPSGCodes[2] = EPSGStdParallel1Lat;
1369
150
        panEPSGCodes[5] = EPSGFalseEasting;
1370
150
        panEPSGCodes[6] = EPSGFalseNorthing;
1371
150
        return TRUE;
1372
1373
0
      case CT_Ext_Mercator_2SP:
1374
0
        panProjParamId[0] = ProjNatOriginLatGeoKey;
1375
0
        panProjParamId[1] = ProjNatOriginLongGeoKey;
1376
0
        panProjParamId[2] = ProjStdParallel1GeoKey;
1377
0
        panProjParamId[5] = ProjFalseEastingGeoKey;
1378
0
        panProjParamId[6] = ProjFalseNorthingGeoKey;
1379
1380
0
        panEPSGCodes[0] = EPSGNatOriginLat;
1381
0
        panEPSGCodes[1] = EPSGNatOriginLong;
1382
0
        panEPSGCodes[2] = EPSGStdParallel1Lat;
1383
0
        panEPSGCodes[5] = EPSGFalseEasting;
1384
0
        panEPSGCodes[6] = EPSGFalseNorthing;
1385
0
        return TRUE;
1386
1387
6.81k
      default:
1388
6.81k
        return FALSE;
1389
32.4k
    }
1390
32.4k
}
1391
1392
/************************************************************************/
1393
/*                         GTIFGetProjTRFInfo()                         */
1394
/*                                                                      */
1395
/*      Transform a PROJECTION_TRF_CODE into a projection method,       */
1396
/*      and a set of parameters.  The parameters identify will          */
1397
/*      depend on the returned method, but they will all have been      */
1398
/*      normalized into degrees and meters.                             */
1399
/************************************************************************/
1400
1401
int GTIFGetProjTRFInfoEx( void* ctxIn,
1402
                          int nProjTRFCode,
1403
                          char **ppszProjTRFName,
1404
                          short * pnProjMethod,
1405
                          double * padfProjParams )
1406
1407
23.7k
{
1408
23.7k
    if ((nProjTRFCode >= Proj_UTM_zone_1N && nProjTRFCode <= Proj_UTM_zone_60N) ||
1409
18.4k
        (nProjTRFCode >= Proj_UTM_zone_1S && nProjTRFCode <= Proj_UTM_zone_60S))
1410
8.82k
    {
1411
8.82k
        int bNorth;
1412
8.82k
        int nZone;
1413
8.82k
        if (nProjTRFCode <= Proj_UTM_zone_60N)
1414
5.28k
        {
1415
5.28k
            bNorth = TRUE;
1416
5.28k
            nZone = nProjTRFCode - Proj_UTM_zone_1N + 1;
1417
5.28k
        }
1418
3.54k
        else
1419
3.54k
        {
1420
3.54k
            bNorth = FALSE;
1421
3.54k
            nZone = nProjTRFCode - Proj_UTM_zone_1S + 1;
1422
3.54k
        }
1423
1424
8.82k
        if (ppszProjTRFName)
1425
0
        {
1426
0
            char szProjTRFName[64];
1427
0
            snprintf(szProjTRFName, sizeof(szProjTRFName), "UTM zone %d%c",
1428
0
                     nZone, (bNorth) ? 'N' : 'S');
1429
0
            *ppszProjTRFName = CPLStrdup(szProjTRFName);
1430
0
        }
1431
1432
8.82k
        if (pnProjMethod)
1433
8.82k
            *pnProjMethod = 9807;
1434
1435
8.82k
        if (padfProjParams)
1436
8.82k
        {
1437
8.82k
            padfProjParams[0] = 0;
1438
8.82k
            padfProjParams[1] = -183 + 6 * nZone;
1439
8.82k
            padfProjParams[2] = 0;
1440
8.82k
            padfProjParams[3] = 0;
1441
8.82k
            padfProjParams[4] = 0.9996;
1442
8.82k
            padfProjParams[5] = 500000;
1443
8.82k
            padfProjParams[6] = (bNorth) ? 0 : 10000000;
1444
8.82k
        }
1445
1446
8.82k
        return TRUE;
1447
8.82k
    }
1448
1449
14.9k
    if( nProjTRFCode == KvUserDefined )
1450
0
        return FALSE;
1451
1452
14.9k
    {
1453
14.9k
        char    szCode[12];
1454
14.9k
        snprintf(szCode, sizeof(szCode), "%d", nProjTRFCode);
1455
14.9k
        PJ_CONTEXT* ctx = (PJ_CONTEXT*)ctxIn;
1456
14.9k
        PJ *transf = proj_create_from_database(
1457
14.9k
            ctx, "EPSG", szCode, PJ_CATEGORY_COORDINATE_OPERATION, 0, NULL);
1458
14.9k
        if( !transf )
1459
1.45k
        {
1460
1.45k
            return FALSE;
1461
1.45k
        }
1462
1463
13.4k
        if( proj_get_type(transf) != PJ_TYPE_CONVERSION )
1464
4.80k
        {
1465
4.80k
            proj_destroy(transf);
1466
4.80k
            return FALSE;
1467
4.80k
        }
1468
1469
        /* Get the projection method code */
1470
8.68k
        const char* pszMethodCode = NULL;
1471
8.68k
        proj_coordoperation_get_method_info(ctx, transf,
1472
8.68k
                                            NULL, /* method name */
1473
8.68k
                                            NULL, /* method auth name (should be EPSG) */
1474
8.68k
                                            &pszMethodCode);
1475
8.68k
        assert( pszMethodCode );
1476
8.68k
        const int nProjMethod = atoi(pszMethodCode);
1477
1478
/* -------------------------------------------------------------------- */
1479
/*      Initialize a definition of what EPSG codes need to be loaded    */
1480
/*      into what fields in adfProjParams.                               */
1481
/* -------------------------------------------------------------------- */
1482
8.68k
        const int nCTProjMethod = EPSGProjMethodToCTProjMethod( nProjMethod, TRUE );
1483
8.68k
        int anEPSGCodes[7];
1484
8.68k
        SetGTParamIds( nCTProjMethod, nProjMethod, NULL, anEPSGCodes );
1485
1486
1487
/* -------------------------------------------------------------------- */
1488
/*      Get the parameters for this projection.                         */
1489
/* -------------------------------------------------------------------- */
1490
1491
8.68k
        double  adfProjParams[7];
1492
1493
69.4k
        for( int i = 0; i < 7; i++ )
1494
60.7k
        {
1495
60.7k
            int nEPSGCode = anEPSGCodes[i];
1496
1497
            /* Establish default */
1498
60.7k
            if( nEPSGCode == EPSGAngleRectifiedToSkewedGrid )
1499
129
                adfProjParams[i] = 90.0;
1500
60.6k
            else if( nEPSGCode == EPSGNatOriginScaleFactor
1501
54.8k
                    || nEPSGCode == EPSGInitialLineScaleFactor
1502
54.7k
                    || nEPSGCode == EPSGPseudoStdParallelScaleFactor )
1503
5.91k
                adfProjParams[i] = 1.0;
1504
54.7k
            else
1505
54.7k
                adfProjParams[i] = 0.0;
1506
1507
            /* If there is no parameter, skip */
1508
60.7k
            if( nEPSGCode == 0 )
1509
16.4k
                continue;
1510
1511
44.2k
            const int nParamCount = proj_coordoperation_get_param_count(ctx, transf);
1512
1513
            /* Find the matching parameter */
1514
44.2k
            const char *pszUOMCategory = NULL;
1515
44.2k
            double  dfValue = 0.0;
1516
44.2k
            double  dfUnitConvFactor = 0.0;
1517
44.2k
            int iEPSG;  /* Used after for */
1518
140k
            for( iEPSG = 0; iEPSG < nParamCount; iEPSG++ )
1519
140k
            {
1520
140k
                const char* pszParamCode = NULL;
1521
140k
                proj_coordoperation_get_param(
1522
140k
                    ctx, transf, iEPSG,
1523
140k
                    NULL, /* name */
1524
140k
                    NULL, /* auth name */
1525
140k
                    &pszParamCode,
1526
140k
                    &dfValue,
1527
140k
                    NULL, /* value (string) */
1528
140k
                    &dfUnitConvFactor, /* unit conv factor */
1529
140k
                    NULL, /* unit name */
1530
140k
                    NULL, /* unit auth name */
1531
140k
                    NULL, /* unit code */
1532
140k
                    &pszUOMCategory /* unit category */);
1533
140k
                assert(pszParamCode);
1534
140k
                if( atoi(pszParamCode) == nEPSGCode )
1535
43.9k
                {
1536
43.9k
                    break;
1537
43.9k
                }
1538
140k
            }
1539
1540
            /* not found, accept the default */
1541
44.2k
            if( iEPSG == nParamCount )
1542
380
            {
1543
                /* for CT_ObliqueMercator try alternate parameter codes first */
1544
                /* because EPSG proj method 9812 uses EPSGFalseXXXXX, but 9815 uses EPSGProjCenterXXXXX */
1545
380
                if ( nCTProjMethod == CT_ObliqueMercator && nEPSGCode == EPSGProjCenterEasting )
1546
33
                    nEPSGCode = EPSGFalseEasting;
1547
347
                else if ( nCTProjMethod == CT_ObliqueMercator && nEPSGCode == EPSGProjCenterNorthing )
1548
33
                    nEPSGCode = EPSGFalseNorthing;
1549
                /* for CT_PolarStereographic try alternate parameter codes first */
1550
                /* because EPSG proj method 9829 uses EPSGLatOfStdParallel instead of EPSGNatOriginLat */
1551
                /* and EPSGOriginLong instead of EPSGNatOriginLong */
1552
314
                else if( nCTProjMethod == CT_PolarStereographic && nEPSGCode == EPSGNatOriginLat )
1553
54
                    nEPSGCode = EPSGLatOfStdParallel;
1554
260
                else if( nCTProjMethod == CT_PolarStereographic && nEPSGCode == EPSGNatOriginLong )
1555
54
                    nEPSGCode = EPSGOriginLong;
1556
206
                else
1557
206
                    continue;
1558
1559
591
                for( iEPSG = 0; iEPSG < nParamCount; iEPSG++ )
1560
591
                {
1561
591
                    const char* pszParamCode = NULL;
1562
591
                    proj_coordoperation_get_param(
1563
591
                        ctx, transf, iEPSG,
1564
591
                        NULL, /* name */
1565
591
                        NULL, /* auth name */
1566
591
                        &pszParamCode,
1567
591
                        &dfValue,
1568
591
                        NULL, /* value (string) */
1569
591
                        &dfUnitConvFactor, /* unit conv factor */
1570
591
                        NULL, /* unit name */
1571
591
                        NULL, /* unit auth name */
1572
591
                        NULL, /* unit code */
1573
591
                        &pszUOMCategory /* unit category */);
1574
591
                    assert(pszParamCode);
1575
591
                    if( atoi(pszParamCode) == nEPSGCode )
1576
174
                    {
1577
174
                        break;
1578
174
                    }
1579
591
                }
1580
1581
174
                if( iEPSG == nParamCount )
1582
0
                    continue;
1583
174
            }
1584
1585
44.2k
            assert(pszUOMCategory);
1586
1587
44.0k
            adfProjParams[i] = dfValue * dfUnitConvFactor;
1588
44.0k
            if( strcmp(pszUOMCategory, "angular") == 0.0 )
1589
21.4k
            {
1590
                /* Convert from radians to degrees */
1591
21.4k
                adfProjParams[i] *= 180 / M_PI;
1592
21.4k
            }
1593
44.0k
        }
1594
1595
/* -------------------------------------------------------------------- */
1596
/*      Get the name, if requested.                                     */
1597
/* -------------------------------------------------------------------- */
1598
8.68k
        if( ppszProjTRFName != NULL )
1599
0
        {
1600
0
            const char* pszName = proj_get_name(transf);
1601
0
            if( !pszName )
1602
0
            {
1603
                // shouldn't happen
1604
0
                proj_destroy(transf);
1605
0
                return FALSE;
1606
0
            }
1607
0
            *ppszProjTRFName = CPLStrdup(pszName);
1608
0
        }
1609
1610
/* -------------------------------------------------------------------- */
1611
/*      Transfer requested data into passed variables.                  */
1612
/* -------------------------------------------------------------------- */
1613
8.68k
        if( pnProjMethod != NULL )
1614
8.68k
            *pnProjMethod = (short) nProjMethod;
1615
1616
8.68k
        if( padfProjParams != NULL )
1617
8.68k
        {
1618
69.4k
            for( int i = 0; i < 7; i++ )
1619
60.7k
                padfProjParams[i] = adfProjParams[i];
1620
8.68k
        }
1621
1622
8.68k
        proj_destroy(transf);
1623
1624
8.68k
        return TRUE;
1625
8.68k
    }
1626
8.68k
}
1627
1628
int GTIFGetProjTRFInfo( /* Conversion code */
1629
                        int nProjTRFCode,
1630
                        char **ppszProjTRFName,
1631
                        short * pnProjMethod,
1632
                        double * padfProjParams )
1633
0
{
1634
0
    PJ_CONTEXT* ctx = proj_context_create();
1635
0
    const int ret = GTIFGetProjTRFInfoEx(
1636
0
        ctx, nProjTRFCode, ppszProjTRFName, pnProjMethod, padfProjParams);
1637
0
    proj_context_destroy(ctx);
1638
0
    return ret;
1639
0
}
1640
1641
/************************************************************************/
1642
/*                         GTIFFetchProjParms()                         */
1643
/*                                                                      */
1644
/*      Fetch the projection parameters for a particular projection     */
1645
/*      from a GeoTIFF file, and fill the GTIFDefn structure out        */
1646
/*      with them.                                                      */
1647
/************************************************************************/
1648
1649
static void GTIFFetchProjParms( GTIF * psGTIF, GTIFDefn * psDefn )
1650
10.9k
{
1651
1652
/* -------------------------------------------------------------------- */
1653
/*      Get the false easting, and northing if available.               */
1654
/* -------------------------------------------------------------------- */
1655
10.9k
    double dfFalseEasting = 0.0;
1656
10.9k
    if( !GTIFKeyGetDOUBLE(psGTIF, ProjFalseEastingGeoKey, &dfFalseEasting, 0, 1)
1657
9.29k
        && !GTIFKeyGetDOUBLE(psGTIF, ProjCenterEastingGeoKey,
1658
9.29k
                       &dfFalseEasting, 0, 1)
1659
9.27k
        && !GTIFKeyGetDOUBLE(psGTIF, ProjFalseOriginEastingGeoKey,
1660
9.27k
                       &dfFalseEasting, 0, 1) )
1661
9.25k
        dfFalseEasting = 0.0;
1662
1663
10.9k
    double dfFalseNorthing = 0.0;
1664
10.9k
    if( !GTIFKeyGetDOUBLE(psGTIF, ProjFalseNorthingGeoKey, &dfFalseNorthing,0,1)
1665
9.02k
        && !GTIFKeyGetDOUBLE(psGTIF, ProjCenterNorthingGeoKey,
1666
9.02k
                       &dfFalseNorthing, 0, 1)
1667
8.99k
        && !GTIFKeyGetDOUBLE(psGTIF, ProjFalseOriginNorthingGeoKey,
1668
8.99k
                       &dfFalseNorthing, 0, 1) )
1669
8.90k
        dfFalseNorthing = 0.0;
1670
1671
10.9k
    double dfNatOriginLong = 0.0, dfNatOriginLat = 0.0, dfRectGridAngle = 0.0;
1672
10.9k
    double dfNatOriginScale = 1.0;
1673
10.9k
    double dfStdParallel1 = 0.0, dfStdParallel2 = 0.0, dfAzimuth = 0.0;
1674
1675
10.9k
    switch( psDefn->CTProjection )
1676
10.9k
    {
1677
/* -------------------------------------------------------------------- */
1678
180
      case CT_Stereographic:
1679
/* -------------------------------------------------------------------- */
1680
180
        if( GTIFKeyGetDOUBLE(psGTIF, ProjNatOriginLongGeoKey,
1681
180
                       &dfNatOriginLong, 0, 1 ) == 0
1682
173
            && GTIFKeyGetDOUBLE(psGTIF, ProjFalseOriginLongGeoKey,
1683
173
                          &dfNatOriginLong, 0, 1 ) == 0
1684
170
            && GTIFKeyGetDOUBLE(psGTIF, ProjCenterLongGeoKey,
1685
170
                          &dfNatOriginLong, 0, 1 ) == 0 )
1686
167
            dfNatOriginLong = 0.0;
1687
1688
180
        if( GTIFKeyGetDOUBLE(psGTIF, ProjNatOriginLatGeoKey,
1689
180
                       &dfNatOriginLat, 0, 1 ) == 0
1690
170
            && GTIFKeyGetDOUBLE(psGTIF, ProjFalseOriginLatGeoKey,
1691
170
                          &dfNatOriginLat, 0, 1 ) == 0
1692
168
            && GTIFKeyGetDOUBLE(psGTIF, ProjCenterLatGeoKey,
1693
168
                          &dfNatOriginLat, 0, 1 ) == 0 )
1694
167
            dfNatOriginLat = 0.0;
1695
1696
180
        if( GTIFKeyGetDOUBLE(psGTIF, ProjScaleAtNatOriginGeoKey,
1697
180
                       &dfNatOriginScale, 0, 1 ) == 0 )
1698
169
            dfNatOriginScale = 1.0;
1699
1700
        /* notdef: should transform to decimal degrees at this point */
1701
1702
180
        psDefn->ProjParm[0] = dfNatOriginLat;
1703
180
        psDefn->ProjParmId[0] = ProjCenterLatGeoKey;
1704
180
        psDefn->ProjParm[1] = dfNatOriginLong;
1705
180
        psDefn->ProjParmId[1] = ProjCenterLongGeoKey;
1706
180
        psDefn->ProjParm[4] = dfNatOriginScale;
1707
180
        psDefn->ProjParmId[4] = ProjScaleAtNatOriginGeoKey;
1708
180
        psDefn->ProjParm[5] = dfFalseEasting;
1709
180
        psDefn->ProjParmId[5] = ProjFalseEastingGeoKey;
1710
180
        psDefn->ProjParm[6] = dfFalseNorthing;
1711
180
        psDefn->ProjParmId[6] = ProjFalseNorthingGeoKey;
1712
1713
180
        psDefn->nParms = 7;
1714
180
        break;
1715
1716
/* -------------------------------------------------------------------- */
1717
2.14k
      case CT_Mercator:
1718
/* -------------------------------------------------------------------- */
1719
2.14k
        if( GTIFKeyGetDOUBLE(psGTIF, ProjNatOriginLongGeoKey,
1720
2.14k
                       &dfNatOriginLong, 0, 1 ) == 0
1721
1.81k
            && GTIFKeyGetDOUBLE(psGTIF, ProjFalseOriginLongGeoKey,
1722
1.81k
                          &dfNatOriginLong, 0, 1 ) == 0
1723
1.81k
            && GTIFKeyGetDOUBLE(psGTIF, ProjCenterLongGeoKey,
1724
1.81k
                          &dfNatOriginLong, 0, 1 ) == 0 )
1725
1.78k
            dfNatOriginLong = 0.0;
1726
1727
2.14k
        if( GTIFKeyGetDOUBLE(psGTIF, ProjNatOriginLatGeoKey,
1728
2.14k
                       &dfNatOriginLat, 0, 1 ) == 0
1729
590
            && GTIFKeyGetDOUBLE(psGTIF, ProjFalseOriginLatGeoKey,
1730
590
                          &dfNatOriginLat, 0, 1 ) == 0
1731
548
            && GTIFKeyGetDOUBLE(psGTIF, ProjCenterLatGeoKey,
1732
548
                          &dfNatOriginLat, 0, 1 ) == 0 )
1733
542
            dfNatOriginLat = 0.0;
1734
1735
1736
2.14k
        const int bHaveSP1 = GTIFKeyGetDOUBLE(psGTIF, ProjStdParallel1GeoKey,
1737
2.14k
                              &dfStdParallel1, 0, 1 );
1738
1739
2.14k
        int bHaveNOS = GTIFKeyGetDOUBLE(psGTIF, ProjScaleAtNatOriginGeoKey,
1740
2.14k
                              &dfNatOriginScale, 0, 1 );
1741
1742
        /* Default scale only if dfStdParallel1 isn't defined either */
1743
2.14k
        if( !bHaveNOS && !bHaveSP1)
1744
649
        {
1745
649
            bHaveNOS = TRUE;
1746
649
            dfNatOriginScale = 1.0;
1747
649
        }
1748
1749
        /* notdef: should transform to decimal degrees at this point */
1750
1751
2.14k
        psDefn->ProjParm[0] = dfNatOriginLat;
1752
2.14k
        psDefn->ProjParmId[0] = ProjNatOriginLatGeoKey;
1753
2.14k
        psDefn->ProjParm[1] = dfNatOriginLong;
1754
2.14k
        psDefn->ProjParmId[1] = ProjNatOriginLongGeoKey;
1755
2.14k
        if( bHaveSP1 )
1756
1.18k
        {
1757
1.18k
            psDefn->ProjParm[2] = dfStdParallel1;
1758
1.18k
            psDefn->ProjParmId[2] = ProjStdParallel1GeoKey;
1759
1.18k
        }
1760
2.14k
        if( bHaveNOS )
1761
1.06k
        {
1762
1.06k
            psDefn->ProjParm[4] = dfNatOriginScale;
1763
1.06k
            psDefn->ProjParmId[4] = ProjScaleAtNatOriginGeoKey;
1764
1.06k
        }
1765
2.14k
        psDefn->ProjParm[5] = dfFalseEasting;
1766
2.14k
        psDefn->ProjParmId[5] = ProjFalseEastingGeoKey;
1767
2.14k
        psDefn->ProjParm[6] = dfFalseNorthing;
1768
2.14k
        psDefn->ProjParmId[6] = ProjFalseNorthingGeoKey;
1769
1770
2.14k
        psDefn->nParms = 7;
1771
2.14k
        break;
1772
1773
/* -------------------------------------------------------------------- */
1774
519
      case CT_LambertConfConic_1SP:
1775
546
      case CT_ObliqueStereographic:
1776
706
      case CT_TransverseMercator:
1777
735
      case CT_TransvMercator_SouthOriented:
1778
/* -------------------------------------------------------------------- */
1779
735
        if( GTIFKeyGetDOUBLE(psGTIF, ProjNatOriginLongGeoKey,
1780
735
                       &dfNatOriginLong, 0, 1 ) == 0
1781
670
            && GTIFKeyGetDOUBLE(psGTIF, ProjFalseOriginLongGeoKey,
1782
670
                          &dfNatOriginLong, 0, 1 ) == 0
1783
667
            && GTIFKeyGetDOUBLE(psGTIF, ProjCenterLongGeoKey,
1784
667
                          &dfNatOriginLong, 0, 1 ) == 0 )
1785
563
            dfNatOriginLong = 0.0;
1786
1787
735
        if( GTIFKeyGetDOUBLE(psGTIF, ProjNatOriginLatGeoKey,
1788
735
                       &dfNatOriginLat, 0, 1 ) == 0
1789
665
            && GTIFKeyGetDOUBLE(psGTIF, ProjFalseOriginLatGeoKey,
1790
665
                          &dfNatOriginLat, 0, 1 ) == 0
1791
659
            && GTIFKeyGetDOUBLE(psGTIF, ProjCenterLatGeoKey,
1792
659
                          &dfNatOriginLat, 0, 1 ) == 0 )
1793
654
            dfNatOriginLat = 0.0;
1794
1795
735
        if( GTIFKeyGetDOUBLE(psGTIF, ProjScaleAtNatOriginGeoKey,
1796
735
                       &dfNatOriginScale, 0, 1 ) == 0
1797
            /* See https://github.com/OSGeo/gdal/files/1665718/lasinfo.txt */
1798
657
            && GTIFKeyGetDOUBLE(psGTIF, ProjScaleAtCenterGeoKey,
1799
657
                       &dfNatOriginScale, 0, 1 ) == 0 )
1800
652
            dfNatOriginScale = 1.0;
1801
1802
        /* notdef: should transform to decimal degrees at this point */
1803
1804
735
        psDefn->ProjParm[0] = dfNatOriginLat;
1805
735
        psDefn->ProjParmId[0] = ProjNatOriginLatGeoKey;
1806
735
        psDefn->ProjParm[1] = dfNatOriginLong;
1807
735
        psDefn->ProjParmId[1] = ProjNatOriginLongGeoKey;
1808
735
        psDefn->ProjParm[4] = dfNatOriginScale;
1809
735
        psDefn->ProjParmId[4] = ProjScaleAtNatOriginGeoKey;
1810
735
        psDefn->ProjParm[5] = dfFalseEasting;
1811
735
        psDefn->ProjParmId[5] = ProjFalseEastingGeoKey;
1812
735
        psDefn->ProjParm[6] = dfFalseNorthing;
1813
735
        psDefn->ProjParmId[6] = ProjFalseNorthingGeoKey;
1814
1815
735
        psDefn->nParms = 7;
1816
735
        break;
1817
1818
/* -------------------------------------------------------------------- */
1819
915
      case CT_ObliqueMercator: /* hotine */
1820
915
      case CT_HotineObliqueMercatorAzimuthCenter:
1821
/* -------------------------------------------------------------------- */
1822
915
        if( GTIFKeyGetDOUBLE(psGTIF, ProjNatOriginLongGeoKey,
1823
915
                       &dfNatOriginLong, 0, 1 ) == 0
1824
904
            && GTIFKeyGetDOUBLE(psGTIF, ProjFalseOriginLongGeoKey,
1825
904
                          &dfNatOriginLong, 0, 1 ) == 0
1826
903
            && GTIFKeyGetDOUBLE(psGTIF, ProjCenterLongGeoKey,
1827
903
                          &dfNatOriginLong, 0, 1 ) == 0 )
1828
901
            dfNatOriginLong = 0.0;
1829
1830
915
        if( GTIFKeyGetDOUBLE(psGTIF, ProjNatOriginLatGeoKey,
1831
915
                       &dfNatOriginLat, 0, 1 ) == 0
1832
835
            && GTIFKeyGetDOUBLE(psGTIF, ProjFalseOriginLatGeoKey,
1833
835
                          &dfNatOriginLat, 0, 1 ) == 0
1834
831
            && GTIFKeyGetDOUBLE(psGTIF, ProjCenterLatGeoKey,
1835
831
                          &dfNatOriginLat, 0, 1 ) == 0 )
1836
830
            dfNatOriginLat = 0.0;
1837
1838
915
        if( GTIFKeyGetDOUBLE(psGTIF, ProjAzimuthAngleGeoKey,
1839
915
                       &dfAzimuth, 0, 1 ) == 0 )
1840
913
            dfAzimuth = 0.0;
1841
1842
915
        if( GTIFKeyGetDOUBLE(psGTIF, ProjRectifiedGridAngleGeoKey,
1843
915
                       &dfRectGridAngle, 0, 1 ) == 0 )
1844
915
            dfRectGridAngle = 90.0;
1845
1846
915
        if( GTIFKeyGetDOUBLE(psGTIF, ProjScaleAtNatOriginGeoKey,
1847
915
                       &dfNatOriginScale, 0, 1 ) == 0
1848
883
            && GTIFKeyGetDOUBLE(psGTIF, ProjScaleAtCenterGeoKey,
1849
883
                          &dfNatOriginScale, 0, 1 ) == 0 )
1850
881
            dfNatOriginScale = 1.0;
1851
1852
        /* notdef: should transform to decimal degrees at this point */
1853
1854
915
        psDefn->ProjParm[0] = dfNatOriginLat;
1855
915
        psDefn->ProjParmId[0] = ProjCenterLatGeoKey;
1856
915
        psDefn->ProjParm[1] = dfNatOriginLong;
1857
915
        psDefn->ProjParmId[1] = ProjCenterLongGeoKey;
1858
915
        psDefn->ProjParm[2] = dfAzimuth;
1859
915
        psDefn->ProjParmId[2] = ProjAzimuthAngleGeoKey;
1860
915
        psDefn->ProjParm[3] = dfRectGridAngle;
1861
915
        psDefn->ProjParmId[3] = ProjRectifiedGridAngleGeoKey;
1862
915
        psDefn->ProjParm[4] = dfNatOriginScale;
1863
915
        psDefn->ProjParmId[4] = ProjScaleAtCenterGeoKey;
1864
915
        psDefn->ProjParm[5] = dfFalseEasting;
1865
915
        psDefn->ProjParmId[5] = ProjFalseEastingGeoKey;
1866
915
        psDefn->ProjParm[6] = dfFalseNorthing;
1867
915
        psDefn->ProjParmId[6] = ProjFalseNorthingGeoKey;
1868
1869
915
        psDefn->nParms = 7;
1870
915
        break;
1871
1872
/* -------------------------------------------------------------------- */
1873
595
      case CT_ObliqueMercator_Laborde:
1874
/* -------------------------------------------------------------------- */
1875
595
        if( GTIFKeyGetDOUBLE(psGTIF, ProjNatOriginLongGeoKey,
1876
595
                       &dfNatOriginLong, 0, 1 ) == 0
1877
586
            && GTIFKeyGetDOUBLE(psGTIF, ProjFalseOriginLongGeoKey,
1878
586
                          &dfNatOriginLong, 0, 1 ) == 0
1879
586
            && GTIFKeyGetDOUBLE(psGTIF, ProjCenterLongGeoKey,
1880
586
                          &dfNatOriginLong, 0, 1 ) == 0 )
1881
585
            dfNatOriginLong = 0.0;
1882
1883
595
        if( GTIFKeyGetDOUBLE(psGTIF, ProjNatOriginLatGeoKey,
1884
595
                       &dfNatOriginLat, 0, 1 ) == 0
1885
583
            && GTIFKeyGetDOUBLE(psGTIF, ProjFalseOriginLatGeoKey,
1886
583
                          &dfNatOriginLat, 0, 1 ) == 0
1887
579
            && GTIFKeyGetDOUBLE(psGTIF, ProjCenterLatGeoKey,
1888
579
                          &dfNatOriginLat, 0, 1 ) == 0 )
1889
577
            dfNatOriginLat = 0.0;
1890
1891
595
        if( GTIFKeyGetDOUBLE(psGTIF, ProjAzimuthAngleGeoKey,
1892
595
                       &dfAzimuth, 0, 1 ) == 0 )
1893
593
            dfAzimuth = 0.0;
1894
1895
595
        if( GTIFKeyGetDOUBLE(psGTIF, ProjScaleAtNatOriginGeoKey,
1896
595
                       &dfNatOriginScale, 0, 1 ) == 0
1897
589
            && GTIFKeyGetDOUBLE(psGTIF, ProjScaleAtCenterGeoKey,
1898
589
                          &dfNatOriginScale, 0, 1 ) == 0 )
1899
588
            dfNatOriginScale = 1.0;
1900
1901
        /* notdef: should transform to decimal degrees at this point */
1902
1903
595
        psDefn->ProjParm[0] = dfNatOriginLat;
1904
595
        psDefn->ProjParmId[0] = ProjCenterLatGeoKey;
1905
595
        psDefn->ProjParm[1] = dfNatOriginLong;
1906
595
        psDefn->ProjParmId[1] = ProjCenterLongGeoKey;
1907
595
        psDefn->ProjParm[2] = dfAzimuth;
1908
595
        psDefn->ProjParmId[2] = ProjAzimuthAngleGeoKey;
1909
595
        psDefn->ProjParm[4] = dfNatOriginScale;
1910
595
        psDefn->ProjParmId[4] = ProjScaleAtCenterGeoKey;
1911
595
        psDefn->ProjParm[5] = dfFalseEasting;
1912
595
        psDefn->ProjParmId[5] = ProjFalseEastingGeoKey;
1913
595
        psDefn->ProjParm[6] = dfFalseNorthing;
1914
595
        psDefn->ProjParmId[6] = ProjFalseNorthingGeoKey;
1915
1916
595
        psDefn->nParms = 7;
1917
595
        break;
1918
1919
/* -------------------------------------------------------------------- */
1920
80
      case CT_CassiniSoldner:
1921
98
      case CT_Polyconic:
1922
/* -------------------------------------------------------------------- */
1923
98
        if( GTIFKeyGetDOUBLE(psGTIF, ProjNatOriginLongGeoKey,
1924
98
                       &dfNatOriginLong, 0, 1 ) == 0
1925
92
            && GTIFKeyGetDOUBLE(psGTIF, ProjFalseOriginLongGeoKey,
1926
92
                          &dfNatOriginLong, 0, 1 ) == 0
1927
90
            && GTIFKeyGetDOUBLE(psGTIF, ProjCenterLongGeoKey,
1928
90
                          &dfNatOriginLong, 0, 1 ) == 0 )
1929
87
            dfNatOriginLong = 0.0;
1930
1931
98
        if( GTIFKeyGetDOUBLE(psGTIF, ProjNatOriginLatGeoKey,
1932
98
                       &dfNatOriginLat, 0, 1 ) == 0
1933
87
            && GTIFKeyGetDOUBLE(psGTIF, ProjFalseOriginLatGeoKey,
1934
87
                          &dfNatOriginLat, 0, 1 ) == 0
1935
77
            && GTIFKeyGetDOUBLE(psGTIF, ProjCenterLatGeoKey,
1936
77
                          &dfNatOriginLat, 0, 1 ) == 0 )
1937
75
            dfNatOriginLat = 0.0;
1938
1939
98
        if( GTIFKeyGetDOUBLE(psGTIF, ProjScaleAtNatOriginGeoKey,
1940
98
                       &dfNatOriginScale, 0, 1 ) == 0
1941
95
            && GTIFKeyGetDOUBLE(psGTIF, ProjScaleAtCenterGeoKey,
1942
95
                          &dfNatOriginScale, 0, 1 ) == 0 )
1943
94
            dfNatOriginScale = 1.0;
1944
1945
        /* notdef: should transform to decimal degrees at this point */
1946
1947
98
        psDefn->ProjParm[0] = dfNatOriginLat;
1948
98
        psDefn->ProjParmId[0] = ProjNatOriginLatGeoKey;
1949
98
        psDefn->ProjParm[1] = dfNatOriginLong;
1950
98
        psDefn->ProjParmId[1] = ProjNatOriginLongGeoKey;
1951
98
        psDefn->ProjParm[4] = dfNatOriginScale;
1952
98
        psDefn->ProjParmId[4] = ProjScaleAtNatOriginGeoKey;
1953
98
        psDefn->ProjParm[5] = dfFalseEasting;
1954
98
        psDefn->ProjParmId[5] = ProjFalseEastingGeoKey;
1955
98
        psDefn->ProjParm[6] = dfFalseNorthing;
1956
98
        psDefn->ProjParmId[6] = ProjFalseNorthingGeoKey;
1957
1958
98
        psDefn->nParms = 7;
1959
98
        break;
1960
1961
/* -------------------------------------------------------------------- */
1962
2.13k
      case CT_AzimuthalEquidistant:
1963
2.18k
      case CT_MillerCylindrical:
1964
2.88k
      case CT_Gnomonic:
1965
3.31k
      case CT_LambertAzimEqualArea:
1966
3.33k
      case CT_Orthographic:
1967
3.47k
      case CT_NewZealandMapGrid:
1968
/* -------------------------------------------------------------------- */
1969
3.47k
        if( GTIFKeyGetDOUBLE(psGTIF, ProjNatOriginLongGeoKey,
1970
3.47k
                       &dfNatOriginLong, 0, 1 ) == 0
1971
3.43k
            && GTIFKeyGetDOUBLE(psGTIF, ProjFalseOriginLongGeoKey,
1972
3.43k
                          &dfNatOriginLong, 0, 1 ) == 0
1973
3.43k
            && GTIFKeyGetDOUBLE(psGTIF, ProjCenterLongGeoKey,
1974
3.43k
                          &dfNatOriginLong, 0, 1 ) == 0 )
1975
3.41k
            dfNatOriginLong = 0.0;
1976
1977
3.47k
        if( GTIFKeyGetDOUBLE(psGTIF, ProjNatOriginLatGeoKey,
1978
3.47k
                       &dfNatOriginLat, 0, 1 ) == 0
1979
3.36k
            && GTIFKeyGetDOUBLE(psGTIF, ProjFalseOriginLatGeoKey,
1980
3.36k
                          &dfNatOriginLat, 0, 1 ) == 0
1981
3.34k
            && GTIFKeyGetDOUBLE(psGTIF, ProjCenterLatGeoKey,
1982
3.34k
                          &dfNatOriginLat, 0, 1 ) == 0 )
1983
3.32k
            dfNatOriginLat = 0.0;
1984
1985
        /* notdef: should transform to decimal degrees at this point */
1986
1987
3.47k
        psDefn->ProjParm[0] = dfNatOriginLat;
1988
3.47k
        psDefn->ProjParmId[0] = ProjCenterLatGeoKey;
1989
3.47k
        psDefn->ProjParm[1] = dfNatOriginLong;
1990
3.47k
        psDefn->ProjParmId[1] = ProjCenterLongGeoKey;
1991
3.47k
        psDefn->ProjParm[5] = dfFalseEasting;
1992
3.47k
        psDefn->ProjParmId[5] = ProjFalseEastingGeoKey;
1993
3.47k
        psDefn->ProjParm[6] = dfFalseNorthing;
1994
3.47k
        psDefn->ProjParmId[6] = ProjFalseNorthingGeoKey;
1995
1996
3.47k
        psDefn->nParms = 7;
1997
3.47k
        break;
1998
1999
/* -------------------------------------------------------------------- */
2000
1.12k
      case CT_Equirectangular:
2001
/* -------------------------------------------------------------------- */
2002
1.12k
        if( GTIFKeyGetDOUBLE(psGTIF, ProjNatOriginLongGeoKey,
2003
1.12k
                       &dfNatOriginLong, 0, 1 ) == 0
2004
1.00k
            && GTIFKeyGetDOUBLE(psGTIF, ProjFalseOriginLongGeoKey,
2005
1.00k
                          &dfNatOriginLong, 0, 1 ) == 0
2006
996
            && GTIFKeyGetDOUBLE(psGTIF, ProjCenterLongGeoKey,
2007
996
                          &dfNatOriginLong, 0, 1 ) == 0 )
2008
988
            dfNatOriginLong = 0.0;
2009
2010
1.12k
        if( GTIFKeyGetDOUBLE(psGTIF, ProjNatOriginLatGeoKey,
2011
1.12k
                       &dfNatOriginLat, 0, 1 ) == 0
2012
831
            && GTIFKeyGetDOUBLE(psGTIF, ProjFalseOriginLatGeoKey,
2013
831
                          &dfNatOriginLat, 0, 1 ) == 0
2014
657
            && GTIFKeyGetDOUBLE(psGTIF, ProjCenterLatGeoKey,
2015
657
                          &dfNatOriginLat, 0, 1 ) == 0 )
2016
271
            dfNatOriginLat = 0.0;
2017
2018
1.12k
        if( GTIFKeyGetDOUBLE(psGTIF, ProjStdParallel1GeoKey,
2019
1.12k
                       &dfStdParallel1, 0, 1 ) == 0 )
2020
1.11k
            dfStdParallel1 = 0.0;
2021
2022
        /* notdef: should transform to decimal degrees at this point */
2023
2024
1.12k
        psDefn->ProjParm[0] = dfNatOriginLat;
2025
1.12k
        psDefn->ProjParmId[0] = ProjCenterLatGeoKey;
2026
1.12k
        psDefn->ProjParm[1] = dfNatOriginLong;
2027
1.12k
        psDefn->ProjParmId[1] = ProjCenterLongGeoKey;
2028
1.12k
        psDefn->ProjParm[2] = dfStdParallel1;
2029
1.12k
        psDefn->ProjParmId[2] = ProjStdParallel1GeoKey;
2030
1.12k
        psDefn->ProjParm[5] = dfFalseEasting;
2031
1.12k
        psDefn->ProjParmId[5] = ProjFalseEastingGeoKey;
2032
1.12k
        psDefn->ProjParm[6] = dfFalseNorthing;
2033
1.12k
        psDefn->ProjParmId[6] = ProjFalseNorthingGeoKey;
2034
2035
1.12k
        psDefn->nParms = 7;
2036
1.12k
        break;
2037
2038
/* -------------------------------------------------------------------- */
2039
50
      case CT_Robinson:
2040
98
      case CT_Sinusoidal:
2041
114
      case CT_VanDerGrinten:
2042
/* -------------------------------------------------------------------- */
2043
114
        if( GTIFKeyGetDOUBLE(psGTIF, ProjNatOriginLongGeoKey,
2044
114
                       &dfNatOriginLong, 0, 1 ) == 0
2045
98
            && GTIFKeyGetDOUBLE(psGTIF, ProjFalseOriginLongGeoKey,
2046
98
                          &dfNatOriginLong, 0, 1 ) == 0
2047
97
            && GTIFKeyGetDOUBLE(psGTIF, ProjCenterLongGeoKey,
2048
97
                          &dfNatOriginLong, 0, 1 ) == 0 )
2049
97
            dfNatOriginLong = 0.0;
2050
2051
        /* notdef: should transform to decimal degrees at this point */
2052
2053
114
        psDefn->ProjParm[1] = dfNatOriginLong;
2054
114
        psDefn->ProjParmId[1] = ProjCenterLongGeoKey;
2055
114
        psDefn->ProjParm[5] = dfFalseEasting;
2056
114
        psDefn->ProjParmId[5] = ProjFalseEastingGeoKey;
2057
114
        psDefn->ProjParm[6] = dfFalseNorthing;
2058
114
        psDefn->ProjParmId[6] = ProjFalseNorthingGeoKey;
2059
2060
114
        psDefn->nParms = 7;
2061
114
        break;
2062
2063
/* -------------------------------------------------------------------- */
2064
136
      case CT_PolarStereographic:
2065
/* -------------------------------------------------------------------- */
2066
136
        if( GTIFKeyGetDOUBLE(psGTIF, ProjStraightVertPoleLongGeoKey,
2067
136
                       &dfNatOriginLong, 0, 1 ) == 0
2068
131
            && GTIFKeyGetDOUBLE(psGTIF, ProjNatOriginLongGeoKey,
2069
131
                          &dfNatOriginLong, 0, 1 ) == 0
2070
85
            && GTIFKeyGetDOUBLE(psGTIF, ProjFalseOriginLongGeoKey,
2071
85
                          &dfNatOriginLong, 0, 1 ) == 0
2072
79
            && GTIFKeyGetDOUBLE(psGTIF, ProjCenterLongGeoKey,
2073
79
                          &dfNatOriginLong, 0, 1 ) == 0 )
2074
70
            dfNatOriginLong = 0.0;
2075
2076
136
        if( GTIFKeyGetDOUBLE(psGTIF, ProjNatOriginLatGeoKey,
2077
136
                       &dfNatOriginLat, 0, 1 ) == 0
2078
63
            && GTIFKeyGetDOUBLE(psGTIF, ProjFalseOriginLatGeoKey,
2079
63
                          &dfNatOriginLat, 0, 1 ) == 0
2080
56
            && GTIFKeyGetDOUBLE(psGTIF, ProjCenterLatGeoKey,
2081
56
                          &dfNatOriginLat, 0, 1 ) == 0 )
2082
55
            dfNatOriginLat = 0.0;
2083
2084
136
        if( GTIFKeyGetDOUBLE(psGTIF, ProjScaleAtNatOriginGeoKey,
2085
136
                       &dfNatOriginScale, 0, 1 ) == 0
2086
85
            && GTIFKeyGetDOUBLE(psGTIF, ProjScaleAtCenterGeoKey,
2087
85
                          &dfNatOriginScale, 0, 1 ) == 0 )
2088
83
            dfNatOriginScale = 1.0;
2089
2090
        /* notdef: should transform to decimal degrees at this point */
2091
2092
136
        psDefn->ProjParm[0] = dfNatOriginLat;
2093
136
        psDefn->ProjParmId[0] = ProjNatOriginLatGeoKey;;
2094
136
        psDefn->ProjParm[1] = dfNatOriginLong;
2095
136
        psDefn->ProjParmId[1] = ProjStraightVertPoleLongGeoKey;
2096
136
        psDefn->ProjParm[4] = dfNatOriginScale;
2097
136
        psDefn->ProjParmId[4] = ProjScaleAtNatOriginGeoKey;
2098
136
        psDefn->ProjParm[5] = dfFalseEasting;
2099
136
        psDefn->ProjParmId[5] = ProjFalseEastingGeoKey;
2100
136
        psDefn->ProjParm[6] = dfFalseNorthing;
2101
136
        psDefn->ProjParmId[6] = ProjFalseNorthingGeoKey;
2102
2103
136
        psDefn->nParms = 7;
2104
136
        break;
2105
2106
/* -------------------------------------------------------------------- */
2107
77
      case CT_LambertConfConic_2SP:
2108
/* -------------------------------------------------------------------- */
2109
77
        if( GTIFKeyGetDOUBLE(psGTIF, ProjStdParallel1GeoKey,
2110
77
                       &dfStdParallel1, 0, 1 ) == 0 )
2111
71
            dfStdParallel1 = 0.0;
2112
2113
77
        if( GTIFKeyGetDOUBLE(psGTIF, ProjStdParallel2GeoKey,
2114
77
                       &dfStdParallel2, 0, 1 ) == 0 )
2115
71
            dfStdParallel1 = 0.0;
2116
2117
77
        if( GTIFKeyGetDOUBLE(psGTIF, ProjNatOriginLongGeoKey,
2118
77
                       &dfNatOriginLong, 0, 1 ) == 0
2119
73
            && GTIFKeyGetDOUBLE(psGTIF, ProjFalseOriginLongGeoKey,
2120
73
                          &dfNatOriginLong, 0, 1 ) == 0
2121
68
            && GTIFKeyGetDOUBLE(psGTIF, ProjCenterLongGeoKey,
2122
68
                          &dfNatOriginLong, 0, 1 ) == 0 )
2123
60
            dfNatOriginLong = 0.0;
2124
2125
77
        if( GTIFKeyGetDOUBLE(psGTIF, ProjNatOriginLatGeoKey,
2126
77
                       &dfNatOriginLat, 0, 1 ) == 0
2127
73
            && GTIFKeyGetDOUBLE(psGTIF, ProjFalseOriginLatGeoKey,
2128
73
                          &dfNatOriginLat, 0, 1 ) == 0
2129
67
            && GTIFKeyGetDOUBLE(psGTIF, ProjCenterLatGeoKey,
2130
67
                          &dfNatOriginLat, 0, 1 ) == 0 )
2131
66
            dfNatOriginLat = 0.0;
2132
2133
        /* notdef: should transform to decimal degrees at this point */
2134
2135
77
        psDefn->ProjParm[0] = dfNatOriginLat;
2136
77
        psDefn->ProjParmId[0] = ProjFalseOriginLatGeoKey;
2137
77
        psDefn->ProjParm[1] = dfNatOriginLong;
2138
77
        psDefn->ProjParmId[1] = ProjFalseOriginLongGeoKey;
2139
77
        psDefn->ProjParm[2] = dfStdParallel1;
2140
77
        psDefn->ProjParmId[2] = ProjStdParallel1GeoKey;
2141
77
        psDefn->ProjParm[3] = dfStdParallel2;
2142
77
        psDefn->ProjParmId[3] = ProjStdParallel2GeoKey;
2143
77
        psDefn->ProjParm[5] = dfFalseEasting;
2144
77
        psDefn->ProjParmId[5] = ProjFalseEastingGeoKey;
2145
77
        psDefn->ProjParm[6] = dfFalseNorthing;
2146
77
        psDefn->ProjParmId[6] = ProjFalseNorthingGeoKey;
2147
2148
77
        psDefn->nParms = 7;
2149
77
        break;
2150
2151
/* -------------------------------------------------------------------- */
2152
75
      case CT_AlbersEqualArea:
2153
211
      case CT_EquidistantConic:
2154
/* -------------------------------------------------------------------- */
2155
211
        if( GTIFKeyGetDOUBLE(psGTIF, ProjStdParallel1GeoKey,
2156
211
                       &dfStdParallel1, 0, 1 ) == 0 )
2157
185
            dfStdParallel1 = 0.0;
2158
2159
211
        if( GTIFKeyGetDOUBLE(psGTIF, ProjStdParallel2GeoKey,
2160
211
                       &dfStdParallel2, 0, 1 ) == 0 )
2161
186
            dfStdParallel2 = 0.0;
2162
2163
211
        if( GTIFKeyGetDOUBLE(psGTIF, ProjNatOriginLongGeoKey,
2164
211
                       &dfNatOriginLong, 0, 1 ) == 0
2165
190
            && GTIFKeyGetDOUBLE(psGTIF, ProjFalseOriginLongGeoKey,
2166
190
                          &dfNatOriginLong, 0, 1 ) == 0
2167
189
            && GTIFKeyGetDOUBLE(psGTIF, ProjCenterLongGeoKey,
2168
189
                          &dfNatOriginLong, 0, 1 ) == 0 )
2169
152
            dfNatOriginLong = 0.0;
2170
2171
211
        if( GTIFKeyGetDOUBLE(psGTIF, ProjNatOriginLatGeoKey,
2172
211
                       &dfNatOriginLat, 0, 1 ) == 0
2173
154
            && GTIFKeyGetDOUBLE(psGTIF, ProjFalseOriginLatGeoKey,
2174
154
                          &dfNatOriginLat, 0, 1 ) == 0
2175
153
            && GTIFKeyGetDOUBLE(psGTIF, ProjCenterLatGeoKey,
2176
153
                          &dfNatOriginLat, 0, 1 ) == 0 )
2177
152
            dfNatOriginLat = 0.0;
2178
2179
        /* notdef: should transform to decimal degrees at this point */
2180
2181
211
        psDefn->ProjParm[0] = dfStdParallel1;
2182
211
        psDefn->ProjParmId[0] = ProjStdParallel1GeoKey;
2183
211
        psDefn->ProjParm[1] = dfStdParallel2;
2184
211
        psDefn->ProjParmId[1] = ProjStdParallel2GeoKey;
2185
211
        psDefn->ProjParm[2] = dfNatOriginLat;
2186
211
        psDefn->ProjParmId[2] = ProjNatOriginLatGeoKey;
2187
211
        psDefn->ProjParm[3] = dfNatOriginLong;
2188
211
        psDefn->ProjParmId[3] = ProjNatOriginLongGeoKey;
2189
211
        psDefn->ProjParm[5] = dfFalseEasting;
2190
211
        psDefn->ProjParmId[5] = ProjFalseEastingGeoKey;
2191
211
        psDefn->ProjParm[6] = dfFalseNorthing;
2192
211
        psDefn->ProjParmId[6] = ProjFalseNorthingGeoKey;
2193
2194
211
        psDefn->nParms = 7;
2195
211
        break;
2196
2197
/* -------------------------------------------------------------------- */
2198
91
      case CT_CylindricalEqualArea:
2199
/* -------------------------------------------------------------------- */
2200
91
        if( GTIFKeyGetDOUBLE(psGTIF, ProjStdParallel1GeoKey,
2201
91
                       &dfStdParallel1, 0, 1 ) == 0 )
2202
89
            dfStdParallel1 = 0.0;
2203
2204
91
        if( GTIFKeyGetDOUBLE(psGTIF, ProjNatOriginLongGeoKey,
2205
91
                       &dfNatOriginLong, 0, 1 ) == 0
2206
88
            && GTIFKeyGetDOUBLE(psGTIF, ProjFalseOriginLongGeoKey,
2207
88
                          &dfNatOriginLong, 0, 1 ) == 0
2208
87
            && GTIFKeyGetDOUBLE(psGTIF, ProjCenterLongGeoKey,
2209
87
                          &dfNatOriginLong, 0, 1 ) == 0 )
2210
87
            dfNatOriginLong = 0.0;
2211
2212
        /* notdef: should transform to decimal degrees at this point */
2213
2214
91
        psDefn->ProjParm[0] = dfStdParallel1;
2215
91
        psDefn->ProjParmId[0] = ProjStdParallel1GeoKey;
2216
91
        psDefn->ProjParm[1] = dfNatOriginLong;
2217
91
        psDefn->ProjParmId[1] = ProjNatOriginLongGeoKey;
2218
91
        psDefn->ProjParm[5] = dfFalseEasting;
2219
91
        psDefn->ProjParmId[5] = ProjFalseEastingGeoKey;
2220
91
        psDefn->ProjParm[6] = dfFalseNorthing;
2221
91
        psDefn->ProjParmId[6] = ProjFalseNorthingGeoKey;
2222
2223
91
        psDefn->nParms = 7;
2224
91
        break;
2225
10.9k
    }
2226
2227
85.0k
    for( int iParam = 0; iParam < psDefn->nParms; iParam++ )
2228
74.1k
    {
2229
74.1k
        switch( psDefn->ProjParmId[iParam] )
2230
74.1k
        {
2231
2232
/* -------------------------------------------------------------------- */
2233
/*      Normalize any linear parameters into meters.  In GeoTIFF        */
2234
/*      the linear projection parameter tags are normally in the        */
2235
/*      units of the coordinate system described.                       */
2236
/* -------------------------------------------------------------------- */
2237
10.2k
          case ProjFalseEastingGeoKey:
2238
20.5k
          case ProjFalseNorthingGeoKey:
2239
20.5k
          case ProjFalseOriginEastingGeoKey:
2240
20.5k
          case ProjFalseOriginNorthingGeoKey:
2241
20.5k
          case ProjCenterEastingGeoKey:
2242
20.5k
          case ProjCenterNorthingGeoKey:
2243
20.5k
            if( psDefn->UOMLengthInMeters != 0
2244
20.5k
                && psDefn->UOMLengthInMeters != 1.0 )
2245
3.01k
            {
2246
3.01k
                psDefn->ProjParm[iParam] *= psDefn->UOMLengthInMeters;
2247
3.01k
            }
2248
20.5k
            break;
2249
2250
/* -------------------------------------------------------------------- */
2251
/*      Normalize any angular parameters into degrees.  In GeoTIFF      */
2252
/*      the angular projection parameter tags are normally in the       */
2253
/*      units of GeogAngularUnit. Note: this conversion is only done    */
2254
/*      since libgeotiff 1.7.4                                          */
2255
/* -------------------------------------------------------------------- */
2256
2257
3.02k
          case ProjStdParallel1GeoKey:
2258
3.91k
          case ProjStdParallel2GeoKey:
2259
7.63k
          case ProjNatOriginLongGeoKey:
2260
11.4k
          case ProjNatOriginLatGeoKey:
2261
11.5k
          case ProjFalseOriginLongGeoKey:
2262
11.5k
          case ProjFalseOriginLatGeoKey:
2263
18.0k
          case ProjCenterLongGeoKey:
2264
24.3k
          case ProjCenterLatGeoKey:
2265
24.4k
          case ProjStraightVertPoleLongGeoKey:
2266
25.3k
          case ProjRectifiedGridAngleGeoKey:
2267
25.3k
            if( psDefn->UOMAngleInDegrees != 0
2268
25.3k
                && psDefn->UOMAngleInDegrees != 1.0 )
2269
13.2k
            {
2270
13.2k
                psDefn->ProjParm[iParam] *= psDefn->UOMAngleInDegrees;
2271
13.2k
            }
2272
25.3k
            break;
2273
2274
28.1k
          default:
2275
28.1k
            break;
2276
74.1k
        }
2277
74.1k
    }
2278
10.9k
}
2279
2280
/************************************************************************/
2281
/*                            GTIFGetDefn()                             */
2282
/************************************************************************/
2283
2284
/**
2285
@param psGTIF GeoTIFF information handle as returned by GTIFNew.
2286
@param psDefn Pointer to an existing GTIFDefn structure allocated by GTIFAllocDefn().
2287
2288
@return TRUE if the function has been successful, otherwise FALSE.
2289
2290
This function reads the coordinate system definition from a GeoTIFF file,
2291
and <i>normalizes</i> it into a set of component information using
2292
definitions from the EPSG database as provided by the PROJ library.
2293
This function is intended to simplify correct support for
2294
reading files with defined PCS (Projected Coordinate System) codes that
2295
wouldn't otherwise be directly known by application software by reducing
2296
it to the underlying projection method, parameters, datum, ellipsoid,
2297
prime meridian and units.<p>
2298
2299
The application should pass a pointer to an existing uninitialized
2300
GTIFDefn structure, and GTIFGetDefn() will fill it in.  The function
2301
currently always returns TRUE but in the future will return FALSE if
2302
the database is not found.  In any event, all geokeys actually found in the
2303
file will be copied into the GTIFDefn.  However, if the database isn't
2304
found, codes implied by other codes will not be set properly.<p>
2305
2306
GTIFGetDefn() will not generally work if the EPSG derived database cannot
2307
be found.<p>
2308
2309
The normalization methodology operates by fetching tags from the GeoTIFF
2310
file, and then setting all other tags implied by them in the structure.  The
2311
implied relationships are worked out by reading definitions from the
2312
various EPSG derived database tables.<p>
2313
2314
For instance, if a PCS (ProjectedCSTypeGeoKey) is found in the GeoTIFF file
2315
this code is used to lookup a record in the database.
2316
For example given the PCS 26746 we can find the name
2317
(NAD27 / California zone VI), the GCS 4257 (NAD27), and the ProjectionCode
2318
10406 (California CS27 zone VI).  The GCS, and ProjectionCode can in turn
2319
be looked up in other tables until all the details of units, ellipsoid,
2320
prime meridian, datum, projection (LambertConfConic_2SP) and projection
2321
parameters are established.  A full listgeo dump of a file
2322
for this result might look like the following, all based on a single PCS
2323
value:<p>
2324
2325
<pre>
2326
% listgeo -norm ~/data/geotiff/pci_eg/spaf27.tif
2327
Geotiff_Information:
2328
   Version: 1
2329
   Key_Revision: 1.0
2330
   Tagged_Information:
2331
      ModelTiepointTag (2,3):
2332
         0                0                0
2333
         1577139.71       634349.176       0
2334
      ModelPixelScaleTag (1,3):
2335
         195.509321       198.32184        0
2336
      End_Of_Tags.
2337
   Keyed_Information:
2338
      GTModelTypeGeoKey (Short,1): ModelTypeProjected
2339
      GTRasterTypeGeoKey (Short,1): RasterPixelIsArea
2340
      ProjectedCSTypeGeoKey (Short,1): PCS_NAD27_California_VI
2341
      End_Of_Keys.
2342
   End_Of_Geotiff.
2343
2344
PCS = 26746 (NAD27 / California zone VI)
2345
Projection = 10406 (California CS27 zone VI)
2346
Projection Method: CT_LambertConfConic_2SP
2347
   ProjStdParallel1GeoKey: 33.883333
2348
   ProjStdParallel2GeoKey: 32.766667
2349
   ProjFalseOriginLatGeoKey: 32.166667
2350
   ProjFalseOriginLongGeoKey: -116.233333
2351
   ProjFalseEastingGeoKey: 609601.219202
2352
   ProjFalseNorthingGeoKey: 0.000000
2353
GCS: 4267/NAD27
2354
Datum: 6267/North American Datum 1927
2355
Ellipsoid: 7008/Clarke 1866 (6378206.40,6356583.80)
2356
Prime Meridian: 8901/Greenwich (0.000000)
2357
Projection Linear Units: 9003/US survey foot (0.304801m)
2358
</pre>
2359
2360
Note that GTIFGetDefn() does not inspect or return the tiepoints and scale.
2361
This must be handled separately as it normally would.  It is intended to
2362
simplify capture and normalization of the coordinate system definition.
2363
Note that GTIFGetDefn() also does the following things:
2364
2365
<ol>
2366
<li> Convert all angular values to decimal degrees.
2367
<li> Convert all linear values to meters.
2368
<li> Return the linear units and conversion to meters for the tiepoints and
2369
scale (though the tiepoints and scale remain in their native units).
2370
<li> When reading projection parameters a variety of differences between
2371
different GeoTIFF generators are handled, and a normalized set of parameters
2372
for each projection are always returned.
2373
</ol>
2374
2375
Code fields in the GTIFDefn are filled with KvUserDefined if there is not
2376
value to assign.  The parameter lists for each of the underlying projection
2377
transform methods can be found at the
2378
<a href="http://www.remotesensing.org/geotiff/proj_list">Projections</a>
2379
page.  Note that nParms will be set based on the maximum parameter used.
2380
Some of the parameters may not be used in which case the
2381
GTIFDefn::ProjParmId[] will
2382
be zero.  This is done to retain correspondence to the EPSG parameter
2383
numbering scheme.<p>
2384
2385
The
2386
<a href="http://www.remotesensing.org/cgi-bin/cvsweb.cgi/~checkout~/osrs/geotiff/libgeotiff/geotiff_proj4.c">geotiff_proj4.c</a> module distributed with libgeotiff can
2387
be used as an example of code that converts a GTIFDefn into another projection
2388
system.<p>
2389
2390
@see GTIFKeySet()
2391
2392
*/
2393
2394
int GTIFGetDefn( GTIF * psGTIF, GTIFDefn * psDefn )
2395
2396
125k
{
2397
125k
    if( !GTIFGetPROJContext(psGTIF, TRUE, NULL) )
2398
0
    {
2399
0
        return FALSE;
2400
0
    }
2401
2402
/* -------------------------------------------------------------------- */
2403
/*      Initially we default all the information we can.                */
2404
/* -------------------------------------------------------------------- */
2405
125k
    psDefn->DefnSet = 1;
2406
125k
    psDefn->Model = KvUserDefined;
2407
125k
    psDefn->PCS = KvUserDefined;
2408
125k
    psDefn->GCS = KvUserDefined;
2409
125k
    psDefn->UOMLength = KvUserDefined;
2410
125k
    psDefn->UOMLengthInMeters = 1.0;
2411
125k
    psDefn->UOMAngle = KvUserDefined;
2412
125k
    psDefn->UOMAngleInDegrees = 1.0;
2413
125k
    psDefn->Datum = KvUserDefined;
2414
125k
    psDefn->Ellipsoid = KvUserDefined;
2415
125k
    psDefn->SemiMajor = 0.0;
2416
125k
    psDefn->SemiMinor = 0.0;
2417
125k
    psDefn->PM = KvUserDefined;
2418
125k
    psDefn->PMLongToGreenwich = 0.0;
2419
125k
#if !defined(GEO_NORMALIZE_DISABLE_TOWGS84)
2420
125k
    psDefn->TOWGS84Count = 0;
2421
125k
    memset( psDefn->TOWGS84, 0, sizeof(psDefn->TOWGS84) );
2422
125k
#endif
2423
2424
125k
    psDefn->ProjCode = KvUserDefined;
2425
125k
    psDefn->Projection = KvUserDefined;
2426
125k
    psDefn->CTProjection = KvUserDefined;
2427
2428
125k
    psDefn->nParms = 0;
2429
1.37M
    for( int i = 0; i < MAX_GTIF_PROJPARMS; i++ )
2430
1.25M
    {
2431
1.25M
        psDefn->ProjParm[i] = 0.0;
2432
1.25M
        psDefn->ProjParmId[i] = 0;
2433
1.25M
    }
2434
2435
125k
    psDefn->MapSys = KvUserDefined;
2436
125k
    psDefn->Zone = 0;
2437
2438
/* -------------------------------------------------------------------- */
2439
/*      Do we have any geokeys?                                         */
2440
/* -------------------------------------------------------------------- */
2441
125k
    {
2442
125k
        int     nKeyCount = 0;
2443
125k
        int     anVersion[3];
2444
125k
        GTIFDirectoryInfo( psGTIF, anVersion, &nKeyCount );
2445
2446
125k
        if( nKeyCount == 0 )
2447
69.4k
        {
2448
69.4k
            psDefn->DefnSet = 0;
2449
69.4k
            return FALSE;
2450
69.4k
        }
2451
125k
    }
2452
2453
/* -------------------------------------------------------------------- */
2454
/*  Try to get the overall model type.        */
2455
/* -------------------------------------------------------------------- */
2456
55.8k
    GTIFKeyGetSSHORT(psGTIF,GTModelTypeGeoKey,&(psDefn->Model));
2457
2458
/* -------------------------------------------------------------------- */
2459
/*  Extract the Geog units.           */
2460
/* -------------------------------------------------------------------- */
2461
55.8k
    short nGeogUOMLinear = 9001; /* Linear_Meter */
2462
55.8k
    if( GTIFKeyGetSSHORT(psGTIF, GeogLinearUnitsGeoKey, &nGeogUOMLinear) == 1 )
2463
297
    {
2464
297
        psDefn->UOMLength = nGeogUOMLinear;
2465
297
    }
2466
2467
/* -------------------------------------------------------------------- */
2468
/*      Try to get a PCS.                                               */
2469
/* -------------------------------------------------------------------- */
2470
55.8k
    if( GTIFKeyGetSSHORT(psGTIF,ProjectedCSTypeGeoKey, &(psDefn->PCS)) == 1
2471
36.7k
        && psDefn->PCS != KvUserDefined )
2472
35.7k
    {
2473
        /*
2474
         * Translate this into useful information.
2475
         */
2476
35.7k
        GTIFGetPCSInfoEx( psGTIF->pj_context,
2477
35.7k
                          psDefn->PCS, NULL, &(psDefn->ProjCode),
2478
35.7k
                          &(psDefn->UOMLength), &(psDefn->GCS) );
2479
35.7k
    }
2480
2481
/* -------------------------------------------------------------------- */
2482
/*       If we have the PCS code, but didn't find it in the database    */
2483
/*      (likely because we can't find them) we will try some ``jiffy    */
2484
/*      rules'' for UTM and state plane.                                */
2485
/* -------------------------------------------------------------------- */
2486
55.8k
    if( psDefn->PCS != KvUserDefined && psDefn->ProjCode == KvUserDefined )
2487
19.3k
    {
2488
19.3k
        int nZone;
2489
19.3k
        int nGCS = psDefn->GCS;
2490
2491
19.3k
        const int nMapSys = GTIFPCSToMapSys( psDefn->PCS, &nGCS, &nZone );
2492
19.3k
        if( nMapSys != KvUserDefined )
2493
4.72k
        {
2494
4.72k
            psDefn->ProjCode = (short) GTIFMapSysToProj( nMapSys, nZone );
2495
4.72k
            psDefn->GCS = (short) nGCS;
2496
4.72k
        }
2497
19.3k
    }
2498
2499
/* -------------------------------------------------------------------- */
2500
/*      If the Proj_ code is specified directly, use that.              */
2501
/* -------------------------------------------------------------------- */
2502
55.8k
    if( psDefn->ProjCode == KvUserDefined )
2503
34.6k
        GTIFKeyGetSSHORT(psGTIF, ProjectionGeoKey, &(psDefn->ProjCode));
2504
2505
55.8k
    if( psDefn->ProjCode != KvUserDefined )
2506
23.7k
    {
2507
        /*
2508
         * We have an underlying projection transformation value.  Look
2509
         * this up.  For a PCS of ``WGS 84 / UTM 11'' the transformation
2510
         * would be Transverse Mercator, with a particular set of options.
2511
         * The nProjTRFCode itself would correspond to the name
2512
         * ``UTM zone 11N'', and doesn't include datum info.
2513
         */
2514
23.7k
        GTIFGetProjTRFInfoEx( psGTIF->pj_context,
2515
23.7k
                              psDefn->ProjCode, NULL, &(psDefn->Projection),
2516
23.7k
                              psDefn->ProjParm );
2517
2518
        /*
2519
         * Set the GeoTIFF identity of the parameters.
2520
         */
2521
23.7k
        psDefn->CTProjection = (short)
2522
23.7k
            EPSGProjMethodToCTProjMethod( psDefn->Projection, FALSE );
2523
2524
23.7k
        SetGTParamIds( EPSGProjMethodToCTProjMethod(psDefn->Projection, TRUE),
2525
23.7k
                      psDefn->Projection,
2526
23.7k
                      psDefn->ProjParmId, NULL);
2527
23.7k
        psDefn->nParms = 7;
2528
23.7k
    }
2529
2530
/* -------------------------------------------------------------------- */
2531
/*      Try to get a GCS.  If found, it will override any implied by    */
2532
/*      the PCS.                                                        */
2533
/* -------------------------------------------------------------------- */
2534
55.8k
    GTIFKeyGetSSHORT(psGTIF, GeographicTypeGeoKey, &(psDefn->GCS));
2535
55.8k
    if( psDefn->GCS < 1 || psDefn->GCS >= KvUserDefined )
2536
32.5k
        psDefn->GCS = KvUserDefined;
2537
2538
/* -------------------------------------------------------------------- */
2539
/*      Derive the datum, and prime meridian from the GCS.              */
2540
/* -------------------------------------------------------------------- */
2541
55.8k
    if( psDefn->GCS != KvUserDefined )
2542
23.2k
    {
2543
23.2k
        GTIFGetGCSInfoEx( psGTIF->pj_context,
2544
23.2k
                          psDefn->GCS, NULL, &(psDefn->Datum), &(psDefn->PM),
2545
23.2k
                          &(psDefn->UOMAngle) );
2546
23.2k
    }
2547
2548
/* -------------------------------------------------------------------- */
2549
/*      Handle the GCS angular units.  GeogAngularUnitsGeoKey           */
2550
/*      overrides the GCS or PCS setting.                               */
2551
/* -------------------------------------------------------------------- */
2552
55.8k
    GTIFKeyGetSSHORT(psGTIF, GeogAngularUnitsGeoKey, &(psDefn->UOMAngle));
2553
55.8k
    if( psDefn->UOMAngle != KvUserDefined )
2554
46.8k
    {
2555
46.8k
        GTIFGetUOMAngleInfoEx( psGTIF->pj_context,
2556
46.8k
                               psDefn->UOMAngle, NULL,
2557
46.8k
                               &(psDefn->UOMAngleInDegrees) );
2558
46.8k
    }
2559
2560
/* -------------------------------------------------------------------- */
2561
/*      Check for a datum setting, and then use the datum to derive     */
2562
/*      an ellipsoid.                                                   */
2563
/* -------------------------------------------------------------------- */
2564
55.8k
    GTIFKeyGetSSHORT(psGTIF, GeogGeodeticDatumGeoKey, &(psDefn->Datum));
2565
2566
55.8k
    if( psDefn->Datum != KvUserDefined )
2567
22.8k
    {
2568
22.8k
        GTIFGetDatumInfoEx( psGTIF->pj_context,
2569
22.8k
                            psDefn->Datum, NULL, &(psDefn->Ellipsoid) );
2570
22.8k
    }
2571
2572
/* -------------------------------------------------------------------- */
2573
/*      Check for an explicit ellipsoid.  Use the ellipsoid to          */
2574
/*      derive the ellipsoid characteristics, if possible.              */
2575
/* -------------------------------------------------------------------- */
2576
55.8k
    GTIFKeyGetSSHORT(psGTIF, GeogEllipsoidGeoKey, &(psDefn->Ellipsoid));
2577
2578
55.8k
    if( psDefn->Ellipsoid != KvUserDefined )
2579
19.5k
    {
2580
19.5k
        GTIFGetEllipsoidInfoEx( psGTIF->pj_context,
2581
19.5k
                                psDefn->Ellipsoid, NULL,
2582
19.5k
                                &(psDefn->SemiMajor), &(psDefn->SemiMinor) );
2583
19.5k
    }
2584
2585
/* -------------------------------------------------------------------- */
2586
/*      Check for overridden ellipsoid parameters.  It would be nice    */
2587
/*      to warn if they conflict with provided information, but for     */
2588
/*      now we just override.                                           */
2589
/* -------------------------------------------------------------------- */
2590
55.8k
    CPL_IGNORE_RET_VAL_INT(GTIFKeyGetDOUBLE(psGTIF, GeogSemiMajorAxisGeoKey, &(psDefn->SemiMajor), 0, 1 ));
2591
55.8k
    CPL_IGNORE_RET_VAL_INT(GTIFKeyGetDOUBLE(psGTIF, GeogSemiMinorAxisGeoKey, &(psDefn->SemiMinor), 0, 1 ));
2592
2593
55.8k
    double dfInvFlattening;
2594
55.8k
    if( GTIFKeyGetDOUBLE(psGTIF, GeogInvFlatteningGeoKey, &dfInvFlattening,
2595
55.8k
                   0, 1 ) == 1 )
2596
1.07k
    {
2597
1.07k
        if( dfInvFlattening != 0.0 )
2598
1.04k
            psDefn->SemiMinor =
2599
1.04k
                psDefn->SemiMajor * (1 - 1.0/dfInvFlattening);
2600
30
        else
2601
30
            psDefn->SemiMinor = psDefn->SemiMajor;
2602
1.07k
    }
2603
2604
/* -------------------------------------------------------------------- */
2605
/*      Get the prime meridian info.                                    */
2606
/* -------------------------------------------------------------------- */
2607
55.8k
    GTIFKeyGetSSHORT(psGTIF, GeogPrimeMeridianGeoKey, &(psDefn->PM));
2608
2609
55.8k
    if( psDefn->PM != KvUserDefined )
2610
19.9k
    {
2611
19.9k
        GTIFGetPMInfoEx( psGTIF->pj_context,
2612
19.9k
                         psDefn->PM, NULL, &(psDefn->PMLongToGreenwich) );
2613
19.9k
    }
2614
35.8k
    else
2615
35.8k
    {
2616
35.8k
        CPL_IGNORE_RET_VAL_INT(GTIFKeyGetDOUBLE(psGTIF, GeogPrimeMeridianLongGeoKey,
2617
35.8k
                   &(psDefn->PMLongToGreenwich), 0, 1 ));
2618
2619
35.8k
        psDefn->PMLongToGreenwich =
2620
35.8k
            GTIFAngleToDD( psDefn->PMLongToGreenwich,
2621
35.8k
                           psDefn->UOMAngle );
2622
35.8k
    }
2623
2624
/* -------------------------------------------------------------------- */
2625
/*      Get the TOWGS84 parameters.                                     */
2626
/* -------------------------------------------------------------------- */
2627
55.8k
#if !defined(GEO_NORMALIZE_DISABLE_TOWGS84)
2628
55.8k
    psDefn->TOWGS84Count =
2629
55.8k
        (short)GTIFKeyGetDOUBLE(psGTIF, GeogTOWGS84GeoKey, psDefn->TOWGS84, 0, 7 );
2630
55.8k
#endif
2631
2632
/* -------------------------------------------------------------------- */
2633
/*      Have the projection units of measure been overridden?  We       */
2634
/*      should likely be doing something about angular units too,       */
2635
/*      but these are very rarely not decimal degrees for actual        */
2636
/*      file coordinates.                                               */
2637
/* -------------------------------------------------------------------- */
2638
55.8k
    GTIFKeyGetSSHORT(psGTIF,ProjLinearUnitsGeoKey,&(psDefn->UOMLength));
2639
2640
55.8k
    if( psDefn->UOMLength != KvUserDefined )
2641
36.7k
    {
2642
36.7k
        GTIFGetUOMLengthInfoEx( psGTIF->pj_context,
2643
36.7k
                                psDefn->UOMLength, NULL,
2644
36.7k
                                &(psDefn->UOMLengthInMeters) );
2645
36.7k
    }
2646
19.0k
    else
2647
19.0k
    {
2648
19.0k
        CPL_IGNORE_RET_VAL_INT(GTIFKeyGetDOUBLE(psGTIF,ProjLinearUnitSizeGeoKey,&(psDefn->UOMLengthInMeters),0,1));
2649
19.0k
    }
2650
2651
/* -------------------------------------------------------------------- */
2652
/*      Handle a variety of user defined transform types.               */
2653
/* -------------------------------------------------------------------- */
2654
55.8k
    if( GTIFKeyGetSSHORT(psGTIF,ProjCoordTransGeoKey,
2655
55.8k
                   &(psDefn->CTProjection)) == 1)
2656
10.9k
    {
2657
10.9k
        GTIFFetchProjParms( psGTIF, psDefn );
2658
10.9k
    }
2659
2660
/* -------------------------------------------------------------------- */
2661
/*      Try to set the zoned map system information.                    */
2662
/* -------------------------------------------------------------------- */
2663
55.8k
    psDefn->MapSys = GTIFProjToMapSys( psDefn->ProjCode, &(psDefn->Zone) );
2664
2665
/* -------------------------------------------------------------------- */
2666
/*      If this is UTM, and we were unable to extract the projection    */
2667
/*      parameters from the database just set them directly now,        */
2668
/*      since it's pretty easy, and a common case.                      */
2669
/* -------------------------------------------------------------------- */
2670
55.8k
    if( (psDefn->MapSys == MapSys_UTM_North
2671
50.5k
         || psDefn->MapSys == MapSys_UTM_South)
2672
8.82k
        && psDefn->CTProjection == KvUserDefined )
2673
9
    {
2674
9
        psDefn->CTProjection = CT_TransverseMercator;
2675
9
        psDefn->nParms = 7;
2676
9
        psDefn->ProjParmId[0] = ProjNatOriginLatGeoKey;
2677
9
        psDefn->ProjParm[0] = 0.0;
2678
2679
9
        psDefn->ProjParmId[1] = ProjNatOriginLongGeoKey;
2680
9
        psDefn->ProjParm[1] = psDefn->Zone*6 - 183.0;
2681
2682
9
        psDefn->ProjParmId[4] = ProjScaleAtNatOriginGeoKey;
2683
9
        psDefn->ProjParm[4] = 0.9996;
2684
2685
9
        psDefn->ProjParmId[5] = ProjFalseEastingGeoKey;
2686
9
        psDefn->ProjParm[5] = 500000.0;
2687
2688
9
        psDefn->ProjParmId[6] = ProjFalseNorthingGeoKey;
2689
2690
9
        if( psDefn->MapSys == MapSys_UTM_North )
2691
1
            psDefn->ProjParm[6] = 0.0;
2692
8
        else
2693
8
            psDefn->ProjParm[6] = 10000000.0;
2694
9
    }
2695
2696
55.8k
    return TRUE;
2697
125k
}
2698
2699
/************************************************************************/
2700
/*                            GTIFDecToDMS()                            */
2701
/*                                                                      */
2702
/*      Convenient function to translate decimal degrees to DMS         */
2703
/*      format for reporting to a user.                                 */
2704
/************************************************************************/
2705
2706
const char *GTIFDecToDMS( double dfAngle, const char * pszAxis,
2707
                          int nPrecision )
2708
0
{
2709
0
    if( !(dfAngle >= -360 && dfAngle <= 360) )
2710
0
        return "";
2711
2712
0
    double dfRound = 0.5/60;
2713
0
    for( int i = 0; i < nPrecision; i++ )
2714
0
        dfRound = dfRound * 0.1;
2715
2716
0
    int nDegrees = (int) ABS(dfAngle);
2717
0
    int nMinutes = (int) ((ABS(dfAngle) - nDegrees) * 60 + dfRound);
2718
0
    if( nMinutes == 60 )
2719
0
    {
2720
0
        nDegrees ++;
2721
0
        nMinutes = 0;
2722
0
    }
2723
0
    const double dfSeconds = ABS((ABS(dfAngle) * 3600 - nDegrees*3600 - nMinutes*60));
2724
2725
0
    const char  *pszHemisphere = NULL;
2726
0
    if( EQUAL(pszAxis,"Long") && dfAngle < 0.0 )
2727
0
        pszHemisphere = "W";
2728
0
    else if( EQUAL(pszAxis,"Long") )
2729
0
        pszHemisphere = "E";
2730
0
    else if( dfAngle < 0.0 )
2731
0
        pszHemisphere = "S";
2732
0
    else
2733
0
        pszHemisphere = "N";
2734
2735
0
    char szFormat[30];
2736
0
    snprintf(szFormat, sizeof(szFormat), "%%3dd%%2d\'%%%d.%df\"%s",
2737
0
             nPrecision+3, nPrecision, pszHemisphere);
2738
0
    static char szBuffer[50];
2739
0
    snprintf(
2740
0
        szBuffer, sizeof(szBuffer), szFormat, nDegrees, nMinutes, dfSeconds);
2741
2742
0
    return szBuffer;
2743
0
}
2744
2745
/************************************************************************/
2746
/*                           GTIFPrintDefn()                            */
2747
/*                                                                      */
2748
/*      Report the contents of a GTIFDefn structure ... mostly for      */
2749
/*      debugging.                                                      */
2750
/************************************************************************/
2751
2752
void GTIFPrintDefnEx( GTIF *psGTIF, GTIFDefn * psDefn, FILE * fp )
2753
2754
0
{
2755
0
    GTIFGetPROJContext(psGTIF, TRUE, NULL);
2756
2757
/* -------------------------------------------------------------------- */
2758
/*      Do we have anything to report?                                  */
2759
/* -------------------------------------------------------------------- */
2760
0
    if( !psDefn->DefnSet )
2761
0
    {
2762
0
        fprintf( fp, "No GeoKeys found.\n" );
2763
0
        return;
2764
0
    }
2765
2766
/* -------------------------------------------------------------------- */
2767
/*      Get the PCS name if possible.                                   */
2768
/* -------------------------------------------------------------------- */
2769
0
    if( psDefn->PCS != KvUserDefined )
2770
0
    {
2771
0
        char  *pszPCSName = NULL;
2772
2773
0
        if( psGTIF->pj_context )
2774
0
        {
2775
0
            GTIFGetPCSInfoEx( psGTIF->pj_context,
2776
0
                              psDefn->PCS, &pszPCSName, NULL, NULL, NULL );
2777
0
        }
2778
0
        if( pszPCSName == NULL )
2779
0
            pszPCSName = CPLStrdup("name unknown");
2780
2781
0
        fprintf( fp, "PCS = %d (%s)\n", psDefn->PCS, pszPCSName );
2782
0
        CPLFree( pszPCSName );
2783
0
    }
2784
2785
/* -------------------------------------------------------------------- */
2786
/*  Dump the projection code if possible.       */
2787
/* -------------------------------------------------------------------- */
2788
0
    if( psDefn->ProjCode != KvUserDefined )
2789
0
    {
2790
0
        char  *pszTRFName = NULL;
2791
2792
0
        if( psGTIF->pj_context )
2793
0
        {
2794
0
            GTIFGetProjTRFInfoEx( psGTIF->pj_context,
2795
0
                                  psDefn->ProjCode, &pszTRFName, NULL, NULL );
2796
0
        }
2797
0
        if( pszTRFName == NULL )
2798
0
            pszTRFName = CPLStrdup("");
2799
2800
0
        fprintf( fp, "Projection = %d (%s)\n",
2801
0
                 psDefn->ProjCode, pszTRFName );
2802
2803
0
        CPLFree( pszTRFName );
2804
0
    }
2805
2806
/* -------------------------------------------------------------------- */
2807
/*      Try to dump the projection method name, and parameters if possible.*/
2808
/* -------------------------------------------------------------------- */
2809
0
    if( psDefn->CTProjection != KvUserDefined )
2810
0
    {
2811
0
        const char *pszProjectionMethodName =
2812
0
            GTIFValueNameEx(psGTIF,
2813
0
                            ProjCoordTransGeoKey,
2814
0
                            psDefn->CTProjection);
2815
2816
0
        if( pszProjectionMethodName == NULL )
2817
0
            pszProjectionMethodName = "(unknown)";
2818
2819
0
        fprintf( fp, "Projection Method: %s\n", pszProjectionMethodName );
2820
2821
0
        for( int i = 0; i < psDefn->nParms; i++ )
2822
0
        {
2823
0
            if( psDefn->ProjParmId[i] == 0 )
2824
0
                continue;
2825
2826
0
            char* pszName = GTIFKeyName((geokey_t) psDefn->ProjParmId[i]);
2827
0
            if( pszName == NULL )
2828
0
                pszName = "(unknown)";
2829
2830
0
            if( i < 4 )
2831
0
            {
2832
0
                char  *pszAxisName;
2833
2834
0
                if( strstr(pszName,"Long") != NULL )
2835
0
                    pszAxisName = "Long";
2836
0
                else if( strstr(pszName,"Lat") != NULL )
2837
0
                    pszAxisName = "Lat";
2838
0
                else
2839
0
                    pszAxisName = "?";
2840
2841
0
                fprintf( fp, "   %s: %f (%s)\n",
2842
0
                         pszName, psDefn->ProjParm[i],
2843
0
                         GTIFDecToDMS( psDefn->ProjParm[i], pszAxisName, 2 ) );
2844
0
            }
2845
0
            else if( i == 4 )
2846
0
                fprintf( fp, "   %s: %f\n", pszName, psDefn->ProjParm[i] );
2847
0
            else
2848
0
                fprintf( fp, "   %s: %f m\n", pszName, psDefn->ProjParm[i] );
2849
0
        }
2850
0
    }
2851
2852
/* -------------------------------------------------------------------- */
2853
/*      Report the GCS name, and number.                                */
2854
/* -------------------------------------------------------------------- */
2855
0
    if( psDefn->GCS != KvUserDefined )
2856
0
    {
2857
0
        char  *pszName = NULL;
2858
2859
0
        if( psGTIF->pj_context )
2860
0
        {
2861
0
            GTIFGetGCSInfoEx( psGTIF->pj_context,
2862
0
                              psDefn->GCS, &pszName, NULL, NULL, NULL );
2863
0
        }
2864
0
        if( pszName == NULL )
2865
0
            pszName = CPLStrdup("(unknown)");
2866
2867
0
        fprintf( fp, "GCS: %d/%s\n", psDefn->GCS, pszName );
2868
0
        CPLFree( pszName );
2869
0
    }
2870
2871
/* -------------------------------------------------------------------- */
2872
/*      Report the datum name.                                          */
2873
/* -------------------------------------------------------------------- */
2874
0
    if( psDefn->Datum != KvUserDefined )
2875
0
    {
2876
0
        char  *pszName = NULL;
2877
2878
0
        if( psGTIF->pj_context )
2879
0
        {
2880
0
            GTIFGetDatumInfoEx( psGTIF->pj_context,
2881
0
                                psDefn->Datum, &pszName, NULL );
2882
0
        }
2883
0
        if( pszName == NULL )
2884
0
            pszName = CPLStrdup("(unknown)");
2885
2886
0
        fprintf( fp, "Datum: %d/%s\n", psDefn->Datum, pszName );
2887
0
        CPLFree( pszName );
2888
0
    }
2889
2890
/* -------------------------------------------------------------------- */
2891
/*      Report the ellipsoid.                                           */
2892
/* -------------------------------------------------------------------- */
2893
0
    if( psDefn->Ellipsoid != KvUserDefined )
2894
0
    {
2895
0
        char  *pszName = NULL;
2896
2897
0
        if( psGTIF->pj_context )
2898
0
        {
2899
0
            GTIFGetEllipsoidInfoEx( psGTIF->pj_context,
2900
0
                                    psDefn->Ellipsoid, &pszName, NULL, NULL );
2901
0
        }
2902
0
        if( pszName == NULL )
2903
0
            pszName = CPLStrdup("(unknown)");
2904
2905
0
        fprintf( fp, "Ellipsoid: %d/%s (%.2f,%.2f)\n",
2906
0
                 psDefn->Ellipsoid, pszName,
2907
0
                 psDefn->SemiMajor, psDefn->SemiMinor );
2908
0
        CPLFree( pszName );
2909
0
    }
2910
2911
/* -------------------------------------------------------------------- */
2912
/*      Report the prime meridian.                                      */
2913
/* -------------------------------------------------------------------- */
2914
0
    if( psDefn->PM != KvUserDefined )
2915
0
    {
2916
0
        char  *pszName = NULL;
2917
2918
0
        if( psGTIF->pj_context )
2919
0
        {
2920
0
            GTIFGetPMInfoEx( psGTIF->pj_context,
2921
0
                             psDefn->PM, &pszName, NULL );
2922
0
        }
2923
2924
0
        if( pszName == NULL )
2925
0
            pszName = CPLStrdup("(unknown)");
2926
2927
0
        fprintf( fp, "Prime Meridian: %d/%s (%f/%s)\n",
2928
0
                 psDefn->PM, pszName,
2929
0
                 psDefn->PMLongToGreenwich,
2930
0
                 GTIFDecToDMS( psDefn->PMLongToGreenwich, "Long", 2 ) );
2931
0
        CPLFree( pszName );
2932
0
    }
2933
2934
/* -------------------------------------------------------------------- */
2935
/*      Report TOWGS84 parameters.                                      */
2936
/* -------------------------------------------------------------------- */
2937
0
#if !defined(GEO_NORMALIZE_DISABLE_TOWGS84)
2938
0
    if( psDefn->TOWGS84Count > 0 )
2939
0
    {
2940
0
        fprintf( fp, "TOWGS84: " );
2941
2942
0
        for( int i = 0; i < psDefn->TOWGS84Count; i++ )
2943
0
        {
2944
0
            if( i > 0 )
2945
0
                fprintf( fp, "," );
2946
0
            fprintf( fp, "%g", psDefn->TOWGS84[i] );
2947
0
        }
2948
2949
0
        fprintf( fp, "\n" );
2950
0
    }
2951
0
#endif
2952
2953
/* -------------------------------------------------------------------- */
2954
/*      Report the projection units of measure (currently just          */
2955
/*      linear).                                                        */
2956
/* -------------------------------------------------------------------- */
2957
0
    if( psDefn->UOMLength != KvUserDefined )
2958
0
    {
2959
0
        char  *pszName = NULL;
2960
2961
0
        if( psGTIF->pj_context )
2962
0
        {
2963
0
            GTIFGetUOMLengthInfoEx(
2964
0
                psGTIF->pj_context, psDefn->UOMLength, &pszName, NULL );
2965
0
        }
2966
0
        if( pszName == NULL )
2967
0
            pszName = CPLStrdup( "(unknown)" );
2968
2969
0
        fprintf( fp, "Projection Linear Units: %d/%s (%fm)\n",
2970
0
                 psDefn->UOMLength, pszName, psDefn->UOMLengthInMeters );
2971
0
        CPLFree( pszName );
2972
0
    }
2973
0
    else
2974
0
    {
2975
0
        fprintf( fp, "Projection Linear Units: User-Defined (%fm)\n",
2976
0
                 psDefn->UOMLengthInMeters );
2977
0
    }
2978
0
}
2979
2980
void GTIFPrintDefn( GTIFDefn * psDefn, FILE * fp )
2981
0
{
2982
0
    GTIF *psGTIF = GTIFNew(NULL);
2983
0
    if( psGTIF )
2984
0
    {
2985
0
        GTIFPrintDefnEx(psGTIF, psDefn, fp);
2986
0
        GTIFFree(psGTIF);
2987
0
    }
2988
0
}
2989
2990
/************************************************************************/
2991
/*                           GTIFFreeMemory()                           */
2992
/*                                                                      */
2993
/*      Externally visible function to free memory allocated within     */
2994
/*      geo_normalize.c.                                                */
2995
/************************************************************************/
2996
2997
void GTIFFreeMemory( char * pMemory )
2998
2999
171k
{
3000
171k
    if( pMemory != NULL )
3001
151k
        VSIFree( pMemory );
3002
171k
}
3003
3004
/************************************************************************/
3005
/*                           GTIFAllocDefn()                            */
3006
/*                                                                      */
3007
/*      This allocates a GTIF structure in such a way that the          */
3008
/*      calling application doesn't need to know the size and           */
3009
/*      initializes it appropriately.                                   */
3010
/************************************************************************/
3011
3012
GTIFDefn *GTIFAllocDefn()
3013
125k
{
3014
125k
    return (GTIFDefn *) CPLCalloc(sizeof(GTIFDefn),1);
3015
125k
}
3016
3017
/************************************************************************/
3018
/*                            GTIFFreeDefn()                            */
3019
/*                                                                      */
3020
/*      Free a GTIF structure allocated by GTIFAllocDefn().             */
3021
/************************************************************************/
3022
3023
void GTIFFreeDefn( GTIFDefn *defn )
3024
125k
{
3025
125k
    VSIFree( defn );
3026
125k
}
3027
3028
/************************************************************************/
3029
/*                       GTIFAttachPROJContext()                        */
3030
/*                                                                      */
3031
/*      Attach an existing PROJ context to the GTIF handle, but         */
3032
/*      ownership of the context remains to the caller.                 */
3033
/************************************************************************/
3034
3035
void GTIFAttachPROJContext( GTIF *psGTIF, void* pjContext )
3036
263k
{
3037
263k
    if( psGTIF->own_pj_context )
3038
0
    {
3039
0
        proj_context_destroy(psGTIF->pj_context);
3040
0
    }
3041
263k
    psGTIF->own_pj_context = FALSE;
3042
263k
    psGTIF->pj_context = (PJ_CONTEXT*) pjContext;
3043
263k
}
3044
3045
/************************************************************************/
3046
/*                         GTIFGetPROJContext()                         */
3047
/*                                                                      */
3048
/*      Return the PROJ context attached to the GTIF handle.            */
3049
/*      If it has not yet been instantiated and instantiateIfNeeded=TRUE*/
3050
/*      then, it will be instantiated (and owned by GTIF handle).       */
3051
/************************************************************************/
3052
3053
void *GTIFGetPROJContext( GTIF *psGTIF, int instantiateIfNeeded,
3054
                          int* out_gtif_own_pj_context )
3055
225k
{
3056
225k
    if( psGTIF->pj_context || !instantiateIfNeeded )
3057
225k
    {
3058
225k
        if( out_gtif_own_pj_context )
3059
0
        {
3060
0
            *out_gtif_own_pj_context = psGTIF->own_pj_context;
3061
0
        }
3062
225k
        return psGTIF->pj_context;
3063
225k
    }
3064
0
    psGTIF->pj_context = proj_context_create();
3065
0
    psGTIF->own_pj_context = psGTIF->pj_context != NULL;
3066
0
    if( out_gtif_own_pj_context )
3067
0
    {
3068
0
        *out_gtif_own_pj_context = psGTIF->own_pj_context;
3069
0
    }
3070
0
    return psGTIF->pj_context;
3071
225k
}
3072
3073
3074
void GTIFDeaccessCSV( void )
3075
0
{
3076
    /* No operation */
3077
0
}
3078
3079
#ifndef GDAL_COMPILATION
3080
void SetCSVFilenameHook( const char *(*CSVFileOverride)(const char *) )
3081
{
3082
    (void)CSVFileOverride;
3083
    /* No operation */
3084
}
3085
#endif