Coverage Report

Created: 2025-06-13 06:29

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