Coverage Report

Created: 2025-06-09 08:44

/src/gdal/frmts/gtiff/libgeotiff/geo_normalize.c
Line
Count
Source (jump to first uncovered line)
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
17.1k
#define EPSGNatOriginLat         8801
39
17.1k
#define EPSGNatOriginLong        8802
40
104k
#define EPSGNatOriginScaleFactor 8805
41
17.1k
#define EPSGFalseEasting         8806
42
17.1k
#define EPSGFalseNorthing        8807
43
74
#define EPSGProjCenterLat        8811
44
74
#define EPSGProjCenterLong       8812
45
74
#define EPSGAzimuth              8813
46
43.8k
#define EPSGAngleRectifiedToSkewedGrid 8814
47
83.0k
#define EPSGInitialLineScaleFactor 8815
48
128
#define EPSGProjCenterEasting    8816
49
98
#define EPSGProjCenterNorthing   8817
50
#define EPSGPseudoStdParallelLat 8818
51
39.1k
#define EPSGPseudoStdParallelScaleFactor 8819
52
2.54k
#define EPSGFalseOriginLat       8821
53
2.54k
#define EPSGFalseOriginLong      8822
54
2.59k
#define EPSGStdParallel1Lat      8823
55
2.54k
#define EPSGStdParallel2Lat      8824
56
2.54k
#define EPSGFalseOriginEasting   8826
57
2.54k
#define EPSGFalseOriginNorthing  8827
58
#define EPSGSphericalOriginLat   8828
59
#define EPSGSphericalOriginLong  8829
60
#define EPSGInitialLongitude     8830
61
#define EPSGZoneWidth            8831
62
32
#define EPSGLatOfStdParallel     8832
63
32
#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
121k
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
434k
{
101
434k
    unsigned short sVal;
102
434k
    if( GTIFKeyGetSHORT(gtif, key, &sVal, 0, 1) == 1 )
103
122k
    {
104
122k
        memcpy(pnVal, &sVal, 2);
105
122k
        return 1;
106
122k
    }
107
312k
    return 0;
108
434k
}
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
41.8k
{
120
41.8k
    int         nDatum;
121
41.8k
    int         nZone;
122
123
    /* Deal with a few well known CRS */
124
41.8k
    int Proj = GTIFPCSToMapSys( nPCSCode, &nDatum, &nZone );
125
41.8k
    if ((Proj == MapSys_UTM_North || Proj == MapSys_UTM_South) &&
126
41.8k
        nDatum != KvUserDefined)
127
14.1k
    {
128
14.1k
        const char* pszDatumName = NULL;
129
14.1k
        switch (nDatum)
130
14.1k
        {
131
615
            case GCS_NAD27: pszDatumName = "NAD27"; break;
132
140
            case GCS_NAD83: pszDatumName = "NAD83"; break;
133
1.79k
            case GCS_WGS_72: pszDatumName = "WGS 72"; break;
134
11.1k
            case GCS_WGS_72BE: pszDatumName = "WGS 72BE"; break;
135
463
            case GCS_WGS_84: pszDatumName = "WGS 84"; break;
136
0
            default: break;
137
14.1k
        }
138
139
14.1k
        if (pszDatumName)
140
14.1k
        {
141
14.1k
            if (ppszEPSGName)
142
7.04k
            {
143
7.04k
                char szEPSGName[64];
144
7.04k
                snprintf(
145
7.04k
                    szEPSGName, sizeof(szEPSGName), "%s / UTM zone %d%c",
146
7.04k
                    pszDatumName, nZone,
147
7.04k
                    (Proj == MapSys_UTM_North) ? 'N' : 'S');
148
7.04k
                *ppszEPSGName = CPLStrdup(szEPSGName);
149
7.04k
            }
150
151
14.1k
            if (pnProjOp)
152
7.15k
                *pnProjOp = (short) (((Proj == MapSys_UTM_North) ? Proj_UTM_zone_1N - 1 : Proj_UTM_zone_1S - 1) + nZone);
153
154
14.1k
            if (pnUOMLengthCode)
155
7.15k
                *pnUOMLengthCode = 9001; /* Linear_Meter */
156
157
14.1k
            if (pnGeogCS)
158
7.15k
                *pnGeogCS = (short) nDatum;
159
160
14.1k
            return TRUE;
161
14.1k
        }
162
14.1k
    }
163
164
27.6k
    if( nPCSCode == KvUserDefined )
165
0
        return FALSE;
166
167
27.6k
    PJ_CONTEXT* ctx = (PJ_CONTEXT*)ctxIn;
168
27.6k
    {
169
27.6k
        char szCode[12];
170
171
27.6k
        snprintf(szCode, sizeof(szCode), "%d", nPCSCode);
172
27.6k
        PJ* proj_crs = proj_create_from_database(
173
27.6k
            ctx, "EPSG", szCode, PJ_CATEGORY_CRS, 0, NULL);
174
27.6k
        if( !proj_crs )
175
13.6k
        {
176
13.6k
            return FALSE;
177
13.6k
        }
178
179
13.9k
        if( proj_get_type(proj_crs) != PJ_TYPE_PROJECTED_CRS )
180
1.52k
        {
181
1.52k
            proj_destroy(proj_crs);
182
1.52k
            return FALSE;
183
1.52k
        }
184
185
12.4k
        if( ppszEPSGName )
186
6.19k
        {
187
6.19k
            const char* pszName = proj_get_name(proj_crs);
188
6.19k
            if( !pszName )
189
0
            {
190
                // shouldn't happen
191
0
                proj_destroy(proj_crs);
192
0
                return FALSE;
193
0
            }
194
6.19k
            *ppszEPSGName = CPLStrdup(pszName);
195
6.19k
        }
196
197
12.4k
        if( pnProjOp )
198
6.27k
        {
199
6.27k
            PJ* conversion = proj_crs_get_coordoperation(
200
6.27k
                ctx, proj_crs);
201
6.27k
            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
6.27k
            {
209
6.27k
                const char* pszConvCode = proj_get_id_code(conversion, 0);
210
6.27k
                assert( pszConvCode );
211
6.27k
                *pnProjOp = (short) atoi(pszConvCode);
212
6.27k
            }
213
214
0
            proj_destroy(conversion);
215
6.27k
        }
216
217
12.4k
        if( pnUOMLengthCode )
218
6.27k
        {
219
6.27k
            PJ* coordSys = proj_crs_get_coordinate_system(
220
6.27k
                ctx, proj_crs);
221
6.27k
            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
6.27k
            {
229
6.27k
                const char* pszUnitCode = NULL;
230
6.27k
                if( !proj_cs_get_axis_info(
231
6.27k
                    ctx, coordSys, 0,
232
6.27k
                    NULL, /* name */
233
6.27k
                    NULL, /* abbreviation*/
234
6.27k
                    NULL, /* direction */
235
6.27k
                    NULL, /* conversion factor */
236
6.27k
                    NULL, /* unit name */
237
6.27k
                    NULL, /* unit auth name (should be EPSG) */
238
6.27k
                    &pszUnitCode) || pszUnitCode == NULL )
239
0
                {
240
0
                    proj_destroy(coordSys);
241
0
                    return FALSE;
242
0
                }
243
6.27k
                *pnUOMLengthCode = (short) atoi(pszUnitCode);
244
6.27k
                proj_destroy(coordSys);
245
6.27k
            }
246
6.27k
        }
247
248
12.4k
        if( pnGeogCS )
249
6.27k
        {
250
6.27k
            PJ* geod_crs = proj_crs_get_geodetic_crs(ctx, proj_crs);
251
6.27k
            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
6.27k
            {
259
6.27k
                const char* pszGeodCode = proj_get_id_code(geod_crs, 0);
260
6.27k
                assert( pszGeodCode );
261
6.27k
                *pnGeogCS = (short) atoi(pszGeodCode);
262
6.27k
            }
263
264
0
            proj_destroy(geod_crs);
265
6.27k
        }
266
267
268
12.4k
        proj_destroy(proj_crs);
269
12.4k
        return TRUE;
270
12.4k
    }
271
12.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
23.5k
{
295
23.5k
    if( nUOMAngle == 9110 )    /* DDD.MMSSsss */
296
54
    {
297
54
        if( dfAngle > -999.9 && dfAngle < 999.9 )
298
48
        {
299
48
            char  szAngleString[32];
300
301
48
            snprintf(szAngleString, sizeof(szAngleString), "%12.7f", dfAngle);
302
48
            dfAngle = GTIFAngleStringToDD( szAngleString, nUOMAngle );
303
48
        }
304
54
    }
305
23.4k
    else if ( nUOMAngle != KvUserDefined )
306
17.3k
    {
307
17.3k
        double    dfInDegrees = 1.0;
308
309
17.3k
        GTIFGetUOMAngleInfo( nUOMAngle, NULL, &dfInDegrees );
310
17.3k
        dfAngle = dfAngle * dfInDegrees;
311
17.3k
    }
312
313
23.5k
    return dfAngle;
314
23.5k
}
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
48
{
325
48
    double  dfAngle;
326
327
48
    if( nUOMAngle == 9110 )    /* DDD.MMSSsss */
328
48
    {
329
48
        char  *pszDecimal;
330
331
48
        dfAngle = ABS(atoi(pszAngle));
332
48
        pszDecimal = strchr(pszAngle,'.');
333
48
        if( pszDecimal != NULL && strlen(pszDecimal) > 1 )
334
48
        {
335
48
            char  szMinutes[3];
336
48
            char  szSeconds[64];
337
338
48
            szMinutes[0] = pszDecimal[1];
339
48
            if( pszDecimal[2] >= '0' && pszDecimal[2] <= '9' )
340
48
                szMinutes[1] = pszDecimal[2];
341
0
            else
342
0
                szMinutes[1] = '0';
343
344
48
            szMinutes[2] = '\0';
345
48
            dfAngle += atoi(szMinutes) / 60.0;
346
347
48
            if( strlen(pszDecimal) > 3 )
348
48
            {
349
48
                szSeconds[0] = pszDecimal[3];
350
48
                if( pszDecimal[4] >= '0' && pszDecimal[4] <= '9' )
351
48
                {
352
48
                    szSeconds[1] = pszDecimal[4];
353
48
                    szSeconds[2] = '.';
354
48
                    strncpy( szSeconds+3, pszDecimal + 5, sizeof(szSeconds) - 3 );
355
48
                    szSeconds[sizeof(szSeconds) - 1] = 0;
356
48
                }
357
0
                else
358
0
                {
359
0
                    szSeconds[1] = '0';
360
0
                    szSeconds[2] = '\0';
361
0
                }
362
48
                dfAngle += GTIFAtof(szSeconds) / 3600.0;
363
48
            }
364
48
        }
365
366
48
        if( pszAngle[0] == '-' )
367
1
            dfAngle *= -1;
368
48
    }
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
48
    return dfAngle;
394
48
}
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
58.7k
{
408
58.7k
    int   nDatum=0;
409
410
/* -------------------------------------------------------------------- */
411
/*      Handle some "well known" GCS codes directly                     */
412
/* -------------------------------------------------------------------- */
413
58.7k
    const char * pszName = NULL;
414
58.7k
    const int nPM = PM_Greenwich;
415
58.7k
    const int nUOMAngle = Angular_DMS_Hemisphere;
416
58.7k
    if( nGCSCode == GCS_NAD27 )
417
3.12k
    {
418
3.12k
        nDatum = Datum_North_American_Datum_1927;
419
3.12k
        pszName = "NAD27";
420
3.12k
    }
421
55.6k
    else if( nGCSCode == GCS_NAD83 )
422
7.11k
    {
423
7.11k
        nDatum = Datum_North_American_Datum_1983;
424
7.11k
        pszName = "NAD83";
425
7.11k
    }
426
48.5k
    else if( nGCSCode == GCS_WGS_84 )
427
1.50k
    {
428
1.50k
        nDatum = Datum_WGS84;
429
1.50k
        pszName = "WGS 84";
430
1.50k
    }
431
47.0k
    else if( nGCSCode == GCS_WGS_72 )
432
1.76k
    {
433
1.76k
        nDatum = Datum_WGS72;
434
1.76k
        pszName = "WGS 72";
435
1.76k
    }
436
45.2k
    else if ( nGCSCode == KvUserDefined )
437
20.5k
    {
438
20.5k
        return FALSE;
439
20.5k
    }
440
441
38.2k
    if (pszName != NULL)
442
13.5k
    {
443
13.5k
        if( ppszName != NULL )
444
6.62k
            *ppszName = CPLStrdup( pszName );
445
13.5k
        if( pnDatum != NULL )
446
6.88k
            *pnDatum = (short) nDatum;
447
13.5k
        if( pnPM != NULL )
448
6.88k
            *pnPM = (short) nPM;
449
13.5k
        if( pnUOMAngle != NULL )
450
6.88k
            *pnUOMAngle = (short) nUOMAngle;
451
452
13.5k
        return TRUE;
453
13.5k
    }
454
455
/* -------------------------------------------------------------------- */
456
/*      Search the database.                                            */
457
/* -------------------------------------------------------------------- */
458
459
24.7k
    PJ_CONTEXT* ctx = (PJ_CONTEXT*)ctxIn;
460
24.7k
    {
461
24.7k
        char szCode[12];
462
463
24.7k
        snprintf(szCode, sizeof(szCode), "%d", nGCSCode);
464
24.7k
        PJ* geod_crs = proj_create_from_database(
465
24.7k
            ctx, "EPSG", szCode, PJ_CATEGORY_CRS, 0, NULL);
466
24.7k
        if( !geod_crs )
467
1.20k
        {
468
1.20k
            return FALSE;
469
1.20k
        }
470
471
23.4k
        {
472
23.4k
            const int objType = proj_get_type(geod_crs);
473
23.4k
            if( objType != PJ_TYPE_GEODETIC_CRS &&
474
23.4k
                objType != PJ_TYPE_GEOCENTRIC_CRS &&
475
23.4k
                objType != PJ_TYPE_GEOGRAPHIC_2D_CRS &&
476
23.4k
                objType != PJ_TYPE_GEOGRAPHIC_3D_CRS )
477
2.58k
            {
478
2.58k
                proj_destroy(geod_crs);
479
2.58k
                return FALSE;
480
2.58k
            }
481
23.4k
        }
482
483
20.9k
        if( ppszName )
484
10.4k
        {
485
10.4k
            pszName = proj_get_name(geod_crs);
486
10.4k
            if( !pszName )
487
0
            {
488
                // shouldn't happen
489
0
                proj_destroy(geod_crs);
490
0
                return FALSE;
491
0
            }
492
10.4k
            *ppszName = CPLStrdup(pszName);
493
10.4k
        }
494
495
20.9k
        if( pnDatum )
496
10.4k
        {
497
10.4k
#if PROJ_VERSION_MAJOR >= 8
498
10.4k
            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
10.4k
            if( !datum )
503
0
            {
504
0
                proj_destroy(geod_crs);
505
0
                return FALSE;
506
0
            }
507
508
10.4k
            {
509
10.4k
                const char* pszDatumCode = proj_get_id_code(datum, 0);
510
10.4k
                assert( pszDatumCode );
511
10.4k
                *pnDatum = (short) atoi(pszDatumCode);
512
10.4k
            }
513
514
0
            proj_destroy(datum);
515
10.4k
        }
516
517
20.9k
        if( pnPM )
518
10.4k
        {
519
10.4k
            PJ* pm = proj_get_prime_meridian(ctx, geod_crs);
520
10.4k
            if( !pm )
521
0
            {
522
0
                proj_destroy(geod_crs);
523
0
                return FALSE;
524
0
            }
525
526
10.4k
            {
527
10.4k
                const char* pszPMCode = proj_get_id_code(pm, 0);
528
10.4k
                assert( pszPMCode );
529
10.4k
                *pnPM = (short) atoi(pszPMCode);
530
10.4k
            }
531
532
0
            proj_destroy(pm);
533
10.4k
        }
534
535
20.9k
        if( pnUOMAngle )
536
10.4k
        {
537
10.4k
            PJ* coordSys = proj_crs_get_coordinate_system(
538
10.4k
                ctx, geod_crs);
539
10.4k
            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
10.4k
            {
547
10.4k
                const char* pszUnitCode = NULL;
548
10.4k
                if( !proj_cs_get_axis_info(
549
10.4k
                    ctx, coordSys, 0,
550
10.4k
                    NULL, /* name */
551
10.4k
                    NULL, /* abbreviation*/
552
10.4k
                    NULL, /* direction */
553
10.4k
                    NULL, /* conversion factor */
554
10.4k
                    NULL, /* unit name */
555
10.4k
                    NULL, /* unit auth name (should be EPSG) */
556
10.4k
                    &pszUnitCode) || pszUnitCode == NULL )
557
0
                {
558
0
                    proj_destroy(coordSys);
559
0
                    return FALSE;
560
0
                }
561
10.4k
                *pnUOMAngle = (short) atoi(pszUnitCode);
562
10.4k
                proj_destroy(coordSys);
563
10.4k
            }
564
10.4k
        }
565
566
20.9k
        proj_destroy(geod_crs);
567
20.9k
        return TRUE;
568
20.9k
    }
569
20.9k
}
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
56.3k
{
595
56.3k
    PJ_CONTEXT* ctx = (PJ_CONTEXT*)ctxIn;
596
/* -------------------------------------------------------------------- */
597
/*      Try some well known ellipsoids.                                 */
598
/* -------------------------------------------------------------------- */
599
56.3k
    double  dfSemiMajor=0.0;
600
56.3k
    double     dfInvFlattening=0.0, dfSemiMinor=0.0;
601
56.3k
    const char *pszName = NULL;
602
603
56.3k
    if( nEllipseCode == Ellipse_Clarke_1866 )
604
3.06k
    {
605
3.06k
        pszName = "Clarke 1866";
606
3.06k
        dfSemiMajor = 6378206.4;
607
3.06k
        dfSemiMinor = 6356583.8;
608
3.06k
        dfInvFlattening = 0.0;
609
3.06k
    }
610
53.3k
    else if( nEllipseCode == Ellipse_GRS_1980 )
611
10.1k
    {
612
10.1k
        pszName = "GRS 1980";
613
10.1k
        dfSemiMajor = 6378137.0;
614
10.1k
        dfSemiMinor = 0.0;
615
10.1k
        dfInvFlattening = 298.257222101;
616
10.1k
    }
617
43.2k
    else if( nEllipseCode == Ellipse_WGS_84 )
618
1.70k
    {
619
1.70k
        pszName = "WGS 84";
620
1.70k
        dfSemiMajor = 6378137.0;
621
1.70k
        dfSemiMinor = 0.0;
622
1.70k
        dfInvFlattening = 298.257223563;
623
1.70k
    }
624
41.5k
    else if( nEllipseCode == 7043 )
625
11.5k
    {
626
11.5k
        pszName = "WGS 72";
627
11.5k
        dfSemiMajor = 6378135.0;
628
11.5k
        dfSemiMinor = 0.0;
629
11.5k
        dfInvFlattening = 298.26;
630
11.5k
    }
631
632
56.3k
    if (pszName != NULL)
633
26.3k
    {
634
26.3k
        if( dfSemiMinor == 0.0 )
635
23.3k
            dfSemiMinor = dfSemiMajor * (1 - 1.0/dfInvFlattening);
636
637
26.3k
        if( pdfSemiMinor != NULL )
638
13.4k
            *pdfSemiMinor = dfSemiMinor;
639
26.3k
        if( pdfSemiMajor != NULL )
640
13.4k
            *pdfSemiMajor = dfSemiMajor;
641
26.3k
        if( ppszName != NULL )
642
12.9k
            *ppszName = CPLStrdup( pszName );
643
644
26.3k
        return TRUE;
645
26.3k
    }
646
647
29.9k
    if( nEllipseCode == KvUserDefined )
648
22.4k
        return FALSE;
649
650
/* -------------------------------------------------------------------- */
651
/*      Search the database.                                            */
652
/* -------------------------------------------------------------------- */
653
7.49k
    {
654
7.49k
        char szCode[12];
655
656
7.49k
        snprintf(szCode, sizeof(szCode), "%d", nEllipseCode);
657
7.49k
        PJ* ellipsoid = proj_create_from_database(
658
7.49k
            ctx, "EPSG", szCode, PJ_CATEGORY_ELLIPSOID, 0, NULL);
659
7.49k
        if( !ellipsoid )
660
1.88k
        {
661
1.88k
            return FALSE;
662
1.88k
        }
663
664
5.61k
        if( ppszName )
665
2.79k
        {
666
2.79k
            pszName = proj_get_name(ellipsoid);
667
2.79k
            if( !pszName )
668
0
            {
669
                // shouldn't happen
670
0
                proj_destroy(ellipsoid);
671
0
                return FALSE;
672
0
            }
673
2.79k
            *ppszName = CPLStrdup(pszName);
674
2.79k
        }
675
676
5.61k
        proj_ellipsoid_get_parameters(
677
5.61k
            ctx, ellipsoid, pdfSemiMajor, pdfSemiMinor, NULL, NULL);
678
679
5.61k
        proj_destroy(ellipsoid);
680
681
5.61k
        return TRUE;
682
5.61k
    }
683
5.61k
}
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
56.4k
{
707
/* -------------------------------------------------------------------- */
708
/*      Use a special short cut for Greenwich, since it is so common.   */
709
/* -------------------------------------------------------------------- */
710
56.4k
    if( nPMCode == PM_Greenwich )
711
34.1k
    {
712
34.1k
        if( pdfOffset != NULL )
713
17.2k
            *pdfOffset = 0.0;
714
34.1k
        if( ppszName != NULL )
715
16.9k
            *ppszName = CPLStrdup( "Greenwich" );
716
34.1k
        return TRUE;
717
34.1k
    }
718
719
720
22.3k
    if( nPMCode == KvUserDefined )
721
21.7k
        return FALSE;
722
723
/* -------------------------------------------------------------------- */
724
/*      Search the database.                                            */
725
/* -------------------------------------------------------------------- */
726
589
    PJ_CONTEXT* ctx = (PJ_CONTEXT*)ctxIn;
727
589
    {
728
589
        char szCode[12];
729
730
589
        snprintf(szCode, sizeof(szCode), "%d", nPMCode);
731
589
        PJ* pm = proj_create_from_database(
732
589
            ctx, "EPSG", szCode, PJ_CATEGORY_PRIME_MERIDIAN, 0, NULL);
733
589
        if( !pm )
734
246
        {
735
246
            return FALSE;
736
246
        }
737
738
343
        if( ppszName )
739
170
        {
740
170
            const char* pszName = proj_get_name(pm);
741
170
            if( !pszName )
742
0
            {
743
                // shouldn't happen
744
0
                proj_destroy(pm);
745
0
                return FALSE;
746
0
            }
747
170
            *ppszName = CPLStrdup(pszName);
748
170
        }
749
750
343
        if( pdfOffset )
751
173
        {
752
173
            double conv_factor = 0;
753
173
            proj_prime_meridian_get_parameters(
754
173
                ctx, pm, pdfOffset, &conv_factor, NULL);
755
173
            *pdfOffset *= conv_factor * 180.0 / M_PI;
756
173
        }
757
758
343
        proj_destroy(pm);
759
760
343
        return TRUE;
761
343
    }
762
343
}
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
55.8k
{
783
55.8k
    const char* pszName = NULL;
784
55.8k
    int   nEllipsoid = 0;
785
786
/* -------------------------------------------------------------------- */
787
/*      Handle a few built-in datums.                                   */
788
/* -------------------------------------------------------------------- */
789
55.8k
    if( nDatumCode == Datum_North_American_Datum_1927 )
790
3.03k
    {
791
3.03k
        nEllipsoid = Ellipse_Clarke_1866;
792
3.03k
        pszName = "North American Datum 1927";
793
3.03k
    }
794
52.8k
    else if( nDatumCode == Datum_North_American_Datum_1983 )
795
6.94k
    {
796
6.94k
        nEllipsoid = Ellipse_GRS_1980;
797
6.94k
        pszName = "North American Datum 1983";
798
6.94k
    }
799
45.8k
    else if( nDatumCode == Datum_WGS84 )
800
1.65k
    {
801
1.65k
        nEllipsoid = Ellipse_WGS_84;
802
1.65k
        pszName = "World Geodetic System 1984";
803
1.65k
    }
804
44.2k
    else if( nDatumCode == Datum_WGS72 )
805
1.71k
    {
806
1.71k
        nEllipsoid = 7043; /* WGS72 */
807
1.71k
        pszName = "World Geodetic System 1972";
808
1.71k
    }
809
810
55.8k
    if (pszName != NULL)
811
13.3k
    {
812
13.3k
        if( pnEllipsoid != NULL )
813
6.85k
            *pnEllipsoid = (short) nEllipsoid;
814
815
13.3k
        if( ppszName != NULL )
816
6.50k
            *ppszName = CPLStrdup( pszName );
817
818
13.3k
        return TRUE;
819
13.3k
    }
820
821
42.5k
    if( nDatumCode == KvUserDefined )
822
17.0k
        return FALSE;
823
824
/* -------------------------------------------------------------------- */
825
/*      Search the database.                                            */
826
/* -------------------------------------------------------------------- */
827
25.4k
    PJ_CONTEXT* ctx = (PJ_CONTEXT*)ctxIn;
828
25.4k
    {
829
25.4k
        char szCode[12];
830
831
25.4k
        snprintf(szCode, sizeof(szCode), "%d", nDatumCode);
832
25.4k
        PJ* datum = proj_create_from_database(
833
25.4k
            ctx, "EPSG", szCode, PJ_CATEGORY_DATUM, 0, NULL);
834
25.4k
        if( !datum )
835
6.74k
        {
836
6.74k
            return FALSE;
837
6.74k
        }
838
839
18.7k
        const PJ_TYPE pjType = proj_get_type(datum);
840
18.7k
        if( pjType != PJ_TYPE_GEODETIC_REFERENCE_FRAME &&
841
18.7k
            pjType != PJ_TYPE_DYNAMIC_GEODETIC_REFERENCE_FRAME )
842
49
        {
843
49
            proj_destroy(datum);
844
49
            return FALSE;
845
49
        }
846
847
18.6k
        if( ppszName )
848
9.29k
        {
849
9.29k
            pszName = proj_get_name(datum);
850
9.29k
            if( !pszName )
851
0
            {
852
                // shouldn't happen
853
0
                proj_destroy(datum);
854
0
                return FALSE;
855
0
            }
856
9.29k
            *ppszName = CPLStrdup(pszName);
857
9.29k
        }
858
859
18.6k
        if( pnEllipsoid )
860
9.38k
        {
861
9.38k
            PJ* ellipsoid = proj_get_ellipsoid(ctx, datum);
862
9.38k
            if( !ellipsoid )
863
0
            {
864
0
                proj_destroy(datum);
865
0
                return FALSE;
866
0
            }
867
868
9.38k
            {
869
9.38k
                const char* pszEllipsoidCode = proj_get_id_code(
870
9.38k
                    ellipsoid, 0);
871
9.38k
                assert( pszEllipsoidCode );
872
9.38k
                *pnEllipsoid = (short) atoi(pszEllipsoidCode);
873
9.38k
            }
874
875
0
            proj_destroy(ellipsoid);
876
9.38k
        }
877
878
18.6k
        proj_destroy(datum);
879
880
18.6k
        return TRUE;
881
18.6k
    }
882
18.6k
}
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
54.1k
{
906
/* -------------------------------------------------------------------- */
907
/*      We short cut meter to save work and avoid failure for missing   */
908
/*      in the most common cases.               */
909
/* -------------------------------------------------------------------- */
910
54.1k
    if( nUOMLengthCode == 9001 )
911
27.7k
    {
912
27.7k
        if( ppszUOMName != NULL )
913
14.3k
            *ppszUOMName = CPLStrdup( "metre" );
914
27.7k
        if( pdfInMeters != NULL )
915
13.4k
            *pdfInMeters = 1.0;
916
917
27.7k
        return TRUE;
918
27.7k
    }
919
920
26.3k
    if( nUOMLengthCode == 9002 )
921
2.35k
    {
922
2.35k
        if( ppszUOMName != NULL )
923
1.40k
            *ppszUOMName = CPLStrdup( "foot" );
924
2.35k
        if( pdfInMeters != NULL )
925
955
            *pdfInMeters = 0.3048;
926
927
2.35k
        return TRUE;
928
2.35k
    }
929
930
24.0k
    if( nUOMLengthCode == 9003 )
931
2.05k
    {
932
2.05k
        if( ppszUOMName != NULL )
933
1.02k
            *ppszUOMName = CPLStrdup( "US survey foot" );
934
2.05k
        if( pdfInMeters != NULL )
935
1.03k
            *pdfInMeters = 12.0 / 39.37;
936
937
2.05k
        return TRUE;
938
2.05k
    }
939
940
21.9k
    if( nUOMLengthCode == KvUserDefined )
941
5.16k
        return FALSE;
942
943
/* -------------------------------------------------------------------- */
944
/*      Search the units database for this unit.  If we don't find      */
945
/*      it return failure.                                              */
946
/* -------------------------------------------------------------------- */
947
16.8k
    char szCode[12];
948
16.8k
    const char* pszName = NULL;
949
950
16.8k
    snprintf(szCode, sizeof(szCode), "%d", nUOMLengthCode);
951
16.8k
    PJ_CONTEXT* ctx = (PJ_CONTEXT*)ctxIn;
952
16.8k
    if( !proj_uom_get_info_from_database(
953
16.8k
        ctx, "EPSG", szCode, &pszName, pdfInMeters,  NULL) )
954
15.9k
    {
955
15.9k
        return FALSE;
956
15.9k
    }
957
830
    if( ppszUOMName )
958
457
    {
959
457
        *ppszUOMName = CPLStrdup(pszName);
960
457
    }
961
830
    return TRUE;
962
16.8k
}
963
964
int GTIFGetUOMLengthInfo( int nUOMLengthCode,
965
                          char **ppszUOMName,
966
                          double * pdfInMeters )
967
968
9.02k
{
969
9.02k
    PJ_CONTEXT* ctx = proj_context_create();
970
9.02k
    const int ret = GTIFGetUOMLengthInfoEx(
971
9.02k
        ctx, nUOMLengthCode, ppszUOMName, pdfInMeters);
972
9.02k
    proj_context_destroy(ctx);
973
9.02k
    return ret;
974
9.02k
}
975
976
/************************************************************************/
977
/*                        GTIFGetUOMAngleInfo()                         */
978
/************************************************************************/
979
980
int GTIFGetUOMAngleInfoEx( void* ctxIn,
981
                           int nUOMAngleCode,
982
                           char **ppszUOMName,
983
                           double * pdfInDegrees )
984
985
91.4k
{
986
91.4k
    const char  *pszUOMName = NULL;
987
91.4k
    double  dfInDegrees = 1.0;
988
989
91.4k
    switch( nUOMAngleCode )
990
91.4k
    {
991
2.66k
      case 9101:
992
2.66k
        pszUOMName = "radian";
993
2.66k
        dfInDegrees = 180.0 / M_PI;
994
2.66k
        break;
995
996
3.64k
      case 9102:
997
3.68k
      case 9107:
998
6.15k
      case 9108:
999
6.26k
      case 9110:
1000
11.3k
      case 9122:
1001
11.3k
        pszUOMName = "degree";
1002
11.3k
        dfInDegrees = 1.0;
1003
11.3k
        break;
1004
1005
37.8k
      case 9103:
1006
37.8k
        pszUOMName = "arc-minute";
1007
37.8k
        dfInDegrees = 1 / 60.0;
1008
37.8k
        break;
1009
1010
1.11k
      case 9104:
1011
1.11k
        pszUOMName = "arc-second";
1012
1.11k
        dfInDegrees = 1 / 3600.0;
1013
1.11k
        break;
1014
1015
1.67k
      case 9105:
1016
1.67k
        pszUOMName = "grad";
1017
1.67k
        dfInDegrees = 180.0 / 200.0;
1018
1.67k
        break;
1019
1020
16.4k
      case 9106:
1021
16.4k
        pszUOMName = "gon";
1022
16.4k
        dfInDegrees = 180.0 / 200.0;
1023
16.4k
        break;
1024
1025
25
      case 9109:
1026
25
        pszUOMName = "microradian";
1027
25
        dfInDegrees = 180.0 / (M_PI * 1000000.0);
1028
25
        break;
1029
1030
20.4k
      default:
1031
20.4k
        break;
1032
91.4k
    }
1033
1034
91.4k
    if (pszUOMName)
1035
71.0k
    {
1036
71.0k
        if( ppszUOMName != NULL )
1037
27.9k
        {
1038
27.9k
            *ppszUOMName = CPLStrdup( pszUOMName );
1039
27.9k
        }
1040
1041
71.0k
        if( pdfInDegrees != NULL )
1042
71.0k
            *pdfInDegrees = dfInDegrees;
1043
1044
71.0k
        return TRUE;
1045
71.0k
    }
1046
1047
20.4k
    if( nUOMAngleCode == KvUserDefined )
1048
5.43k
        return FALSE;
1049
1050
/* -------------------------------------------------------------------- */
1051
/*      Search the units database for this unit.  If we don't find      */
1052
/*      it return failure.                                              */
1053
/* -------------------------------------------------------------------- */
1054
14.9k
    char szCode[12];
1055
14.9k
    const char* pszName = NULL;
1056
14.9k
    double dfConvFactorToRadians = 0;
1057
1058
14.9k
    snprintf(szCode, sizeof(szCode), "%d", nUOMAngleCode);
1059
14.9k
    PJ_CONTEXT* ctx = (PJ_CONTEXT*)ctxIn;
1060
14.9k
    if( !proj_uom_get_info_from_database(
1061
14.9k
        ctx, "EPSG", szCode, &pszName, &dfConvFactorToRadians, NULL) )
1062
11.9k
    {
1063
11.9k
        return FALSE;
1064
11.9k
    }
1065
3.01k
    if( ppszUOMName )
1066
1.13k
    {
1067
1.13k
        *ppszUOMName = CPLStrdup(pszName);
1068
1.13k
    }
1069
3.01k
    if( pdfInDegrees )
1070
3.01k
    {
1071
3.01k
        *pdfInDegrees = dfConvFactorToRadians * 180.0 / M_PI;
1072
3.01k
    }
1073
3.01k
    return TRUE;
1074
14.9k
}
1075
1076
int GTIFGetUOMAngleInfo( int nUOMAngleCode,
1077
                         char **ppszUOMName,
1078
                         double * pdfInDegrees )
1079
1080
17.3k
{
1081
17.3k
    PJ_CONTEXT* ctx = proj_context_create();
1082
17.3k
    const int ret = GTIFGetUOMAngleInfoEx(
1083
17.3k
        ctx, nUOMAngleCode, ppszUOMName, pdfInDegrees);
1084
17.3k
    proj_context_destroy(ctx);
1085
17.3k
    return ret;
1086
17.3k
}
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
44.8k
{
1098
44.8k
    switch( nEPSG )
1099
44.8k
    {
1100
666
      case 9801:
1101
666
        return CT_LambertConfConic_1SP;
1102
1103
3.57k
      case 9802:
1104
3.57k
        return CT_LambertConfConic_2SP;
1105
1106
0
      case 9803:
1107
0
        return CT_LambertConfConic_2SP; /* Belgian variant not supported */
1108
1109
18
      case 9804:
1110
18
        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
21
      case 9841:
1120
21
        return CT_Mercator;  /* 1SP and 2SP not differentiated */
1121
1122
      /* Google Mercator For EPSG:3857 */
1123
123
      case 1024:
1124
123
        return CT_Mercator;  /* 1SP and 2SP not differentiated */
1125
1126
87
      case 9806:
1127
87
        return CT_CassiniSoldner;
1128
1129
27.7k
      case 9807:
1130
27.7k
        return CT_TransverseMercator;
1131
1132
135
      case 9808:
1133
135
        return CT_TransvMercator_SouthOriented;
1134
1135
72
      case 9809:
1136
72
        return CT_ObliqueStereographic;
1137
1138
333
      case 9810:
1139
429
      case 9829: /* variant B not quite the same - not sure how to handle */
1140
429
        return CT_PolarStereographic;
1141
1142
0
      case 9811:
1143
0
        return CT_NewZealandMapGrid;
1144
1145
90
      case 9812:
1146
90
        return CT_ObliqueMercator; /* is hotine actually different? */
1147
1148
9
      case 9813:
1149
9
        return CT_ObliqueMercator_Laborde;
1150
1151
0
      case 9814:
1152
0
        return CT_ObliqueMercator_Rosenmund; /* swiss  */
1153
1154
12
      case 9815:
1155
12
        return CT_HotineObliqueMercatorAzimuthCenter;
1156
1157
0
      case 9816: /* tunesia mining grid has no counterpart */
1158
0
        return KvUserDefined;
1159
1160
9
      case 9818:
1161
9
        return CT_Polyconic;
1162
1163
114
      case 9820:
1164
117
      case 1027:
1165
117
        return CT_LambertAzimEqualArea;
1166
1167
240
      case 9822:
1168
240
        return CT_AlbersEqualArea;
1169
1170
9
      case 9834:
1171
9
        return CT_CylindricalEqualArea;
1172
1173
6
      case 1028:
1174
30
      case 1029:
1175
63
      case 9823: /* spherical */
1176
75
      case 9842: /* elliptical */
1177
75
        return CT_Equirectangular;
1178
1179
11.3k
      default: /* use the EPSG code for other methods */
1180
11.3k
        return nEPSG;
1181
44.8k
    }
1182
44.8k
}
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
25.5k
{
1198
25.5k
    int anWorkingDummy[7];
1199
1200
25.5k
    if( panEPSGCodes == NULL )
1201
19.2k
        panEPSGCodes = anWorkingDummy;
1202
25.5k
    if( panProjParamId == NULL )
1203
6.25k
        panProjParamId = anWorkingDummy;
1204
1205
25.5k
    memset( panEPSGCodes, 0, sizeof(int) * 7 );
1206
1207
    /* psDefn->nParms = 7; */
1208
1209
25.5k
    switch( nCTProjection )
1210
25.5k
    {
1211
58
      case CT_CassiniSoldner:
1212
58
      case CT_NewZealandMapGrid:
1213
64
      case CT_Polyconic:
1214
64
        panProjParamId[0] = ProjNatOriginLatGeoKey;
1215
64
        panProjParamId[1] = ProjNatOriginLongGeoKey;
1216
64
        panProjParamId[5] = ProjFalseEastingGeoKey;
1217
64
        panProjParamId[6] = ProjFalseNorthingGeoKey;
1218
1219
64
        panEPSGCodes[0] = EPSGNatOriginLat;
1220
64
        panEPSGCodes[1] = EPSGNatOriginLong;
1221
64
        panEPSGCodes[5] = EPSGFalseEasting;
1222
64
        panEPSGCodes[6] = EPSGFalseNorthing;
1223
64
        return TRUE;
1224
1225
60
      case CT_ObliqueMercator:
1226
68
      case CT_HotineObliqueMercatorAzimuthCenter:
1227
68
        panProjParamId[0] = ProjCenterLatGeoKey;
1228
68
        panProjParamId[1] = ProjCenterLongGeoKey;
1229
68
        panProjParamId[2] = ProjAzimuthAngleGeoKey;
1230
68
        panProjParamId[3] = ProjRectifiedGridAngleGeoKey;
1231
68
        panProjParamId[4] = ProjScaleAtCenterGeoKey;
1232
68
        panProjParamId[5] = ProjFalseEastingGeoKey;
1233
68
        panProjParamId[6] = ProjFalseNorthingGeoKey;
1234
1235
68
        panEPSGCodes[0] = EPSGProjCenterLat;
1236
68
        panEPSGCodes[1] = EPSGProjCenterLong;
1237
68
        panEPSGCodes[2] = EPSGAzimuth;
1238
68
        panEPSGCodes[3] = EPSGAngleRectifiedToSkewedGrid;
1239
68
        panEPSGCodes[4] = EPSGInitialLineScaleFactor;
1240
68
        panEPSGCodes[5] = EPSGProjCenterEasting; /* EPSG proj method 9812 uses EPSGFalseEasting, but 9815 uses EPSGProjCenterEasting */
1241
68
        panEPSGCodes[6] = EPSGProjCenterNorthing; /* EPSG proj method 9812 uses EPSGFalseNorthing, but 9815 uses EPSGProjCenterNorthing */
1242
68
        return TRUE;
1243
1244
6
      case CT_ObliqueMercator_Laborde:
1245
6
        panProjParamId[0] = ProjCenterLatGeoKey;
1246
6
        panProjParamId[1] = ProjCenterLongGeoKey;
1247
6
        panProjParamId[2] = ProjAzimuthAngleGeoKey;
1248
6
        panProjParamId[4] = ProjScaleAtCenterGeoKey;
1249
6
        panProjParamId[5] = ProjFalseEastingGeoKey;
1250
6
        panProjParamId[6] = ProjFalseNorthingGeoKey;
1251
1252
6
        panEPSGCodes[0] = EPSGProjCenterLat;
1253
6
        panEPSGCodes[1] = EPSGProjCenterLong;
1254
6
        panEPSGCodes[2] = EPSGAzimuth;
1255
6
        panEPSGCodes[4] = EPSGInitialLineScaleFactor;
1256
6
        panEPSGCodes[5] = EPSGFalseEasting;
1257
6
        panEPSGCodes[6] = EPSGFalseNorthing;
1258
6
        return TRUE;
1259
1260
444
      case CT_LambertConfConic_1SP:
1261
552
      case CT_Mercator:
1262
600
      case CT_ObliqueStereographic:
1263
886
      case CT_PolarStereographic:
1264
16.8k
      case CT_TransverseMercator:
1265
16.9k
      case CT_TransvMercator_SouthOriented:
1266
16.9k
        panProjParamId[0] = ProjNatOriginLatGeoKey;
1267
16.9k
        if( nCTProjection == CT_PolarStereographic )
1268
286
        {
1269
286
            panProjParamId[1] = ProjStraightVertPoleLongGeoKey;
1270
286
        }
1271
16.6k
        else
1272
16.6k
        {
1273
16.6k
            panProjParamId[1] = ProjNatOriginLongGeoKey;
1274
16.6k
        }
1275
16.9k
        if( nEPSGProjMethod == 9805 ) /* Mercator_2SP */
1276
0
        {
1277
0
            panProjParamId[2] = ProjStdParallel1GeoKey;
1278
0
        }
1279
16.9k
        panProjParamId[4] = ProjScaleAtNatOriginGeoKey;
1280
16.9k
        panProjParamId[5] = ProjFalseEastingGeoKey;
1281
16.9k
        panProjParamId[6] = ProjFalseNorthingGeoKey;
1282
1283
16.9k
        panEPSGCodes[0] = EPSGNatOriginLat;
1284
16.9k
        panEPSGCodes[1] = EPSGNatOriginLong;
1285
16.9k
        if( nEPSGProjMethod == 9805 ) /* Mercator_2SP */
1286
0
        {
1287
0
            panEPSGCodes[2] = EPSGStdParallel1Lat;
1288
0
        }
1289
16.9k
        panEPSGCodes[4] = EPSGNatOriginScaleFactor;
1290
16.9k
        panEPSGCodes[5] = EPSGFalseEasting;
1291
16.9k
        panEPSGCodes[6] = EPSGFalseNorthing;
1292
16.9k
        return TRUE;
1293
1294
2.38k
      case CT_LambertConfConic_2SP:
1295
2.38k
        panProjParamId[0] = ProjFalseOriginLatGeoKey;
1296
2.38k
        panProjParamId[1] = ProjFalseOriginLongGeoKey;
1297
2.38k
        panProjParamId[2] = ProjStdParallel1GeoKey;
1298
2.38k
        panProjParamId[3] = ProjStdParallel2GeoKey;
1299
2.38k
        panProjParamId[5] = ProjFalseEastingGeoKey;
1300
2.38k
        panProjParamId[6] = ProjFalseNorthingGeoKey;
1301
1302
2.38k
        panEPSGCodes[0] = EPSGFalseOriginLat;
1303
2.38k
        panEPSGCodes[1] = EPSGFalseOriginLong;
1304
2.38k
        panEPSGCodes[2] = EPSGStdParallel1Lat;
1305
2.38k
        panEPSGCodes[3] = EPSGStdParallel2Lat;
1306
2.38k
        panEPSGCodes[5] = EPSGFalseOriginEasting;
1307
2.38k
        panEPSGCodes[6] = EPSGFalseOriginNorthing;
1308
2.38k
        return TRUE;
1309
1310
160
      case CT_AlbersEqualArea:
1311
160
        panProjParamId[0] = ProjStdParallel1GeoKey;
1312
160
        panProjParamId[1] = ProjStdParallel2GeoKey;
1313
160
        panProjParamId[2] = ProjNatOriginLatGeoKey;
1314
160
        panProjParamId[3] = ProjNatOriginLongGeoKey;
1315
160
        panProjParamId[5] = ProjFalseEastingGeoKey;
1316
160
        panProjParamId[6] = ProjFalseNorthingGeoKey;
1317
1318
160
        panEPSGCodes[0] = EPSGStdParallel1Lat;
1319
160
        panEPSGCodes[1] = EPSGStdParallel2Lat;
1320
160
        panEPSGCodes[2] = EPSGFalseOriginLat;
1321
160
        panEPSGCodes[3] = EPSGFalseOriginLong;
1322
160
        panEPSGCodes[5] = EPSGFalseOriginEasting;
1323
160
        panEPSGCodes[6] = EPSGFalseOriginNorthing;
1324
160
        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
78
      case CT_LambertAzimEqualArea:
1336
78
        panProjParamId[0] = ProjCenterLatGeoKey;
1337
78
        panProjParamId[1] = ProjCenterLongGeoKey;
1338
78
        panProjParamId[5] = ProjFalseEastingGeoKey;
1339
78
        panProjParamId[6] = ProjFalseNorthingGeoKey;
1340
1341
78
        panEPSGCodes[0] = EPSGNatOriginLat;
1342
78
        panEPSGCodes[1] = EPSGNatOriginLong;
1343
78
        panEPSGCodes[5] = EPSGFalseEasting;
1344
78
        panEPSGCodes[6] = EPSGFalseNorthing;
1345
78
        return TRUE;
1346
1347
6
      case CT_CylindricalEqualArea:
1348
6
        panProjParamId[0] = ProjStdParallel1GeoKey;
1349
6
        panProjParamId[1] = ProjNatOriginLongGeoKey;
1350
6
        panProjParamId[5] = ProjFalseEastingGeoKey;
1351
6
        panProjParamId[6] = ProjFalseNorthingGeoKey;
1352
1353
6
        panEPSGCodes[0] = EPSGStdParallel1Lat;
1354
6
        panEPSGCodes[1] = EPSGFalseOriginLong;
1355
6
        panEPSGCodes[5] = EPSGFalseOriginEasting;
1356
6
        panEPSGCodes[6] = EPSGFalseOriginNorthing;
1357
6
        return TRUE;
1358
1359
50
      case CT_Equirectangular:
1360
50
        panProjParamId[0] = ProjCenterLatGeoKey;
1361
50
        panProjParamId[1] = ProjCenterLongGeoKey;
1362
50
        panProjParamId[2] = ProjStdParallel1GeoKey;
1363
50
        panProjParamId[5] = ProjFalseEastingGeoKey;
1364
50
        panProjParamId[6] = ProjFalseNorthingGeoKey;
1365
1366
50
        panEPSGCodes[0] = EPSGNatOriginLat;
1367
50
        panEPSGCodes[1] = EPSGNatOriginLong;
1368
50
        panEPSGCodes[2] = EPSGStdParallel1Lat;
1369
50
        panEPSGCodes[5] = EPSGFalseEasting;
1370
50
        panEPSGCodes[6] = EPSGFalseNorthing;
1371
50
        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
5.83k
      default:
1388
5.83k
        return FALSE;
1389
25.5k
    }
1390
25.5k
}
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
19.2k
{
1408
19.2k
    if ((nProjTRFCode >= Proj_UTM_zone_1N && nProjTRFCode <= Proj_UTM_zone_60N) ||
1409
19.2k
        (nProjTRFCode >= Proj_UTM_zone_1S && nProjTRFCode <= Proj_UTM_zone_60S))
1410
7.76k
    {
1411
7.76k
        int bNorth;
1412
7.76k
        int nZone;
1413
7.76k
        if (nProjTRFCode <= Proj_UTM_zone_60N)
1414
4.22k
        {
1415
4.22k
            bNorth = TRUE;
1416
4.22k
            nZone = nProjTRFCode - Proj_UTM_zone_1N + 1;
1417
4.22k
        }
1418
3.53k
        else
1419
3.53k
        {
1420
3.53k
            bNorth = FALSE;
1421
3.53k
            nZone = nProjTRFCode - Proj_UTM_zone_1S + 1;
1422
3.53k
        }
1423
1424
7.76k
        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
7.76k
        if (pnProjMethod)
1433
7.76k
            *pnProjMethod = 9807;
1434
1435
7.76k
        if (padfProjParams)
1436
7.76k
        {
1437
7.76k
            padfProjParams[0] = 0;
1438
7.76k
            padfProjParams[1] = -183 + 6 * nZone;
1439
7.76k
            padfProjParams[2] = 0;
1440
7.76k
            padfProjParams[3] = 0;
1441
7.76k
            padfProjParams[4] = 0.9996;
1442
7.76k
            padfProjParams[5] = 500000;
1443
7.76k
            padfProjParams[6] = (bNorth) ? 0 : 10000000;
1444
7.76k
        }
1445
1446
7.76k
        return TRUE;
1447
7.76k
    }
1448
1449
11.5k
    if( nProjTRFCode == KvUserDefined )
1450
0
        return FALSE;
1451
1452
11.5k
    {
1453
11.5k
        char    szCode[12];
1454
11.5k
        snprintf(szCode, sizeof(szCode), "%d", nProjTRFCode);
1455
11.5k
        PJ_CONTEXT* ctx = (PJ_CONTEXT*)ctxIn;
1456
11.5k
        PJ *transf = proj_create_from_database(
1457
11.5k
            ctx, "EPSG", szCode, PJ_CATEGORY_COORDINATE_OPERATION, 0, NULL);
1458
11.5k
        if( !transf )
1459
4.04k
        {
1460
4.04k
            return FALSE;
1461
4.04k
        }
1462
1463
7.49k
        if( proj_get_type(transf) != PJ_TYPE_CONVERSION )
1464
1.23k
        {
1465
1.23k
            proj_destroy(transf);
1466
1.23k
            return FALSE;
1467
1.23k
        }
1468
1469
        /* Get the projection method code */
1470
6.25k
        const char* pszMethodCode = NULL;
1471
6.25k
        proj_coordoperation_get_method_info(ctx, transf,
1472
6.25k
                                            NULL, /* method name */
1473
6.25k
                                            NULL, /* method auth name (should be EPSG) */
1474
6.25k
                                            &pszMethodCode);
1475
6.25k
        assert( pszMethodCode );
1476
6.25k
        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
6.25k
        const int nCTProjMethod = EPSGProjMethodToCTProjMethod( nProjMethod, TRUE );
1483
6.25k
        int anEPSGCodes[7];
1484
6.25k
        SetGTParamIds( nCTProjMethod, nProjMethod, NULL, anEPSGCodes );
1485
1486
1487
/* -------------------------------------------------------------------- */
1488
/*      Get the parameters for this projection.                         */
1489
/* -------------------------------------------------------------------- */
1490
1491
6.25k
        double  adfProjParams[7];
1492
1493
50.0k
        for( int i = 0; i < 7; i++ )
1494
43.7k
        {
1495
43.7k
            int nEPSGCode = anEPSGCodes[i];
1496
1497
            /* Establish default */
1498
43.7k
            if( nEPSGCode == EPSGAngleRectifiedToSkewedGrid )
1499
34
                adfProjParams[i] = 90.0;
1500
43.7k
            else if( nEPSGCode == EPSGNatOriginScaleFactor
1501
43.7k
                    || nEPSGCode == EPSGInitialLineScaleFactor
1502
43.7k
                    || nEPSGCode == EPSGPseudoStdParallelScaleFactor )
1503
4.60k
                adfProjParams[i] = 1.0;
1504
39.1k
            else
1505
39.1k
                adfProjParams[i] = 0.0;
1506
1507
            /* If there is no parameter, skip */
1508
43.7k
            if( nEPSGCode == 0 )
1509
12.6k
                continue;
1510
1511
31.1k
            const int nParamCount = proj_coordoperation_get_param_count(ctx, transf);
1512
1513
            /* Find the matching parameter */
1514
31.1k
            const char *pszUOMCategory = NULL;
1515
31.1k
            double  dfValue = 0.0;
1516
31.1k
            double  dfUnitConvFactor = 0.0;
1517
31.1k
            int iEPSG;  /* Used after for */
1518
97.7k
            for( iEPSG = 0; iEPSG < nParamCount; iEPSG++ )
1519
97.4k
            {
1520
97.4k
                const char* pszParamCode = NULL;
1521
97.4k
                proj_coordoperation_get_param(
1522
97.4k
                    ctx, transf, iEPSG,
1523
97.4k
                    NULL, /* name */
1524
97.4k
                    NULL, /* auth name */
1525
97.4k
                    &pszParamCode,
1526
97.4k
                    &dfValue,
1527
97.4k
                    NULL, /* value (string) */
1528
97.4k
                    &dfUnitConvFactor, /* unit conv factor */
1529
97.4k
                    NULL, /* unit name */
1530
97.4k
                    NULL, /* unit auth name */
1531
97.4k
                    NULL, /* unit code */
1532
97.4k
                    &pszUOMCategory /* unit category */);
1533
97.4k
                assert(pszParamCode);
1534
97.4k
                if( atoi(pszParamCode) == nEPSGCode )
1535
30.9k
                {
1536
30.9k
                    break;
1537
30.9k
                }
1538
97.4k
            }
1539
1540
            /* not found, accept the default */
1541
31.1k
            if( iEPSG == nParamCount )
1542
235
            {
1543
                /* for CT_ObliqueMercator try alternate parameter codes first */
1544
                /* because EPSG proj method 9812 uses EPSGFalseXXXXX, but 9815 uses EPSGProjCenterXXXXX */
1545
235
                if ( nCTProjMethod == CT_ObliqueMercator && nEPSGCode == EPSGProjCenterEasting )
1546
30
                    nEPSGCode = EPSGFalseEasting;
1547
205
                else if ( nCTProjMethod == CT_ObliqueMercator && nEPSGCode == EPSGProjCenterNorthing )
1548
30
                    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
175
                else if( nCTProjMethod == CT_PolarStereographic && nEPSGCode == EPSGNatOriginLat )
1553
32
                    nEPSGCode = EPSGLatOfStdParallel;
1554
143
                else if( nCTProjMethod == CT_PolarStereographic && nEPSGCode == EPSGNatOriginLong )
1555
32
                    nEPSGCode = EPSGOriginLong;
1556
111
                else
1557
111
                    continue;
1558
1559
486
                for( iEPSG = 0; iEPSG < nParamCount; iEPSG++ )
1560
486
                {
1561
486
                    const char* pszParamCode = NULL;
1562
486
                    proj_coordoperation_get_param(
1563
486
                        ctx, transf, iEPSG,
1564
486
                        NULL, /* name */
1565
486
                        NULL, /* auth name */
1566
486
                        &pszParamCode,
1567
486
                        &dfValue,
1568
486
                        NULL, /* value (string) */
1569
486
                        &dfUnitConvFactor, /* unit conv factor */
1570
486
                        NULL, /* unit name */
1571
486
                        NULL, /* unit auth name */
1572
486
                        NULL, /* unit code */
1573
486
                        &pszUOMCategory /* unit category */);
1574
486
                    assert(pszParamCode);
1575
486
                    if( atoi(pszParamCode) == nEPSGCode )
1576
124
                    {
1577
124
                        break;
1578
124
                    }
1579
486
                }
1580
1581
124
                if( iEPSG == nParamCount )
1582
0
                    continue;
1583
124
            }
1584
1585
31.0k
            assert(pszUOMCategory);
1586
1587
31.0k
            adfProjParams[i] = dfValue * dfUnitConvFactor;
1588
31.0k
            if( strcmp(pszUOMCategory, "angular") == 0.0 )
1589
14.5k
            {
1590
                /* Convert from radians to degrees */
1591
14.5k
                adfProjParams[i] *= 180 / M_PI;
1592
14.5k
            }
1593
31.0k
        }
1594
1595
/* -------------------------------------------------------------------- */
1596
/*      Get the name, if requested.                                     */
1597
/* -------------------------------------------------------------------- */
1598
6.25k
        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
6.25k
        if( pnProjMethod != NULL )
1614
6.25k
            *pnProjMethod = (short) nProjMethod;
1615
1616
6.25k
        if( padfProjParams != NULL )
1617
6.25k
        {
1618
50.0k
            for( int i = 0; i < 7; i++ )
1619
43.7k
                padfProjParams[i] = adfProjParams[i];
1620
6.25k
        }
1621
1622
6.25k
        proj_destroy(transf);
1623
1624
6.25k
        return TRUE;
1625
6.25k
    }
1626
6.25k
}
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
9.29k
{
1651
1652
/* -------------------------------------------------------------------- */
1653
/*      Get the false easting, and northing if available.               */
1654
/* -------------------------------------------------------------------- */
1655
9.29k
    double dfFalseEasting = 0.0;
1656
9.29k
    if( !GTIFKeyGetDOUBLE(psGTIF, ProjFalseEastingGeoKey, &dfFalseEasting, 0, 1)
1657
9.29k
        && !GTIFKeyGetDOUBLE(psGTIF, ProjCenterEastingGeoKey,
1658
7.20k
                       &dfFalseEasting, 0, 1)
1659
9.29k
        && !GTIFKeyGetDOUBLE(psGTIF, ProjFalseOriginEastingGeoKey,
1660
7.16k
                       &dfFalseEasting, 0, 1) )
1661
7.16k
        dfFalseEasting = 0.0;
1662
1663
9.29k
    double dfFalseNorthing = 0.0;
1664
9.29k
    if( !GTIFKeyGetDOUBLE(psGTIF, ProjFalseNorthingGeoKey, &dfFalseNorthing,0,1)
1665
9.29k
        && !GTIFKeyGetDOUBLE(psGTIF, ProjCenterNorthingGeoKey,
1666
7.08k
                       &dfFalseNorthing, 0, 1)
1667
9.29k
        && !GTIFKeyGetDOUBLE(psGTIF, ProjFalseOriginNorthingGeoKey,
1668
7.05k
                       &dfFalseNorthing, 0, 1) )
1669
7.00k
        dfFalseNorthing = 0.0;
1670
1671
9.29k
    double dfNatOriginLong = 0.0, dfNatOriginLat = 0.0, dfRectGridAngle = 0.0;
1672
9.29k
    double dfNatOriginScale = 1.0;
1673
9.29k
    double dfStdParallel1 = 0.0, dfStdParallel2 = 0.0, dfAzimuth = 0.0;
1674
1675
9.29k
    switch( psDefn->CTProjection )
1676
9.29k
    {
1677
/* -------------------------------------------------------------------- */
1678
276
      case CT_Stereographic:
1679
/* -------------------------------------------------------------------- */
1680
276
        if( GTIFKeyGetDOUBLE(psGTIF, ProjNatOriginLongGeoKey,
1681
276
                       &dfNatOriginLong, 0, 1 ) == 0
1682
276
            && GTIFKeyGetDOUBLE(psGTIF, ProjFalseOriginLongGeoKey,
1683
269
                          &dfNatOriginLong, 0, 1 ) == 0
1684
276
            && GTIFKeyGetDOUBLE(psGTIF, ProjCenterLongGeoKey,
1685
266
                          &dfNatOriginLong, 0, 1 ) == 0 )
1686
264
            dfNatOriginLong = 0.0;
1687
1688
276
        if( GTIFKeyGetDOUBLE(psGTIF, ProjNatOriginLatGeoKey,
1689
276
                       &dfNatOriginLat, 0, 1 ) == 0
1690
276
            && GTIFKeyGetDOUBLE(psGTIF, ProjFalseOriginLatGeoKey,
1691
272
                          &dfNatOriginLat, 0, 1 ) == 0
1692
276
            && GTIFKeyGetDOUBLE(psGTIF, ProjCenterLatGeoKey,
1693
257
                          &dfNatOriginLat, 0, 1 ) == 0 )
1694
256
            dfNatOriginLat = 0.0;
1695
1696
276
        if( GTIFKeyGetDOUBLE(psGTIF, ProjScaleAtNatOriginGeoKey,
1697
276
                       &dfNatOriginScale, 0, 1 ) == 0 )
1698
262
            dfNatOriginScale = 1.0;
1699
1700
        /* notdef: should transform to decimal degrees at this point */
1701
1702
276
        psDefn->ProjParm[0] = dfNatOriginLat;
1703
276
        psDefn->ProjParmId[0] = ProjCenterLatGeoKey;
1704
276
        psDefn->ProjParm[1] = dfNatOriginLong;
1705
276
        psDefn->ProjParmId[1] = ProjCenterLongGeoKey;
1706
276
        psDefn->ProjParm[4] = dfNatOriginScale;
1707
276
        psDefn->ProjParmId[4] = ProjScaleAtNatOriginGeoKey;
1708
276
        psDefn->ProjParm[5] = dfFalseEasting;
1709
276
        psDefn->ProjParmId[5] = ProjFalseEastingGeoKey;
1710
276
        psDefn->ProjParm[6] = dfFalseNorthing;
1711
276
        psDefn->ProjParmId[6] = ProjFalseNorthingGeoKey;
1712
1713
276
        psDefn->nParms = 7;
1714
276
        break;
1715
1716
/* -------------------------------------------------------------------- */
1717
1.04k
      case CT_Mercator:
1718
/* -------------------------------------------------------------------- */
1719
1.04k
        if( GTIFKeyGetDOUBLE(psGTIF, ProjNatOriginLongGeoKey,
1720
1.04k
                       &dfNatOriginLong, 0, 1 ) == 0
1721
1.04k
            && GTIFKeyGetDOUBLE(psGTIF, ProjFalseOriginLongGeoKey,
1722
809
                          &dfNatOriginLong, 0, 1 ) == 0
1723
1.04k
            && GTIFKeyGetDOUBLE(psGTIF, ProjCenterLongGeoKey,
1724
806
                          &dfNatOriginLong, 0, 1 ) == 0 )
1725
788
            dfNatOriginLong = 0.0;
1726
1727
1.04k
        if( GTIFKeyGetDOUBLE(psGTIF, ProjNatOriginLatGeoKey,
1728
1.04k
                       &dfNatOriginLat, 0, 1 ) == 0
1729
1.04k
            && GTIFKeyGetDOUBLE(psGTIF, ProjFalseOriginLatGeoKey,
1730
265
                          &dfNatOriginLat, 0, 1 ) == 0
1731
1.04k
            && GTIFKeyGetDOUBLE(psGTIF, ProjCenterLatGeoKey,
1732
155
                          &dfNatOriginLat, 0, 1 ) == 0 )
1733
153
            dfNatOriginLat = 0.0;
1734
1735
1736
1.04k
        const int bHaveSP1 = GTIFKeyGetDOUBLE(psGTIF, ProjStdParallel1GeoKey,
1737
1.04k
                              &dfStdParallel1, 0, 1 );
1738
1739
1.04k
        int bHaveNOS = GTIFKeyGetDOUBLE(psGTIF, ProjScaleAtNatOriginGeoKey,
1740
1.04k
                              &dfNatOriginScale, 0, 1 );
1741
1742
        /* Default scale only if dfStdParallel1 isn't defined either */
1743
1.04k
        if( !bHaveNOS && !bHaveSP1)
1744
209
        {
1745
209
            bHaveNOS = TRUE;
1746
209
            dfNatOriginScale = 1.0;
1747
209
        }
1748
1749
        /* notdef: should transform to decimal degrees at this point */
1750
1751
1.04k
        psDefn->ProjParm[0] = dfNatOriginLat;
1752
1.04k
        psDefn->ProjParmId[0] = ProjNatOriginLatGeoKey;
1753
1.04k
        psDefn->ProjParm[1] = dfNatOriginLong;
1754
1.04k
        psDefn->ProjParmId[1] = ProjNatOriginLongGeoKey;
1755
1.04k
        if( bHaveSP1 )
1756
650
        {
1757
650
            psDefn->ProjParm[2] = dfStdParallel1;
1758
650
            psDefn->ProjParmId[2] = ProjStdParallel1GeoKey;
1759
650
        }
1760
1.04k
        if( bHaveNOS )
1761
394
        {
1762
394
            psDefn->ProjParm[4] = dfNatOriginScale;
1763
394
            psDefn->ProjParmId[4] = ProjScaleAtNatOriginGeoKey;
1764
394
        }
1765
1.04k
        psDefn->ProjParm[5] = dfFalseEasting;
1766
1.04k
        psDefn->ProjParmId[5] = ProjFalseEastingGeoKey;
1767
1.04k
        psDefn->ProjParm[6] = dfFalseNorthing;
1768
1.04k
        psDefn->ProjParmId[6] = ProjFalseNorthingGeoKey;
1769
1770
1.04k
        psDefn->nParms = 7;
1771
1.04k
        break;
1772
1773
/* -------------------------------------------------------------------- */
1774
502
      case CT_LambertConfConic_1SP:
1775
527
      case CT_ObliqueStereographic:
1776
812
      case CT_TransverseMercator:
1777
842
      case CT_TransvMercator_SouthOriented:
1778
/* -------------------------------------------------------------------- */
1779
842
        if( GTIFKeyGetDOUBLE(psGTIF, ProjNatOriginLongGeoKey,
1780
842
                       &dfNatOriginLong, 0, 1 ) == 0
1781
842
            && GTIFKeyGetDOUBLE(psGTIF, ProjFalseOriginLongGeoKey,
1782
585
                          &dfNatOriginLong, 0, 1 ) == 0
1783
842
            && GTIFKeyGetDOUBLE(psGTIF, ProjCenterLongGeoKey,
1784
580
                          &dfNatOriginLong, 0, 1 ) == 0 )
1785
512
            dfNatOriginLong = 0.0;
1786
1787
842
        if( GTIFKeyGetDOUBLE(psGTIF, ProjNatOriginLatGeoKey,
1788
842
                       &dfNatOriginLat, 0, 1 ) == 0
1789
842
            && GTIFKeyGetDOUBLE(psGTIF, ProjFalseOriginLatGeoKey,
1790
612
                          &dfNatOriginLat, 0, 1 ) == 0
1791
842
            && GTIFKeyGetDOUBLE(psGTIF, ProjCenterLatGeoKey,
1792
609
                          &dfNatOriginLat, 0, 1 ) == 0 )
1793
604
            dfNatOriginLat = 0.0;
1794
1795
842
        if( GTIFKeyGetDOUBLE(psGTIF, ProjScaleAtNatOriginGeoKey,
1796
842
                       &dfNatOriginScale, 0, 1 ) == 0
1797
            /* See https://github.com/OSGeo/gdal/files/1665718/lasinfo.txt */
1798
842
            && GTIFKeyGetDOUBLE(psGTIF, ProjScaleAtCenterGeoKey,
1799
546
                       &dfNatOriginScale, 0, 1 ) == 0 )
1800
545
            dfNatOriginScale = 1.0;
1801
1802
        /* notdef: should transform to decimal degrees at this point */
1803
1804
842
        psDefn->ProjParm[0] = dfNatOriginLat;
1805
842
        psDefn->ProjParmId[0] = ProjNatOriginLatGeoKey;
1806
842
        psDefn->ProjParm[1] = dfNatOriginLong;
1807
842
        psDefn->ProjParmId[1] = ProjNatOriginLongGeoKey;
1808
842
        psDefn->ProjParm[4] = dfNatOriginScale;
1809
842
        psDefn->ProjParmId[4] = ProjScaleAtNatOriginGeoKey;
1810
842
        psDefn->ProjParm[5] = dfFalseEasting;
1811
842
        psDefn->ProjParmId[5] = ProjFalseEastingGeoKey;
1812
842
        psDefn->ProjParm[6] = dfFalseNorthing;
1813
842
        psDefn->ProjParmId[6] = ProjFalseNorthingGeoKey;
1814
1815
842
        psDefn->nParms = 7;
1816
842
        break;
1817
1818
/* -------------------------------------------------------------------- */
1819
814
      case CT_ObliqueMercator: /* hotine */
1820
814
      case CT_HotineObliqueMercatorAzimuthCenter:
1821
/* -------------------------------------------------------------------- */
1822
814
        if( GTIFKeyGetDOUBLE(psGTIF, ProjNatOriginLongGeoKey,
1823
814
                       &dfNatOriginLong, 0, 1 ) == 0
1824
814
            && GTIFKeyGetDOUBLE(psGTIF, ProjFalseOriginLongGeoKey,
1825
809
                          &dfNatOriginLong, 0, 1 ) == 0
1826
814
            && GTIFKeyGetDOUBLE(psGTIF, ProjCenterLongGeoKey,
1827
808
                          &dfNatOriginLong, 0, 1 ) == 0 )
1828
807
            dfNatOriginLong = 0.0;
1829
1830
814
        if( GTIFKeyGetDOUBLE(psGTIF, ProjNatOriginLatGeoKey,
1831
814
                       &dfNatOriginLat, 0, 1 ) == 0
1832
814
            && GTIFKeyGetDOUBLE(psGTIF, ProjFalseOriginLatGeoKey,
1833
756
                          &dfNatOriginLat, 0, 1 ) == 0
1834
814
            && GTIFKeyGetDOUBLE(psGTIF, ProjCenterLatGeoKey,
1835
744
                          &dfNatOriginLat, 0, 1 ) == 0 )
1836
739
            dfNatOriginLat = 0.0;
1837
1838
814
        if( GTIFKeyGetDOUBLE(psGTIF, ProjAzimuthAngleGeoKey,
1839
814
                       &dfAzimuth, 0, 1 ) == 0 )
1840
811
            dfAzimuth = 0.0;
1841
1842
814
        if( GTIFKeyGetDOUBLE(psGTIF, ProjRectifiedGridAngleGeoKey,
1843
814
                       &dfRectGridAngle, 0, 1 ) == 0 )
1844
812
            dfRectGridAngle = 90.0;
1845
1846
814
        if( GTIFKeyGetDOUBLE(psGTIF, ProjScaleAtNatOriginGeoKey,
1847
814
                       &dfNatOriginScale, 0, 1 ) == 0
1848
814
            && GTIFKeyGetDOUBLE(psGTIF, ProjScaleAtCenterGeoKey,
1849
811
                          &dfNatOriginScale, 0, 1 ) == 0 )
1850
809
            dfNatOriginScale = 1.0;
1851
1852
        /* notdef: should transform to decimal degrees at this point */
1853
1854
814
        psDefn->ProjParm[0] = dfNatOriginLat;
1855
814
        psDefn->ProjParmId[0] = ProjCenterLatGeoKey;
1856
814
        psDefn->ProjParm[1] = dfNatOriginLong;
1857
814
        psDefn->ProjParmId[1] = ProjCenterLongGeoKey;
1858
814
        psDefn->ProjParm[2] = dfAzimuth;
1859
814
        psDefn->ProjParmId[2] = ProjAzimuthAngleGeoKey;
1860
814
        psDefn->ProjParm[3] = dfRectGridAngle;
1861
814
        psDefn->ProjParmId[3] = ProjRectifiedGridAngleGeoKey;
1862
814
        psDefn->ProjParm[4] = dfNatOriginScale;
1863
814
        psDefn->ProjParmId[4] = ProjScaleAtCenterGeoKey;
1864
814
        psDefn->ProjParm[5] = dfFalseEasting;
1865
814
        psDefn->ProjParmId[5] = ProjFalseEastingGeoKey;
1866
814
        psDefn->ProjParm[6] = dfFalseNorthing;
1867
814
        psDefn->ProjParmId[6] = ProjFalseNorthingGeoKey;
1868
1869
814
        psDefn->nParms = 7;
1870
814
        break;
1871
1872
/* -------------------------------------------------------------------- */
1873
875
      case CT_ObliqueMercator_Laborde:
1874
/* -------------------------------------------------------------------- */
1875
875
        if( GTIFKeyGetDOUBLE(psGTIF, ProjNatOriginLongGeoKey,
1876
875
                       &dfNatOriginLong, 0, 1 ) == 0
1877
875
            && GTIFKeyGetDOUBLE(psGTIF, ProjFalseOriginLongGeoKey,
1878
866
                          &dfNatOriginLong, 0, 1 ) == 0
1879
875
            && GTIFKeyGetDOUBLE(psGTIF, ProjCenterLongGeoKey,
1880
865
                          &dfNatOriginLong, 0, 1 ) == 0 )
1881
864
            dfNatOriginLong = 0.0;
1882
1883
875
        if( GTIFKeyGetDOUBLE(psGTIF, ProjNatOriginLatGeoKey,
1884
875
                       &dfNatOriginLat, 0, 1 ) == 0
1885
875
            && GTIFKeyGetDOUBLE(psGTIF, ProjFalseOriginLatGeoKey,
1886
853
                          &dfNatOriginLat, 0, 1 ) == 0
1887
875
            && GTIFKeyGetDOUBLE(psGTIF, ProjCenterLatGeoKey,
1888
851
                          &dfNatOriginLat, 0, 1 ) == 0 )
1889
848
            dfNatOriginLat = 0.0;
1890
1891
875
        if( GTIFKeyGetDOUBLE(psGTIF, ProjAzimuthAngleGeoKey,
1892
875
                       &dfAzimuth, 0, 1 ) == 0 )
1893
875
            dfAzimuth = 0.0;
1894
1895
875
        if( GTIFKeyGetDOUBLE(psGTIF, ProjScaleAtNatOriginGeoKey,
1896
875
                       &dfNatOriginScale, 0, 1 ) == 0
1897
875
            && GTIFKeyGetDOUBLE(psGTIF, ProjScaleAtCenterGeoKey,
1898
872
                          &dfNatOriginScale, 0, 1 ) == 0 )
1899
871
            dfNatOriginScale = 1.0;
1900
1901
        /* notdef: should transform to decimal degrees at this point */
1902
1903
875
        psDefn->ProjParm[0] = dfNatOriginLat;
1904
875
        psDefn->ProjParmId[0] = ProjCenterLatGeoKey;
1905
875
        psDefn->ProjParm[1] = dfNatOriginLong;
1906
875
        psDefn->ProjParmId[1] = ProjCenterLongGeoKey;
1907
875
        psDefn->ProjParm[2] = dfAzimuth;
1908
875
        psDefn->ProjParmId[2] = ProjAzimuthAngleGeoKey;
1909
875
        psDefn->ProjParm[4] = dfNatOriginScale;
1910
875
        psDefn->ProjParmId[4] = ProjScaleAtCenterGeoKey;
1911
875
        psDefn->ProjParm[5] = dfFalseEasting;
1912
875
        psDefn->ProjParmId[5] = ProjFalseEastingGeoKey;
1913
875
        psDefn->ProjParm[6] = dfFalseNorthing;
1914
875
        psDefn->ProjParmId[6] = ProjFalseNorthingGeoKey;
1915
1916
875
        psDefn->nParms = 7;
1917
875
        break;
1918
1919
/* -------------------------------------------------------------------- */
1920
121
      case CT_CassiniSoldner:
1921
140
      case CT_Polyconic:
1922
/* -------------------------------------------------------------------- */
1923
140
        if( GTIFKeyGetDOUBLE(psGTIF, ProjNatOriginLongGeoKey,
1924
140
                       &dfNatOriginLong, 0, 1 ) == 0
1925
140
            && GTIFKeyGetDOUBLE(psGTIF, ProjFalseOriginLongGeoKey,
1926
138
                          &dfNatOriginLong, 0, 1 ) == 0
1927
140
            && GTIFKeyGetDOUBLE(psGTIF, ProjCenterLongGeoKey,
1928
137
                          &dfNatOriginLong, 0, 1 ) == 0 )
1929
111
            dfNatOriginLong = 0.0;
1930
1931
140
        if( GTIFKeyGetDOUBLE(psGTIF, ProjNatOriginLatGeoKey,
1932
140
                       &dfNatOriginLat, 0, 1 ) == 0
1933
140
            && GTIFKeyGetDOUBLE(psGTIF, ProjFalseOriginLatGeoKey,
1934
138
                          &dfNatOriginLat, 0, 1 ) == 0
1935
140
            && GTIFKeyGetDOUBLE(psGTIF, ProjCenterLatGeoKey,
1936
126
                          &dfNatOriginLat, 0, 1 ) == 0 )
1937
124
            dfNatOriginLat = 0.0;
1938
1939
140
        if( GTIFKeyGetDOUBLE(psGTIF, ProjScaleAtNatOriginGeoKey,
1940
140
                       &dfNatOriginScale, 0, 1 ) == 0
1941
140
            && GTIFKeyGetDOUBLE(psGTIF, ProjScaleAtCenterGeoKey,
1942
138
                          &dfNatOriginScale, 0, 1 ) == 0 )
1943
138
            dfNatOriginScale = 1.0;
1944
1945
        /* notdef: should transform to decimal degrees at this point */
1946
1947
140
        psDefn->ProjParm[0] = dfNatOriginLat;
1948
140
        psDefn->ProjParmId[0] = ProjNatOriginLatGeoKey;
1949
140
        psDefn->ProjParm[1] = dfNatOriginLong;
1950
140
        psDefn->ProjParmId[1] = ProjNatOriginLongGeoKey;
1951
140
        psDefn->ProjParm[4] = dfNatOriginScale;
1952
140
        psDefn->ProjParmId[4] = ProjScaleAtNatOriginGeoKey;
1953
140
        psDefn->ProjParm[5] = dfFalseEasting;
1954
140
        psDefn->ProjParmId[5] = ProjFalseEastingGeoKey;
1955
140
        psDefn->ProjParm[6] = dfFalseNorthing;
1956
140
        psDefn->ProjParmId[6] = ProjFalseNorthingGeoKey;
1957
1958
140
        psDefn->nParms = 7;
1959
140
        break;
1960
1961
/* -------------------------------------------------------------------- */
1962
1.97k
      case CT_AzimuthalEquidistant:
1963
2.00k
      case CT_MillerCylindrical:
1964
2.27k
      case CT_Gnomonic:
1965
2.84k
      case CT_LambertAzimEqualArea:
1966
2.85k
      case CT_Orthographic:
1967
2.86k
      case CT_NewZealandMapGrid:
1968
/* -------------------------------------------------------------------- */
1969
2.86k
        if( GTIFKeyGetDOUBLE(psGTIF, ProjNatOriginLongGeoKey,
1970
2.86k
                       &dfNatOriginLong, 0, 1 ) == 0
1971
2.86k
            && GTIFKeyGetDOUBLE(psGTIF, ProjFalseOriginLongGeoKey,
1972
2.82k
                          &dfNatOriginLong, 0, 1 ) == 0
1973
2.86k
            && GTIFKeyGetDOUBLE(psGTIF, ProjCenterLongGeoKey,
1974
2.82k
                          &dfNatOriginLong, 0, 1 ) == 0 )
1975
2.80k
            dfNatOriginLong = 0.0;
1976
1977
2.86k
        if( GTIFKeyGetDOUBLE(psGTIF, ProjNatOriginLatGeoKey,
1978
2.86k
                       &dfNatOriginLat, 0, 1 ) == 0
1979
2.86k
            && GTIFKeyGetDOUBLE(psGTIF, ProjFalseOriginLatGeoKey,
1980
2.78k
                          &dfNatOriginLat, 0, 1 ) == 0
1981
2.86k
            && GTIFKeyGetDOUBLE(psGTIF, ProjCenterLatGeoKey,
1982
2.75k
                          &dfNatOriginLat, 0, 1 ) == 0 )
1983
2.75k
            dfNatOriginLat = 0.0;
1984
1985
        /* notdef: should transform to decimal degrees at this point */
1986
1987
2.86k
        psDefn->ProjParm[0] = dfNatOriginLat;
1988
2.86k
        psDefn->ProjParmId[0] = ProjCenterLatGeoKey;
1989
2.86k
        psDefn->ProjParm[1] = dfNatOriginLong;
1990
2.86k
        psDefn->ProjParmId[1] = ProjCenterLongGeoKey;
1991
2.86k
        psDefn->ProjParm[5] = dfFalseEasting;
1992
2.86k
        psDefn->ProjParmId[5] = ProjFalseEastingGeoKey;
1993
2.86k
        psDefn->ProjParm[6] = dfFalseNorthing;
1994
2.86k
        psDefn->ProjParmId[6] = ProjFalseNorthingGeoKey;
1995
1996
2.86k
        psDefn->nParms = 7;
1997
2.86k
        break;
1998
1999
/* -------------------------------------------------------------------- */
2000
737
      case CT_Equirectangular:
2001
/* -------------------------------------------------------------------- */
2002
737
        if( GTIFKeyGetDOUBLE(psGTIF, ProjNatOriginLongGeoKey,
2003
737
                       &dfNatOriginLong, 0, 1 ) == 0
2004
737
            && GTIFKeyGetDOUBLE(psGTIF, ProjFalseOriginLongGeoKey,
2005
615
                          &dfNatOriginLong, 0, 1 ) == 0
2006
737
            && GTIFKeyGetDOUBLE(psGTIF, ProjCenterLongGeoKey,
2007
591
                          &dfNatOriginLong, 0, 1 ) == 0 )
2008
588
            dfNatOriginLong = 0.0;
2009
2010
737
        if( GTIFKeyGetDOUBLE(psGTIF, ProjNatOriginLatGeoKey,
2011
737
                       &dfNatOriginLat, 0, 1 ) == 0
2012
737
            && GTIFKeyGetDOUBLE(psGTIF, ProjFalseOriginLatGeoKey,
2013
533
                          &dfNatOriginLat, 0, 1 ) == 0
2014
737
            && GTIFKeyGetDOUBLE(psGTIF, ProjCenterLatGeoKey,
2015
404
                          &dfNatOriginLat, 0, 1 ) == 0 )
2016
201
            dfNatOriginLat = 0.0;
2017
2018
737
        if( GTIFKeyGetDOUBLE(psGTIF, ProjStdParallel1GeoKey,
2019
737
                       &dfStdParallel1, 0, 1 ) == 0 )
2020
733
            dfStdParallel1 = 0.0;
2021
2022
        /* notdef: should transform to decimal degrees at this point */
2023
2024
737
        psDefn->ProjParm[0] = dfNatOriginLat;
2025
737
        psDefn->ProjParmId[0] = ProjCenterLatGeoKey;
2026
737
        psDefn->ProjParm[1] = dfNatOriginLong;
2027
737
        psDefn->ProjParmId[1] = ProjCenterLongGeoKey;
2028
737
        psDefn->ProjParm[2] = dfStdParallel1;
2029
737
        psDefn->ProjParmId[2] = ProjStdParallel1GeoKey;
2030
737
        psDefn->ProjParm[5] = dfFalseEasting;
2031
737
        psDefn->ProjParmId[5] = ProjFalseEastingGeoKey;
2032
737
        psDefn->ProjParm[6] = dfFalseNorthing;
2033
737
        psDefn->ProjParmId[6] = ProjFalseNorthingGeoKey;
2034
2035
737
        psDefn->nParms = 7;
2036
737
        break;
2037
2038
/* -------------------------------------------------------------------- */
2039
26
      case CT_Robinson:
2040
150
      case CT_Sinusoidal:
2041
162
      case CT_VanDerGrinten:
2042
/* -------------------------------------------------------------------- */
2043
162
        if( GTIFKeyGetDOUBLE(psGTIF, ProjNatOriginLongGeoKey,
2044
162
                       &dfNatOriginLong, 0, 1 ) == 0
2045
162
            && GTIFKeyGetDOUBLE(psGTIF, ProjFalseOriginLongGeoKey,
2046
152
                          &dfNatOriginLong, 0, 1 ) == 0
2047
162
            && GTIFKeyGetDOUBLE(psGTIF, ProjCenterLongGeoKey,
2048
151
                          &dfNatOriginLong, 0, 1 ) == 0 )
2049
146
            dfNatOriginLong = 0.0;
2050
2051
        /* notdef: should transform to decimal degrees at this point */
2052
2053
162
        psDefn->ProjParm[1] = dfNatOriginLong;
2054
162
        psDefn->ProjParmId[1] = ProjCenterLongGeoKey;
2055
162
        psDefn->ProjParm[5] = dfFalseEasting;
2056
162
        psDefn->ProjParmId[5] = ProjFalseEastingGeoKey;
2057
162
        psDefn->ProjParm[6] = dfFalseNorthing;
2058
162
        psDefn->ProjParmId[6] = ProjFalseNorthingGeoKey;
2059
2060
162
        psDefn->nParms = 7;
2061
162
        break;
2062
2063
/* -------------------------------------------------------------------- */
2064
219
      case CT_PolarStereographic:
2065
/* -------------------------------------------------------------------- */
2066
219
        if( GTIFKeyGetDOUBLE(psGTIF, ProjStraightVertPoleLongGeoKey,
2067
219
                       &dfNatOriginLong, 0, 1 ) == 0
2068
219
            && GTIFKeyGetDOUBLE(psGTIF, ProjNatOriginLongGeoKey,
2069
131
                          &dfNatOriginLong, 0, 1 ) == 0
2070
219
            && GTIFKeyGetDOUBLE(psGTIF, ProjFalseOriginLongGeoKey,
2071
103
                          &dfNatOriginLong, 0, 1 ) == 0
2072
219
            && GTIFKeyGetDOUBLE(psGTIF, ProjCenterLongGeoKey,
2073
99
                          &dfNatOriginLong, 0, 1 ) == 0 )
2074
81
            dfNatOriginLong = 0.0;
2075
2076
219
        if( GTIFKeyGetDOUBLE(psGTIF, ProjNatOriginLatGeoKey,
2077
219
                       &dfNatOriginLat, 0, 1 ) == 0
2078
219
            && GTIFKeyGetDOUBLE(psGTIF, ProjFalseOriginLatGeoKey,
2079
77
                          &dfNatOriginLat, 0, 1 ) == 0
2080
219
            && GTIFKeyGetDOUBLE(psGTIF, ProjCenterLatGeoKey,
2081
56
                          &dfNatOriginLat, 0, 1 ) == 0 )
2082
55
            dfNatOriginLat = 0.0;
2083
2084
219
        if( GTIFKeyGetDOUBLE(psGTIF, ProjScaleAtNatOriginGeoKey,
2085
219
                       &dfNatOriginScale, 0, 1 ) == 0
2086
219
            && GTIFKeyGetDOUBLE(psGTIF, ProjScaleAtCenterGeoKey,
2087
108
                          &dfNatOriginScale, 0, 1 ) == 0 )
2088
104
            dfNatOriginScale = 1.0;
2089
2090
        /* notdef: should transform to decimal degrees at this point */
2091
2092
219
        psDefn->ProjParm[0] = dfNatOriginLat;
2093
219
        psDefn->ProjParmId[0] = ProjNatOriginLatGeoKey;;
2094
219
        psDefn->ProjParm[1] = dfNatOriginLong;
2095
219
        psDefn->ProjParmId[1] = ProjStraightVertPoleLongGeoKey;
2096
219
        psDefn->ProjParm[4] = dfNatOriginScale;
2097
219
        psDefn->ProjParmId[4] = ProjScaleAtNatOriginGeoKey;
2098
219
        psDefn->ProjParm[5] = dfFalseEasting;
2099
219
        psDefn->ProjParmId[5] = ProjFalseEastingGeoKey;
2100
219
        psDefn->ProjParm[6] = dfFalseNorthing;
2101
219
        psDefn->ProjParmId[6] = ProjFalseNorthingGeoKey;
2102
2103
219
        psDefn->nParms = 7;
2104
219
        break;
2105
2106
/* -------------------------------------------------------------------- */
2107
42
      case CT_LambertConfConic_2SP:
2108
/* -------------------------------------------------------------------- */
2109
42
        if( GTIFKeyGetDOUBLE(psGTIF, ProjStdParallel1GeoKey,
2110
42
                       &dfStdParallel1, 0, 1 ) == 0 )
2111
40
            dfStdParallel1 = 0.0;
2112
2113
42
        if( GTIFKeyGetDOUBLE(psGTIF, ProjStdParallel2GeoKey,
2114
42
                       &dfStdParallel2, 0, 1 ) == 0 )
2115
37
            dfStdParallel1 = 0.0;
2116
2117
42
        if( GTIFKeyGetDOUBLE(psGTIF, ProjNatOriginLongGeoKey,
2118
42
                       &dfNatOriginLong, 0, 1 ) == 0
2119
42
            && GTIFKeyGetDOUBLE(psGTIF, ProjFalseOriginLongGeoKey,
2120
38
                          &dfNatOriginLong, 0, 1 ) == 0
2121
42
            && GTIFKeyGetDOUBLE(psGTIF, ProjCenterLongGeoKey,
2122
36
                          &dfNatOriginLong, 0, 1 ) == 0 )
2123
33
            dfNatOriginLong = 0.0;
2124
2125
42
        if( GTIFKeyGetDOUBLE(psGTIF, ProjNatOriginLatGeoKey,
2126
42
                       &dfNatOriginLat, 0, 1 ) == 0
2127
42
            && GTIFKeyGetDOUBLE(psGTIF, ProjFalseOriginLatGeoKey,
2128
39
                          &dfNatOriginLat, 0, 1 ) == 0
2129
42
            && GTIFKeyGetDOUBLE(psGTIF, ProjCenterLatGeoKey,
2130
31
                          &dfNatOriginLat, 0, 1 ) == 0 )
2131
30
            dfNatOriginLat = 0.0;
2132
2133
        /* notdef: should transform to decimal degrees at this point */
2134
2135
42
        psDefn->ProjParm[0] = dfNatOriginLat;
2136
42
        psDefn->ProjParmId[0] = ProjFalseOriginLatGeoKey;
2137
42
        psDefn->ProjParm[1] = dfNatOriginLong;
2138
42
        psDefn->ProjParmId[1] = ProjFalseOriginLongGeoKey;
2139
42
        psDefn->ProjParm[2] = dfStdParallel1;
2140
42
        psDefn->ProjParmId[2] = ProjStdParallel1GeoKey;
2141
42
        psDefn->ProjParm[3] = dfStdParallel2;
2142
42
        psDefn->ProjParmId[3] = ProjStdParallel2GeoKey;
2143
42
        psDefn->ProjParm[5] = dfFalseEasting;
2144
42
        psDefn->ProjParmId[5] = ProjFalseEastingGeoKey;
2145
42
        psDefn->ProjParm[6] = dfFalseNorthing;
2146
42
        psDefn->ProjParmId[6] = ProjFalseNorthingGeoKey;
2147
2148
42
        psDefn->nParms = 7;
2149
42
        break;
2150
2151
/* -------------------------------------------------------------------- */
2152
55
      case CT_AlbersEqualArea:
2153
387
      case CT_EquidistantConic:
2154
/* -------------------------------------------------------------------- */
2155
387
        if( GTIFKeyGetDOUBLE(psGTIF, ProjStdParallel1GeoKey,
2156
387
                       &dfStdParallel1, 0, 1 ) == 0 )
2157
250
            dfStdParallel1 = 0.0;
2158
2159
387
        if( GTIFKeyGetDOUBLE(psGTIF, ProjStdParallel2GeoKey,
2160
387
                       &dfStdParallel2, 0, 1 ) == 0 )
2161
252
            dfStdParallel2 = 0.0;
2162
2163
387
        if( GTIFKeyGetDOUBLE(psGTIF, ProjNatOriginLongGeoKey,
2164
387
                       &dfNatOriginLong, 0, 1 ) == 0
2165
387
            && GTIFKeyGetDOUBLE(psGTIF, ProjFalseOriginLongGeoKey,
2166
253
                          &dfNatOriginLong, 0, 1 ) == 0
2167
387
            && GTIFKeyGetDOUBLE(psGTIF, ProjCenterLongGeoKey,
2168
251
                          &dfNatOriginLong, 0, 1 ) == 0 )
2169
167
            dfNatOriginLong = 0.0;
2170
2171
387
        if( GTIFKeyGetDOUBLE(psGTIF, ProjNatOriginLatGeoKey,
2172
387
                       &dfNatOriginLat, 0, 1 ) == 0
2173
387
            && GTIFKeyGetDOUBLE(psGTIF, ProjFalseOriginLatGeoKey,
2174
226
                          &dfNatOriginLat, 0, 1 ) == 0
2175
387
            && GTIFKeyGetDOUBLE(psGTIF, ProjCenterLatGeoKey,
2176
224
                          &dfNatOriginLat, 0, 1 ) == 0 )
2177
223
            dfNatOriginLat = 0.0;
2178
2179
        /* notdef: should transform to decimal degrees at this point */
2180
2181
387
        psDefn->ProjParm[0] = dfStdParallel1;
2182
387
        psDefn->ProjParmId[0] = ProjStdParallel1GeoKey;
2183
387
        psDefn->ProjParm[1] = dfStdParallel2;
2184
387
        psDefn->ProjParmId[1] = ProjStdParallel2GeoKey;
2185
387
        psDefn->ProjParm[2] = dfNatOriginLat;
2186
387
        psDefn->ProjParmId[2] = ProjNatOriginLatGeoKey;
2187
387
        psDefn->ProjParm[3] = dfNatOriginLong;
2188
387
        psDefn->ProjParmId[3] = ProjNatOriginLongGeoKey;
2189
387
        psDefn->ProjParm[5] = dfFalseEasting;
2190
387
        psDefn->ProjParmId[5] = ProjFalseEastingGeoKey;
2191
387
        psDefn->ProjParm[6] = dfFalseNorthing;
2192
387
        psDefn->ProjParmId[6] = ProjFalseNorthingGeoKey;
2193
2194
387
        psDefn->nParms = 7;
2195
387
        break;
2196
2197
/* -------------------------------------------------------------------- */
2198
96
      case CT_CylindricalEqualArea:
2199
/* -------------------------------------------------------------------- */
2200
96
        if( GTIFKeyGetDOUBLE(psGTIF, ProjStdParallel1GeoKey,
2201
96
                       &dfStdParallel1, 0, 1 ) == 0 )
2202
91
            dfStdParallel1 = 0.0;
2203
2204
96
        if( GTIFKeyGetDOUBLE(psGTIF, ProjNatOriginLongGeoKey,
2205
96
                       &dfNatOriginLong, 0, 1 ) == 0
2206
96
            && GTIFKeyGetDOUBLE(psGTIF, ProjFalseOriginLongGeoKey,
2207
95
                          &dfNatOriginLong, 0, 1 ) == 0
2208
96
            && GTIFKeyGetDOUBLE(psGTIF, ProjCenterLongGeoKey,
2209
92
                          &dfNatOriginLong, 0, 1 ) == 0 )
2210
92
            dfNatOriginLong = 0.0;
2211
2212
        /* notdef: should transform to decimal degrees at this point */
2213
2214
96
        psDefn->ProjParm[0] = dfStdParallel1;
2215
96
        psDefn->ProjParmId[0] = ProjStdParallel1GeoKey;
2216
96
        psDefn->ProjParm[1] = dfNatOriginLong;
2217
96
        psDefn->ProjParmId[1] = ProjNatOriginLongGeoKey;
2218
96
        psDefn->ProjParm[5] = dfFalseEasting;
2219
96
        psDefn->ProjParmId[5] = ProjFalseEastingGeoKey;
2220
96
        psDefn->ProjParm[6] = dfFalseNorthing;
2221
96
        psDefn->ProjParmId[6] = ProjFalseNorthingGeoKey;
2222
2223
96
        psDefn->nParms = 7;
2224
96
        break;
2225
9.29k
    }
2226
2227
72.7k
    for( int iParam = 0; iParam < psDefn->nParms; iParam++ )
2228
63.4k
    {
2229
63.4k
        switch( psDefn->ProjParmId[iParam] )
2230
63.4k
        {
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
8.83k
          case ProjFalseEastingGeoKey:
2238
17.6k
          case ProjFalseNorthingGeoKey:
2239
17.6k
          case ProjFalseOriginEastingGeoKey:
2240
17.6k
          case ProjFalseOriginNorthingGeoKey:
2241
17.6k
          case ProjCenterEastingGeoKey:
2242
17.6k
          case ProjCenterNorthingGeoKey:
2243
17.6k
            if( psDefn->UOMLengthInMeters != 0
2244
17.6k
                && psDefn->UOMLengthInMeters != 1.0 )
2245
3.00k
            {
2246
3.00k
                psDefn->ProjParm[iParam] *= psDefn->UOMLengthInMeters;
2247
3.00k
            }
2248
17.6k
            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
2.09k
          case ProjStdParallel1GeoKey:
2258
3.08k
          case ProjStdParallel2GeoKey:
2259
5.92k
          case ProjNatOriginLongGeoKey:
2260
8.90k
          case ProjNatOriginLatGeoKey:
2261
8.97k
          case ProjFalseOriginLongGeoKey:
2262
9.04k
          case ProjFalseOriginLatGeoKey:
2263
14.7k
          case ProjCenterLongGeoKey:
2264
20.3k
          case ProjCenterLatGeoKey:
2265
20.5k
          case ProjStraightVertPoleLongGeoKey:
2266
21.3k
          case ProjRectifiedGridAngleGeoKey:
2267
21.3k
            if( psDefn->UOMAngleInDegrees != 0
2268
21.3k
                && psDefn->UOMAngleInDegrees != 1.0 )
2269
12.1k
            {
2270
12.1k
                psDefn->ProjParm[iParam] *= psDefn->UOMAngleInDegrees;
2271
12.1k
            }
2272
21.3k
            break;
2273
2274
24.4k
          default:
2275
24.4k
            break;
2276
63.4k
        }
2277
63.4k
    }
2278
9.29k
}
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
114k
{
2397
114k
    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
114k
    psDefn->DefnSet = 1;
2406
114k
    psDefn->Model = KvUserDefined;
2407
114k
    psDefn->PCS = KvUserDefined;
2408
114k
    psDefn->GCS = KvUserDefined;
2409
114k
    psDefn->UOMLength = KvUserDefined;
2410
114k
    psDefn->UOMLengthInMeters = 1.0;
2411
114k
    psDefn->UOMAngle = KvUserDefined;
2412
114k
    psDefn->UOMAngleInDegrees = 1.0;
2413
114k
    psDefn->Datum = KvUserDefined;
2414
114k
    psDefn->Ellipsoid = KvUserDefined;
2415
114k
    psDefn->SemiMajor = 0.0;
2416
114k
    psDefn->SemiMinor = 0.0;
2417
114k
    psDefn->PM = KvUserDefined;
2418
114k
    psDefn->PMLongToGreenwich = 0.0;
2419
114k
#if !defined(GEO_NORMALIZE_DISABLE_TOWGS84)
2420
114k
    psDefn->TOWGS84Count = 0;
2421
114k
    memset( psDefn->TOWGS84, 0, sizeof(psDefn->TOWGS84) );
2422
114k
#endif
2423
2424
114k
    psDefn->ProjCode = KvUserDefined;
2425
114k
    psDefn->Projection = KvUserDefined;
2426
114k
    psDefn->CTProjection = KvUserDefined;
2427
2428
114k
    psDefn->nParms = 0;
2429
1.26M
    for( int i = 0; i < MAX_GTIF_PROJPARMS; i++ )
2430
1.14M
    {
2431
1.14M
        psDefn->ProjParm[i] = 0.0;
2432
1.14M
        psDefn->ProjParmId[i] = 0;
2433
1.14M
    }
2434
2435
114k
    psDefn->MapSys = KvUserDefined;
2436
114k
    psDefn->Zone = 0;
2437
2438
/* -------------------------------------------------------------------- */
2439
/*      Do we have any geokeys?                                         */
2440
/* -------------------------------------------------------------------- */
2441
114k
    {
2442
114k
        int     nKeyCount = 0;
2443
114k
        int     anVersion[3];
2444
114k
        GTIFDirectoryInfo( psGTIF, anVersion, &nKeyCount );
2445
2446
114k
        if( nKeyCount == 0 )
2447
73.8k
        {
2448
73.8k
            psDefn->DefnSet = 0;
2449
73.8k
            return FALSE;
2450
73.8k
        }
2451
114k
    }
2452
2453
/* -------------------------------------------------------------------- */
2454
/*  Try to get the overall model type.        */
2455
/* -------------------------------------------------------------------- */
2456
41.0k
    GTIFKeyGetSSHORT(psGTIF,GTModelTypeGeoKey,&(psDefn->Model));
2457
2458
/* -------------------------------------------------------------------- */
2459
/*  Extract the Geog units.           */
2460
/* -------------------------------------------------------------------- */
2461
41.0k
    short nGeogUOMLinear = 9001; /* Linear_Meter */
2462
41.0k
    if( GTIFKeyGetSSHORT(psGTIF, GeogLinearUnitsGeoKey, &nGeogUOMLinear) == 1 )
2463
102
    {
2464
102
        psDefn->UOMLength = nGeogUOMLinear;
2465
102
    }
2466
2467
/* -------------------------------------------------------------------- */
2468
/*      Try to get a PCS.                                               */
2469
/* -------------------------------------------------------------------- */
2470
41.0k
    if( GTIFKeyGetSSHORT(psGTIF,ProjectedCSTypeGeoKey, &(psDefn->PCS)) == 1
2471
41.0k
        && psDefn->PCS != KvUserDefined )
2472
22.4k
    {
2473
        /*
2474
         * Translate this into useful information.
2475
         */
2476
22.4k
        GTIFGetPCSInfoEx( psGTIF->pj_context,
2477
22.4k
                          psDefn->PCS, NULL, &(psDefn->ProjCode),
2478
22.4k
                          &(psDefn->UOMLength), &(psDefn->GCS) );
2479
22.4k
    }
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
41.0k
    if( psDefn->PCS != KvUserDefined && psDefn->ProjCode == KvUserDefined )
2487
9.01k
    {
2488
9.01k
        int nZone;
2489
9.01k
        int nGCS = psDefn->GCS;
2490
2491
9.01k
        const int nMapSys = GTIFPCSToMapSys( psDefn->PCS, &nGCS, &nZone );
2492
9.01k
        if( nMapSys != KvUserDefined )
2493
3.77k
        {
2494
3.77k
            psDefn->ProjCode = (short) GTIFMapSysToProj( nMapSys, nZone );
2495
3.77k
            psDefn->GCS = (short) nGCS;
2496
3.77k
        }
2497
9.01k
    }
2498
2499
/* -------------------------------------------------------------------- */
2500
/*      If the Proj_ code is specified directly, use that.              */
2501
/* -------------------------------------------------------------------- */
2502
41.0k
    if( psDefn->ProjCode == KvUserDefined )
2503
23.8k
        GTIFKeyGetSSHORT(psGTIF, ProjectionGeoKey, &(psDefn->ProjCode));
2504
2505
41.0k
    if( psDefn->ProjCode != KvUserDefined )
2506
19.2k
    {
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
19.2k
        GTIFGetProjTRFInfoEx( psGTIF->pj_context,
2515
19.2k
                              psDefn->ProjCode, NULL, &(psDefn->Projection),
2516
19.2k
                              psDefn->ProjParm );
2517
2518
        /*
2519
         * Set the GeoTIFF identity of the parameters.
2520
         */
2521
19.2k
        psDefn->CTProjection = (short)
2522
19.2k
            EPSGProjMethodToCTProjMethod( psDefn->Projection, FALSE );
2523
2524
19.2k
        SetGTParamIds( EPSGProjMethodToCTProjMethod(psDefn->Projection, TRUE),
2525
19.2k
                      psDefn->Projection,
2526
19.2k
                      psDefn->ProjParmId, NULL);
2527
19.2k
        psDefn->nParms = 7;
2528
19.2k
    }
2529
2530
/* -------------------------------------------------------------------- */
2531
/*      Try to get a GCS.  If found, it will override any implied by    */
2532
/*      the PCS.                                                        */
2533
/* -------------------------------------------------------------------- */
2534
41.0k
    GTIFKeyGetSSHORT(psGTIF, GeographicTypeGeoKey, &(psDefn->GCS));
2535
41.0k
    if( psDefn->GCS < 1 || psDefn->GCS >= KvUserDefined )
2536
21.7k
        psDefn->GCS = KvUserDefined;
2537
2538
/* -------------------------------------------------------------------- */
2539
/*      Derive the datum, and prime meridian from the GCS.              */
2540
/* -------------------------------------------------------------------- */
2541
41.0k
    if( psDefn->GCS != KvUserDefined )
2542
19.2k
    {
2543
19.2k
        GTIFGetGCSInfoEx( psGTIF->pj_context,
2544
19.2k
                          psDefn->GCS, NULL, &(psDefn->Datum), &(psDefn->PM),
2545
19.2k
                          &(psDefn->UOMAngle) );
2546
19.2k
    }
2547
2548
/* -------------------------------------------------------------------- */
2549
/*      Handle the GCS angular units.  GeogAngularUnitsGeoKey           */
2550
/*      overrides the GCS or PCS setting.                               */
2551
/* -------------------------------------------------------------------- */
2552
41.0k
    GTIFKeyGetSSHORT(psGTIF, GeogAngularUnitsGeoKey, &(psDefn->UOMAngle));
2553
41.0k
    if( psDefn->UOMAngle != KvUserDefined )
2554
34.9k
    {
2555
34.9k
        GTIFGetUOMAngleInfoEx( psGTIF->pj_context,
2556
34.9k
                               psDefn->UOMAngle, NULL,
2557
34.9k
                               &(psDefn->UOMAngleInDegrees) );
2558
34.9k
    }
2559
2560
/* -------------------------------------------------------------------- */
2561
/*      Check for a datum setting, and then use the datum to derive     */
2562
/*      an ellipsoid.                                                   */
2563
/* -------------------------------------------------------------------- */
2564
41.0k
    GTIFKeyGetSSHORT(psGTIF, GeogGeodeticDatumGeoKey, &(psDefn->Datum));
2565
2566
41.0k
    if( psDefn->Datum != KvUserDefined )
2567
20.0k
    {
2568
20.0k
        GTIFGetDatumInfoEx( psGTIF->pj_context,
2569
20.0k
                            psDefn->Datum, NULL, &(psDefn->Ellipsoid) );
2570
20.0k
    }
2571
2572
/* -------------------------------------------------------------------- */
2573
/*      Check for an explicit ellipsoid.  Use the ellipsoid to          */
2574
/*      derive the ellipsoid characteristics, if possible.              */
2575
/* -------------------------------------------------------------------- */
2576
41.0k
    GTIFKeyGetSSHORT(psGTIF, GeogEllipsoidGeoKey, &(psDefn->Ellipsoid));
2577
2578
41.0k
    if( psDefn->Ellipsoid != KvUserDefined )
2579
17.2k
    {
2580
17.2k
        GTIFGetEllipsoidInfoEx( psGTIF->pj_context,
2581
17.2k
                                psDefn->Ellipsoid, NULL,
2582
17.2k
                                &(psDefn->SemiMajor), &(psDefn->SemiMinor) );
2583
17.2k
    }
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
41.0k
    CPL_IGNORE_RET_VAL_INT(GTIFKeyGetDOUBLE(psGTIF, GeogSemiMajorAxisGeoKey, &(psDefn->SemiMajor), 0, 1 ));
2591
41.0k
    CPL_IGNORE_RET_VAL_INT(GTIFKeyGetDOUBLE(psGTIF, GeogSemiMinorAxisGeoKey, &(psDefn->SemiMinor), 0, 1 ));
2592
2593
41.0k
    double dfInvFlattening;
2594
41.0k
    if( GTIFKeyGetDOUBLE(psGTIF, GeogInvFlatteningGeoKey, &dfInvFlattening,
2595
41.0k
                   0, 1 ) == 1 )
2596
1.49k
    {
2597
1.49k
        if( dfInvFlattening != 0.0 )
2598
1.48k
            psDefn->SemiMinor =
2599
1.48k
                psDefn->SemiMajor * (1 - 1.0/dfInvFlattening);
2600
12
        else
2601
12
            psDefn->SemiMinor = psDefn->SemiMajor;
2602
1.49k
    }
2603
2604
/* -------------------------------------------------------------------- */
2605
/*      Get the prime meridian info.                                    */
2606
/* -------------------------------------------------------------------- */
2607
41.0k
    GTIFKeyGetSSHORT(psGTIF, GeogPrimeMeridianGeoKey, &(psDefn->PM));
2608
2609
41.0k
    if( psDefn->PM != KvUserDefined )
2610
17.5k
    {
2611
17.5k
        GTIFGetPMInfoEx( psGTIF->pj_context,
2612
17.5k
                         psDefn->PM, NULL, &(psDefn->PMLongToGreenwich) );
2613
17.5k
    }
2614
23.5k
    else
2615
23.5k
    {
2616
23.5k
        CPL_IGNORE_RET_VAL_INT(GTIFKeyGetDOUBLE(psGTIF, GeogPrimeMeridianLongGeoKey,
2617
23.5k
                   &(psDefn->PMLongToGreenwich), 0, 1 ));
2618
2619
23.5k
        psDefn->PMLongToGreenwich =
2620
23.5k
            GTIFAngleToDD( psDefn->PMLongToGreenwich,
2621
23.5k
                           psDefn->UOMAngle );
2622
23.5k
    }
2623
2624
/* -------------------------------------------------------------------- */
2625
/*      Get the TOWGS84 parameters.                                     */
2626
/* -------------------------------------------------------------------- */
2627
41.0k
#if !defined(GEO_NORMALIZE_DISABLE_TOWGS84)
2628
41.0k
    psDefn->TOWGS84Count =
2629
41.0k
        (short)GTIFKeyGetDOUBLE(psGTIF, GeogTOWGS84GeoKey, psDefn->TOWGS84, 0, 7 );
2630
41.0k
#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
41.0k
    GTIFKeyGetSSHORT(psGTIF,ProjLinearUnitsGeoKey,&(psDefn->UOMLength));
2639
2640
41.0k
    if( psDefn->UOMLength != KvUserDefined )
2641
25.6k
    {
2642
25.6k
        GTIFGetUOMLengthInfoEx( psGTIF->pj_context,
2643
25.6k
                                psDefn->UOMLength, NULL,
2644
25.6k
                                &(psDefn->UOMLengthInMeters) );
2645
25.6k
    }
2646
15.4k
    else
2647
15.4k
    {
2648
15.4k
        CPL_IGNORE_RET_VAL_INT(GTIFKeyGetDOUBLE(psGTIF,ProjLinearUnitSizeGeoKey,&(psDefn->UOMLengthInMeters),0,1));
2649
15.4k
    }
2650
2651
/* -------------------------------------------------------------------- */
2652
/*      Handle a variety of user defined transform types.               */
2653
/* -------------------------------------------------------------------- */
2654
41.0k
    if( GTIFKeyGetSSHORT(psGTIF,ProjCoordTransGeoKey,
2655
41.0k
                   &(psDefn->CTProjection)) == 1)
2656
9.29k
    {
2657
9.29k
        GTIFFetchProjParms( psGTIF, psDefn );
2658
9.29k
    }
2659
2660
/* -------------------------------------------------------------------- */
2661
/*      Try to set the zoned map system information.                    */
2662
/* -------------------------------------------------------------------- */
2663
41.0k
    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
41.0k
    if( (psDefn->MapSys == MapSys_UTM_North
2671
41.0k
         || psDefn->MapSys == MapSys_UTM_South)
2672
41.0k
        && psDefn->CTProjection == KvUserDefined )
2673
0
    {
2674
0
        psDefn->CTProjection = CT_TransverseMercator;
2675
0
        psDefn->nParms = 7;
2676
0
        psDefn->ProjParmId[0] = ProjNatOriginLatGeoKey;
2677
0
        psDefn->ProjParm[0] = 0.0;
2678
2679
0
        psDefn->ProjParmId[1] = ProjNatOriginLongGeoKey;
2680
0
        psDefn->ProjParm[1] = psDefn->Zone*6 - 183.0;
2681
2682
0
        psDefn->ProjParmId[4] = ProjScaleAtNatOriginGeoKey;
2683
0
        psDefn->ProjParm[4] = 0.9996;
2684
2685
0
        psDefn->ProjParmId[5] = ProjFalseEastingGeoKey;
2686
0
        psDefn->ProjParm[5] = 500000.0;
2687
2688
0
        psDefn->ProjParmId[6] = ProjFalseNorthingGeoKey;
2689
2690
0
        if( psDefn->MapSys == MapSys_UTM_North )
2691
0
            psDefn->ProjParm[6] = 0.0;
2692
0
        else
2693
0
            psDefn->ProjParm[6] = 10000000.0;
2694
0
    }
2695
2696
41.0k
    return TRUE;
2697
114k
}
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
142k
{
3000
142k
    if( pMemory != NULL )
3001
125k
        VSIFree( pMemory );
3002
142k
}
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
114k
{
3014
114k
    return (GTIFDefn *) CPLCalloc(sizeof(GTIFDefn),1);
3015
114k
}
3016
3017
/************************************************************************/
3018
/*                            GTIFFreeDefn()                            */
3019
/*                                                                      */
3020
/*      Free a GTIF structure allocated by GTIFAllocDefn().             */
3021
/************************************************************************/
3022
3023
void GTIFFreeDefn( GTIFDefn *defn )
3024
114k
{
3025
114k
    VSIFree( defn );
3026
114k
}
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
245k
{
3037
245k
    if( psGTIF->own_pj_context )
3038
0
    {
3039
0
        proj_context_destroy(psGTIF->pj_context);
3040
0
    }
3041
245k
    psGTIF->own_pj_context = FALSE;
3042
245k
    psGTIF->pj_context = (PJ_CONTEXT*) pjContext;
3043
245k
}
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
189k
{
3056
189k
    if( psGTIF->pj_context || !instantiateIfNeeded )
3057
189k
    {
3058
189k
        if( out_gtif_own_pj_context )
3059
0
        {
3060
0
            *out_gtif_own_pj_context = psGTIF->own_pj_context;
3061
0
        }
3062
189k
        return psGTIF->pj_context;
3063
189k
    }
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
189k
}
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