Coverage Report

Created: 2025-06-13 06:29

/src/gdal/frmts/gtiff/gt_wkt_srs.cpp
Line
Count
Source (jump to first uncovered line)
1
/******************************************************************************
2
 *
3
 * Project:  GeoTIFF Driver
4
 * Purpose:  Implements translation between GeoTIFF normalized projection
5
 *           definitions and OpenGIS WKT SRS format.
6
 * Author:   Frank Warmerdam, warmerdam@pobox.com
7
 *
8
 ******************************************************************************
9
 * Copyright (c) 1999, Frank Warmerdam
10
 * Copyright (c) 2008-2014, Even Rouault <even dot rouault at spatialys.com>
11
 *
12
 * SPDX-License-Identifier: MIT
13
 ****************************************************************************/
14
15
#include "cpl_port.h"
16
17
#include "gt_wkt_srs.h"
18
19
#include <cmath>
20
#include <cstdio>
21
#include <cstdlib>
22
#include <cstring>
23
24
#include <algorithm>
25
#include <mutex>
26
27
#include "cpl_conv.h"
28
#include "cpl_error.h"
29
#include "cpl_string.h"
30
#include "cpl_vsi.h"
31
#include "gt_citation.h"
32
#include "gt_wkt_srs_for_gdal.h"
33
#include "gt_wkt_srs_priv.h"
34
#include "gtiff.h"
35
#include "gdal.h"
36
#include "geokeys.h"
37
#include "geovalues.h"
38
#include "ogr_core.h"
39
#include "ogr_spatialref.h"
40
#include "ogr_srs_api.h"
41
#include "ogr_proj_p.h"
42
#include "tiff.h"
43
#include "tiffio.h"
44
#include "tifvsi.h"
45
#include "xtiffio.h"
46
47
#include "proj.h"
48
49
static const geokey_t ProjLinearUnitsInterpCorrectGeoKey =
50
    static_cast<geokey_t>(3059);
51
52
#ifndef CT_HotineObliqueMercatorAzimuthCenter
53
0
#define CT_HotineObliqueMercatorAzimuthCenter 9815
54
#endif
55
56
#if !defined(GTIFAtof)
57
0
#define GTIFAtof CPLAtof
58
#endif
59
60
// To remind myself not to use CPLString in this file!
61
#define CPLString Please_do_not_use_CPLString_in_this_file
62
63
static const char *const papszDatumEquiv[] = {
64
    "Militar_Geographische_Institut",
65
    "Militar_Geographische_Institute",
66
    "World_Geodetic_System_1984",
67
    "WGS_1984",
68
    "WGS_72_Transit_Broadcast_Ephemeris",
69
    "WGS_1972_Transit_Broadcast_Ephemeris",
70
    "World_Geodetic_System_1972",
71
    "WGS_1972",
72
    "European_Terrestrial_Reference_System_89",
73
    "European_Reference_System_1989",
74
    "D_North_American_1927",
75
    "North_American_Datum_1927",  // #6863
76
    nullptr};
77
78
// Older libgeotiff's won't list this.
79
#ifndef CT_CylindricalEqualArea
80
0
#define CT_CylindricalEqualArea 28
81
#endif
82
83
#if LIBGEOTIFF_VERSION < 1700
84
constexpr geokey_t CoordinateEpochGeoKey = static_cast<geokey_t>(5120);
85
#endif
86
87
// Exists since 8.0.1
88
#ifndef PROJ_AT_LEAST_VERSION
89
#define PROJ_COMPUTE_VERSION(maj, min, patch)                                  \
90
    ((maj)*10000 + (min)*100 + (patch))
91
#define PROJ_VERSION_NUMBER                                                    \
92
    PROJ_COMPUTE_VERSION(PROJ_VERSION_MAJOR, PROJ_VERSION_MINOR,               \
93
                         PROJ_VERSION_PATCH)
94
#define PROJ_AT_LEAST_VERSION(maj, min, patch)                                 \
95
    (PROJ_VERSION_NUMBER >= PROJ_COMPUTE_VERSION(maj, min, patch))
96
#endif
97
98
/************************************************************************/
99
/*                       LibgeotiffOneTimeInit()                        */
100
/************************************************************************/
101
102
static std::mutex oDeleteMutex;
103
104
void LibgeotiffOneTimeInit()
105
0
{
106
0
    std::lock_guard<std::mutex> oLock(oDeleteMutex);
107
108
0
    static bool bOneTimeInitDone = false;
109
110
0
    if (bOneTimeInitDone)
111
0
        return;
112
113
0
    bOneTimeInitDone = true;
114
115
    // This isn't thread-safe, so better do it now
116
0
    XTIFFInitialize();
117
0
}
118
119
/************************************************************************/
120
/*                       GTIFToCPLRecyleString()                        */
121
/*                                                                      */
122
/*      This changes a string from the libgeotiff heap to the GDAL      */
123
/*      heap.                                                           */
124
/************************************************************************/
125
126
static void GTIFToCPLRecycleString(char **ppszTarget)
127
128
0
{
129
0
    if (*ppszTarget == nullptr)
130
0
        return;
131
132
0
    char *pszTempString = CPLStrdup(*ppszTarget);
133
0
    GTIFFreeMemory(*ppszTarget);
134
0
    *ppszTarget = pszTempString;
135
0
}
136
137
/************************************************************************/
138
/*                          WKTMassageDatum()                           */
139
/*                                                                      */
140
/*      Massage an EPSG datum name into WMT format.  Also transform     */
141
/*      specific exception cases into WKT versions.                     */
142
/************************************************************************/
143
144
static void WKTMassageDatum(char **ppszDatum)
145
146
0
{
147
0
    char *pszDatum = *ppszDatum;
148
0
    if (!pszDatum || pszDatum[0] == '\0')
149
0
        return;
150
151
    /* -------------------------------------------------------------------- */
152
    /*      Translate non-alphanumeric values to underscores.               */
153
    /* -------------------------------------------------------------------- */
154
0
    for (int i = 0; pszDatum[i] != '\0'; i++)
155
0
    {
156
0
        if (pszDatum[i] != '+' && !(pszDatum[i] >= 'A' && pszDatum[i] <= 'Z') &&
157
0
            !(pszDatum[i] >= 'a' && pszDatum[i] <= 'z') &&
158
0
            !(pszDatum[i] >= '0' && pszDatum[i] <= '9'))
159
0
        {
160
0
            pszDatum[i] = '_';
161
0
        }
162
0
    }
163
164
    /* -------------------------------------------------------------------- */
165
    /*      Remove repeated and trailing underscores.                       */
166
    /* -------------------------------------------------------------------- */
167
0
    int j = 0;  // Used after for.
168
0
    for (int i = 1; pszDatum[i] != '\0'; i++)
169
0
    {
170
0
        if (pszDatum[j] == '_' && pszDatum[i] == '_')
171
0
            continue;
172
173
0
        pszDatum[++j] = pszDatum[i];
174
0
    }
175
0
    if (pszDatum[j] == '_')
176
0
        pszDatum[j] = '\0';
177
0
    else
178
0
        pszDatum[j + 1] = '\0';
179
180
    /* -------------------------------------------------------------------- */
181
    /*      Search for datum equivalences.  Specific massaged names get     */
182
    /*      mapped to OpenGIS specified names.                              */
183
    /* -------------------------------------------------------------------- */
184
0
    for (int i = 0; papszDatumEquiv[i] != nullptr; i += 2)
185
0
    {
186
0
        if (EQUAL(*ppszDatum, papszDatumEquiv[i]))
187
0
        {
188
0
            CPLFree(*ppszDatum);
189
0
            *ppszDatum = CPLStrdup(papszDatumEquiv[i + 1]);
190
0
            return;
191
0
        }
192
0
    }
193
0
}
194
195
/************************************************************************/
196
/*                      GTIFCleanupImageineNames()                      */
197
/*                                                                      */
198
/*      Erdas Imagine sometimes emits big copyright messages, and       */
199
/*      other stuff into citations.  These can be pretty messy when     */
200
/*      turned into WKT, so we try to trim and clean the strings        */
201
/*      somewhat.                                                       */
202
/************************************************************************/
203
204
/* For example:
205
   GTCitationGeoKey (Ascii,215): "IMAGINE GeoTIFF Support\nCopyright 1991 - 2001
206
   by ERDAS, Inc. All Rights Reserved\n@(#)$RCSfile$ $Revision: 34309 $ $Date:
207
   2016-05-29 11:29:40 -0700 (Sun, 29 May 2016) $\nProjection Name = UTM\nUnits
208
   = meters\nGeoTIFF Units = meters"
209
210
   GeogCitationGeoKey (Ascii,267): "IMAGINE GeoTIFF Support\nCopyright 1991 -
211
   2001 by ERDAS, Inc. All Rights Reserved\n@(#)$RCSfile$ $Revision: 34309 $
212
   $Date: 2016-05-29 11:29:40 -0700 (Sun, 29 May 2016) $\nUnable to match
213
   Ellipsoid (Datum) to a GeographicTypeGeoKey value\nEllipsoid = Clarke
214
   1866\nDatum = NAD27 (CONUS)"
215
216
   PCSCitationGeoKey (Ascii,214): "IMAGINE GeoTIFF Support\nCopyright 1991 -
217
   2001 by ERDAS, Inc. All Rights Reserved\n@(#)$RCSfile$ $Revision: 34309 $
218
   $Date: 2016-05-29 11:29:40 -0700 (Sun, 29 May 2016) $\nUTM Zone
219
   10N\nEllipsoid = Clarke 1866\nDatum = NAD27 (CONUS)"
220
*/
221
222
static void GTIFCleanupImagineNames(char *pszCitation)
223
224
0
{
225
0
    if (strstr(pszCitation, "IMAGINE GeoTIFF") == nullptr)
226
0
        return;
227
228
    /* -------------------------------------------------------------------- */
229
    /*      First, we skip past all the copyright, and RCS stuff.  We       */
230
    /*      assume that this will have a "$" at the end of it all.          */
231
    /* -------------------------------------------------------------------- */
232
0
    char *pszSkip = pszCitation + strlen(pszCitation) - 1;
233
234
0
    for (; pszSkip != pszCitation && *pszSkip != '$'; pszSkip--)
235
0
    {
236
0
    }
237
238
0
    if (*pszSkip == '$')
239
0
        pszSkip++;
240
0
    if (*pszSkip == '\n')
241
0
        pszSkip++;
242
243
0
    memmove(pszCitation, pszSkip, strlen(pszSkip) + 1);
244
245
    /* -------------------------------------------------------------------- */
246
    /*      Convert any newlines into spaces, they really gum up the        */
247
    /*      WKT.                                                            */
248
    /* -------------------------------------------------------------------- */
249
0
    for (int i = 0; pszCitation[i] != '\0'; i++)
250
0
    {
251
0
        if (pszCitation[i] == '\n')
252
0
            pszCitation[i] = ' ';
253
0
    }
254
0
}
255
256
#if LIBGEOTIFF_VERSION < 1600
257
258
/************************************************************************/
259
/*                       GDALGTIFKeyGet()                               */
260
/************************************************************************/
261
262
static int GDALGTIFKeyGet(GTIF *hGTIF, geokey_t key, void *pData, int nIndex,
263
                          int nCount, tagtype_t expected_tagtype)
264
{
265
    tagtype_t tagtype = TYPE_UNKNOWN;
266
    if (!GTIFKeyInfo(hGTIF, key, nullptr, &tagtype))
267
        return 0;
268
    if (tagtype != expected_tagtype)
269
    {
270
        CPLError(CE_Warning, CPLE_AppDefined,
271
                 "Expected key %s to be of type %s. Got %s", GTIFKeyName(key),
272
                 GTIFTypeName(expected_tagtype), GTIFTypeName(tagtype));
273
        return 0;
274
    }
275
    return GTIFKeyGet(hGTIF, key, pData, nIndex, nCount);
276
}
277
278
/************************************************************************/
279
/*                       GDALGTIFKeyGetASCII()                          */
280
/************************************************************************/
281
282
int GDALGTIFKeyGetASCII(GTIF *hGTIF, geokey_t key, char *szStr, int szStrMaxLen)
283
{
284
    return GDALGTIFKeyGet(hGTIF, key, szStr, 0, szStrMaxLen, TYPE_ASCII);
285
}
286
287
/************************************************************************/
288
/*                       GDALGTIFKeyGetSHORT()                          */
289
/************************************************************************/
290
291
int GDALGTIFKeyGetSHORT(GTIF *hGTIF, geokey_t key, unsigned short *pnVal,
292
                        int nIndex, int nCount)
293
{
294
    return GDALGTIFKeyGet(hGTIF, key, pnVal, nIndex, nCount, TYPE_SHORT);
295
}
296
297
/************************************************************************/
298
/*                        GDALGTIFKeyGetDOUBLE()                        */
299
/************************************************************************/
300
301
int GDALGTIFKeyGetDOUBLE(GTIF *hGTIF, geokey_t key, double *pdfVal, int nIndex,
302
                         int nCount)
303
{
304
    return GDALGTIFKeyGet(hGTIF, key, pdfVal, nIndex, nCount, TYPE_DOUBLE);
305
}
306
307
#endif
308
309
/************************************************************************/
310
/*                    FillCompoundCRSWithManualVertCS()                 */
311
/************************************************************************/
312
313
static void FillCompoundCRSWithManualVertCS(GTIF *hGTIF,
314
                                            OGRSpatialReference &oSRS,
315
                                            const char *pszVertCSName,
316
                                            int verticalDatum,
317
                                            int verticalUnits)
318
0
{
319
    /* -------------------------------------------------------------------- */
320
    /*      Setup VERT_CS with citation if present.                         */
321
    /* -------------------------------------------------------------------- */
322
0
    oSRS.SetNode("COMPD_CS|VERT_CS", pszVertCSName);
323
324
    /* -------------------------------------------------------------------- */
325
    /*      Setup the vertical datum.                                       */
326
    /* -------------------------------------------------------------------- */
327
0
    std::string osVDatumName = "unknown";
328
0
    const char *pszVDatumType = "2005";  // CS_VD_GeoidModelDerived
329
0
    std::string osVDatumAuthName;
330
0
    int nVDatumCode = 0;
331
332
0
    if (verticalDatum > 0 && verticalDatum != KvUserDefined)
333
0
    {
334
0
        osVDatumAuthName = "EPSG";
335
0
        nVDatumCode = verticalDatum;
336
337
0
        char szCode[12];
338
0
        snprintf(szCode, sizeof(szCode), "%d", verticalDatum);
339
0
        auto ctx =
340
0
            static_cast<PJ_CONTEXT *>(GTIFGetPROJContext(hGTIF, true, nullptr));
341
0
        auto datum = proj_create_from_database(ctx, "EPSG", szCode,
342
0
                                               PJ_CATEGORY_DATUM, 0, nullptr);
343
0
        if (datum)
344
0
        {
345
0
            const char *pszName = proj_get_name(datum);
346
0
            if (pszName)
347
0
            {
348
0
                osVDatumName = pszName;
349
0
            }
350
0
            proj_destroy(datum);
351
0
        }
352
0
    }
353
0
    else if (verticalDatum == KvUserDefined)
354
0
    {
355
        // If the vertical datum is unknown, try to find the vertical CRS
356
        // from the database, and extra the datum information from it.
357
0
        auto ctx =
358
0
            static_cast<PJ_CONTEXT *>(GTIFGetPROJContext(hGTIF, true, nullptr));
359
0
        const auto type = PJ_TYPE_VERTICAL_CRS;
360
0
        auto list = proj_create_from_name(ctx, nullptr, pszVertCSName, &type, 1,
361
0
                                          true,  // exact match
362
0
                                          -1,    // result set limit size,
363
0
                                          nullptr);
364
0
        if (list)
365
0
        {
366
            // If we have several matches, check they all refer to the
367
            // same datum
368
0
            bool bGoOn = true;
369
0
            int ncount = proj_list_get_count(list);
370
0
            for (int i = 0; bGoOn && i < ncount; ++i)
371
0
            {
372
0
                auto crs = proj_list_get(ctx, list, i);
373
0
                if (crs)
374
0
                {
375
0
                    auto datum = proj_crs_get_datum(ctx, crs);
376
0
                    if (datum)
377
0
                    {
378
0
                        osVDatumName = proj_get_name(datum);
379
0
                        const char *pszAuthName =
380
0
                            proj_get_id_auth_name(datum, 0);
381
0
                        const char *pszCode = proj_get_id_code(datum, 0);
382
0
                        if (pszCode && atoi(pszCode) && pszAuthName)
383
0
                        {
384
0
                            if (osVDatumAuthName.empty())
385
0
                            {
386
0
                                osVDatumAuthName = pszAuthName;
387
0
                                nVDatumCode = atoi(pszCode);
388
0
                            }
389
0
                            else if (osVDatumAuthName != pszAuthName ||
390
0
                                     nVDatumCode != atoi(pszCode))
391
0
                            {
392
0
                                osVDatumAuthName.clear();
393
0
                                nVDatumCode = 0;
394
0
                                bGoOn = false;
395
0
                            }
396
0
                        }
397
0
                        proj_destroy(datum);
398
0
                    }
399
0
                    proj_destroy(crs);
400
0
                }
401
0
            }
402
0
        }
403
0
        proj_list_destroy(list);
404
0
    }
405
406
0
    oSRS.SetNode("COMPD_CS|VERT_CS|VERT_DATUM", osVDatumName.c_str());
407
0
    oSRS.GetAttrNode("COMPD_CS|VERT_CS|VERT_DATUM")
408
0
        ->AddChild(new OGR_SRSNode(pszVDatumType));
409
0
    if (!osVDatumAuthName.empty())
410
0
        oSRS.SetAuthority("COMPD_CS|VERT_CS|VERT_DATUM",
411
0
                          osVDatumAuthName.c_str(), nVDatumCode);
412
413
    /* -------------------------------------------------------------------- */
414
    /*      Set the vertical units.                                         */
415
    /* -------------------------------------------------------------------- */
416
0
    if (verticalUnits > 0 && verticalUnits != KvUserDefined &&
417
0
        verticalUnits != 9001)
418
0
    {
419
0
        char szCode[12];
420
0
        snprintf(szCode, sizeof(szCode), "%d", verticalUnits);
421
0
        auto ctx =
422
0
            static_cast<PJ_CONTEXT *>(GTIFGetPROJContext(hGTIF, true, nullptr));
423
0
        const char *pszName = nullptr;
424
0
        double dfInMeters = 0.0;
425
0
        if (proj_uom_get_info_from_database(ctx, "EPSG", szCode, &pszName,
426
0
                                            &dfInMeters, nullptr))
427
0
        {
428
0
            if (pszName)
429
0
                oSRS.SetNode("COMPD_CS|VERT_CS|UNIT", pszName);
430
431
0
            char szInMeters[128] = {};
432
0
            CPLsnprintf(szInMeters, sizeof(szInMeters), "%.16g", dfInMeters);
433
0
            oSRS.GetAttrNode("COMPD_CS|VERT_CS|UNIT")
434
0
                ->AddChild(new OGR_SRSNode(szInMeters));
435
0
        }
436
437
0
        oSRS.SetAuthority("COMPD_CS|VERT_CS|UNIT", "EPSG", verticalUnits);
438
0
    }
439
0
    else
440
0
    {
441
0
        oSRS.SetNode("COMPD_CS|VERT_CS|UNIT", "metre");
442
0
        oSRS.GetAttrNode("COMPD_CS|VERT_CS|UNIT")
443
0
            ->AddChild(new OGR_SRSNode("1.0"));
444
0
        oSRS.SetAuthority("COMPD_CS|VERT_CS|UNIT", "EPSG", 9001);
445
0
    }
446
447
    /* -------------------------------------------------------------------- */
448
    /*      Set the axis and VERT_CS authority.                             */
449
    /* -------------------------------------------------------------------- */
450
0
    oSRS.SetNode("COMPD_CS|VERT_CS|AXIS", "Up");
451
0
    oSRS.GetAttrNode("COMPD_CS|VERT_CS|AXIS")->AddChild(new OGR_SRSNode("UP"));
452
0
}
453
454
/************************************************************************/
455
/*                    GTIFGetEPSGOfficialName()                         */
456
/************************************************************************/
457
458
static char *GTIFGetEPSGOfficialName(GTIF *hGTIF, PJ_TYPE searchType,
459
                                     const char *pszName)
460
0
{
461
0
    char *pszRet = nullptr;
462
    /* Search in database the corresponding EPSG 'official' name */
463
0
    auto ctx =
464
0
        static_cast<PJ_CONTEXT *>(GTIFGetPROJContext(hGTIF, true, nullptr));
465
0
    auto list = proj_create_from_name(ctx, "EPSG", pszName, &searchType, 1,
466
0
                                      false, /* approximate match */
467
0
                                      1, nullptr);
468
0
    if (list)
469
0
    {
470
0
        const auto listSize = proj_list_get_count(list);
471
0
        if (listSize == 1)
472
0
        {
473
0
            auto obj = proj_list_get(ctx, list, 0);
474
0
            if (obj)
475
0
            {
476
0
                const char *pszOfficialName = proj_get_name(obj);
477
0
                if (pszOfficialName)
478
0
                {
479
0
                    pszRet = CPLStrdup(pszOfficialName);
480
0
                }
481
0
            }
482
0
            proj_destroy(obj);
483
0
        }
484
0
        proj_list_destroy(list);
485
0
    }
486
0
    return pszRet;
487
0
}
488
489
/************************************************************************/
490
/*                      GTIFGetOGISDefnAsOSR()                          */
491
/************************************************************************/
492
493
OGRSpatialReferenceH GTIFGetOGISDefnAsOSR(GTIF *hGTIF, GTIFDefn *psDefn)
494
495
0
{
496
0
    OGRSpatialReference oSRS;
497
498
0
    LibgeotiffOneTimeInit();
499
500
0
#if LIBGEOTIFF_VERSION >= 1600
501
0
    void *projContext = GTIFGetPROJContext(hGTIF, FALSE, nullptr);
502
0
#endif
503
504
    /* -------------------------------------------------------------------- */
505
    /*  Handle non-standard coordinate systems where GTModelTypeGeoKey      */
506
    /*  is not defined, but ProjectedCSTypeGeoKey is defined (ticket #3019) */
507
    /* -------------------------------------------------------------------- */
508
0
    if (psDefn->Model == KvUserDefined && psDefn->PCS != KvUserDefined)
509
0
    {
510
0
        psDefn->Model = ModelTypeProjected;
511
0
    }
512
513
    /* ==================================================================== */
514
    /*      Read keys related to vertical component.                        */
515
    /* ==================================================================== */
516
0
    unsigned short verticalCSType = 0;
517
0
    unsigned short verticalDatum = 0;
518
0
    unsigned short verticalUnits = 0;
519
520
0
    GDALGTIFKeyGetSHORT(hGTIF, VerticalCSTypeGeoKey, &verticalCSType, 0, 1);
521
0
    GDALGTIFKeyGetSHORT(hGTIF, VerticalDatumGeoKey, &verticalDatum, 0, 1);
522
0
    GDALGTIFKeyGetSHORT(hGTIF, VerticalUnitsGeoKey, &verticalUnits, 0, 1);
523
524
0
    if (verticalCSType != 0 || verticalDatum != 0 || verticalUnits != 0)
525
0
    {
526
0
        int versions[3];
527
0
        GTIFDirectoryInfo(hGTIF, versions, nullptr);
528
        // GeoTIFF 1.0
529
0
        if (versions[0] == 1 && versions[1] == 1 && versions[2] == 0)
530
0
        {
531
            /* --------------------------------------------------------------------
532
             */
533
            /*      The original geotiff specification appears to have */
534
            /*      misconstrued the EPSG codes 5101 to 5106 to be vertical */
535
            /*      coordinate system codes, when in fact they are vertical */
536
            /*      datum codes.  So if these are found in the */
537
            /*      VerticalCSTypeGeoKey move them to the VerticalDatumGeoKey */
538
            /*      and insert the "normal" corresponding VerticalCSTypeGeoKey
539
             */
540
            /*      value. */
541
            /* --------------------------------------------------------------------
542
             */
543
0
            if ((verticalCSType >= 5101 && verticalCSType <= 5112) &&
544
0
                verticalDatum == 0)
545
0
            {
546
0
                verticalDatum = verticalCSType;
547
0
                verticalCSType = verticalDatum + 600;
548
0
            }
549
550
            /* --------------------------------------------------------------------
551
             */
552
            /*      This addresses another case where the EGM96 Vertical Datum
553
             * code */
554
            /*      is misused as a Vertical CS code (#4922). */
555
            /* --------------------------------------------------------------------
556
             */
557
0
            if (verticalCSType == 5171)
558
0
            {
559
0
                verticalDatum = 5171;
560
0
                verticalCSType = 5773;
561
0
            }
562
0
        }
563
564
        /* --------------------------------------------------------------------
565
         */
566
        /*      Somewhat similarly, codes 5001 to 5033 were treated as */
567
        /*      vertical coordinate systems based on ellipsoidal heights. */
568
        /*      We use the corresponding geodetic datum as the vertical */
569
        /*      datum and clear the vertical coordinate system code since */
570
        /*      there isn't one in EPSG. */
571
        /* --------------------------------------------------------------------
572
         */
573
0
        if ((verticalCSType >= 5001 && verticalCSType <= 5033) &&
574
0
            verticalDatum == 0)
575
0
        {
576
0
            verticalDatum = verticalCSType + 1000;
577
0
            verticalCSType = 0;
578
0
        }
579
0
    }
580
581
    /* -------------------------------------------------------------------- */
582
    /*      Handle non-standard coordinate systems as LOCAL_CS.             */
583
    /* -------------------------------------------------------------------- */
584
0
    if (psDefn->Model != ModelTypeProjected &&
585
0
        psDefn->Model != ModelTypeGeographic &&
586
0
        psDefn->Model != ModelTypeGeocentric)
587
0
    {
588
0
        char szPeStr[2400] = {'\0'};
589
590
        /** check if there is a pe string citation key **/
591
0
        if (GDALGTIFKeyGetASCII(hGTIF, PCSCitationGeoKey, szPeStr,
592
0
                                sizeof(szPeStr)) &&
593
0
            strstr(szPeStr, "ESRI PE String = "))
594
0
        {
595
0
            const char *pszWKT = szPeStr + strlen("ESRI PE String = ");
596
0
            oSRS.importFromWkt(pszWKT);
597
598
0
            if (strstr(pszWKT,
599
0
                       "PROJCS[\"WGS_1984_Web_Mercator_Auxiliary_Sphere\""))
600
0
            {
601
0
                oSRS.SetExtension(
602
0
                    "PROJCS", "PROJ4",
603
0
                    "+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 "
604
0
                    "+x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null "
605
0
                    "+wktext  +no_defs");  // TODO(schwehr): Why 2 spaces?
606
0
            }
607
608
0
            return OGRSpatialReference::ToHandle(oSRS.Clone());
609
0
        }
610
0
        else
611
0
        {
612
0
            char *pszUnitsName = nullptr;
613
0
            char szPCSName[300] = {'\0'};
614
0
            int nKeyCount = 0;
615
0
            int anVersion[3] = {0};
616
617
0
            GTIFDirectoryInfo(hGTIF, anVersion, &nKeyCount);
618
619
0
            if (nKeyCount > 0)  // Use LOCAL_CS if we have any geokeys at all.
620
0
            {
621
                // Handle citation.
622
0
                strcpy(szPCSName, "unnamed");
623
0
                if (!GDALGTIFKeyGetASCII(hGTIF, GTCitationGeoKey, szPCSName,
624
0
                                         sizeof(szPCSName)))
625
0
                    GDALGTIFKeyGetASCII(hGTIF, GeogCitationGeoKey, szPCSName,
626
0
                                        sizeof(szPCSName));
627
628
0
                GTIFCleanupImagineNames(szPCSName);
629
0
                oSRS.SetLocalCS(szPCSName);
630
631
                // Handle units
632
0
                if (psDefn->UOMLength != KvUserDefined)
633
0
                {
634
0
#if LIBGEOTIFF_VERSION >= 1600
635
0
                    GTIFGetUOMLengthInfoEx(projContext,
636
#else
637
                    GTIFGetUOMLengthInfo(
638
#endif
639
0
                                           psDefn->UOMLength, &pszUnitsName,
640
0
                                           nullptr);
641
0
                }
642
643
0
                if (pszUnitsName != nullptr)
644
0
                {
645
0
                    char szUOMLength[12];
646
0
                    snprintf(szUOMLength, sizeof(szUOMLength), "%d",
647
0
                             psDefn->UOMLength);
648
0
                    oSRS.SetTargetLinearUnits(nullptr, pszUnitsName,
649
0
                                              psDefn->UOMLengthInMeters, "EPSG",
650
0
                                              szUOMLength);
651
0
                }
652
0
                else
653
0
                    oSRS.SetLinearUnits("unknown", psDefn->UOMLengthInMeters);
654
655
0
                if (verticalUnits != 0)
656
0
                {
657
0
                    char szVertCSCitation[2048] = {0};
658
0
                    if (GDALGTIFKeyGetASCII(hGTIF, VerticalCitationGeoKey,
659
0
                                            szVertCSCitation,
660
0
                                            sizeof(szVertCSCitation)))
661
0
                    {
662
0
                        if (STARTS_WITH_CI(szVertCSCitation, "VCS Name = "))
663
0
                        {
664
0
                            memmove(szVertCSCitation,
665
0
                                    szVertCSCitation + strlen("VCS Name = "),
666
0
                                    strlen(szVertCSCitation +
667
0
                                           strlen("VCS Name = ")) +
668
0
                                        1);
669
0
                            char *pszPipeChar = strchr(szVertCSCitation, '|');
670
0
                            if (pszPipeChar)
671
0
                                *pszPipeChar = '\0';
672
0
                        }
673
0
                    }
674
0
                    else
675
0
                    {
676
0
                        strcpy(szVertCSCitation, "unknown");
677
0
                    }
678
679
0
                    const char *pszHorizontalName = oSRS.GetName();
680
0
                    const std::string osHorizontalName(
681
0
                        pszHorizontalName ? pszHorizontalName : "unnamed");
682
                    /* --------------------------------------------------------------------
683
                     */
684
                    /*      Promote to being a compound coordinate system. */
685
                    /* --------------------------------------------------------------------
686
                     */
687
0
                    OGR_SRSNode *poOldRoot = oSRS.GetRoot()->Clone();
688
689
0
                    oSRS.Clear();
690
691
                    /* --------------------------------------------------------------------
692
                     */
693
                    /*      Set COMPD_CS name. */
694
                    /* --------------------------------------------------------------------
695
                     */
696
0
                    char szCTString[512];
697
0
                    szCTString[0] = '\0';
698
0
                    if (GDALGTIFKeyGetASCII(hGTIF, GTCitationGeoKey, szCTString,
699
0
                                            sizeof(szCTString)) &&
700
0
                        strstr(szCTString, " = ") == nullptr)
701
0
                    {
702
0
                        oSRS.SetNode("COMPD_CS", szCTString);
703
0
                    }
704
0
                    else
705
0
                    {
706
0
                        oSRS.SetNode("COMPD_CS", (osHorizontalName + " + " +
707
0
                                                  szVertCSCitation)
708
0
                                                     .c_str());
709
0
                    }
710
711
0
                    oSRS.GetRoot()->AddChild(poOldRoot);
712
713
0
                    FillCompoundCRSWithManualVertCS(
714
0
                        hGTIF, oSRS, szVertCSCitation, verticalDatum,
715
0
                        verticalUnits);
716
0
                }
717
718
0
                GTIFFreeMemory(pszUnitsName);
719
0
            }
720
0
            return OGRSpatialReference::ToHandle(oSRS.Clone());
721
0
        }
722
0
    }
723
724
    /* -------------------------------------------------------------------- */
725
    /*      Handle Geocentric coordinate systems.                           */
726
    /* -------------------------------------------------------------------- */
727
0
    if (psDefn->Model == ModelTypeGeocentric)
728
0
    {
729
0
        char szName[300] = {'\0'};
730
731
0
        strcpy(szName, "unnamed");
732
0
        if (!GDALGTIFKeyGetASCII(hGTIF, GTCitationGeoKey, szName,
733
0
                                 sizeof(szName)))
734
0
            GDALGTIFKeyGetASCII(hGTIF, GeogCitationGeoKey, szName,
735
0
                                sizeof(szName));
736
737
0
        oSRS.SetGeocCS(szName);
738
739
0
        char *pszUnitsName = nullptr;
740
741
0
        if (psDefn->UOMLength != KvUserDefined)
742
0
        {
743
0
#if LIBGEOTIFF_VERSION >= 1600
744
0
            GTIFGetUOMLengthInfoEx(projContext,
745
#else
746
            GTIFGetUOMLengthInfo(
747
#endif
748
0
                                   psDefn->UOMLength, &pszUnitsName, nullptr);
749
0
        }
750
751
0
        if (pszUnitsName != nullptr)
752
0
        {
753
0
            char szUOMLength[12];
754
0
            snprintf(szUOMLength, sizeof(szUOMLength), "%d", psDefn->UOMLength);
755
0
            oSRS.SetTargetLinearUnits(nullptr, pszUnitsName,
756
0
                                      psDefn->UOMLengthInMeters, "EPSG",
757
0
                                      szUOMLength);
758
0
        }
759
0
        else
760
0
            oSRS.SetLinearUnits("unknown", psDefn->UOMLengthInMeters);
761
762
0
        GTIFFreeMemory(pszUnitsName);
763
0
    }
764
765
    /* -------------------------------------------------------------------- */
766
    /*      #3901: In libgeotiff 1.3.0 and earlier we incorrectly           */
767
    /*      interpreted linear projection parameter geokeys (false          */
768
    /*      easting/northing) as being in meters instead of the             */
769
    /*      coordinate system of the file.   The following code attempts    */
770
    /*      to provide mechanisms for fixing the issue if we are linked     */
771
    /*      with an older version of libgeotiff.                            */
772
    /* -------------------------------------------------------------------- */
773
0
    const char *pszLinearUnits =
774
0
        CPLGetConfigOption("GTIFF_LINEAR_UNITS", "DEFAULT");
775
776
    /* -------------------------------------------------------------------- */
777
    /*      #3901: If folks have broken GeoTIFF files generated with        */
778
    /*      older versions of GDAL+libgeotiff, then they may need a         */
779
    /*      hack to allow them to be read properly.  This is that           */
780
    /*      hack.  We basically try to undue the conversion applied by      */
781
    /*      libgeotiff to meters (or above) to simulate the old             */
782
    /*      behavior.                                                       */
783
    /* -------------------------------------------------------------------- */
784
0
    unsigned short bLinearUnitsMarkedCorrect = FALSE;
785
786
0
    GDALGTIFKeyGetSHORT(hGTIF, ProjLinearUnitsInterpCorrectGeoKey,
787
0
                        &bLinearUnitsMarkedCorrect, 0, 1);
788
789
0
    if (EQUAL(pszLinearUnits, "BROKEN") &&
790
0
        psDefn->Projection == KvUserDefined && !bLinearUnitsMarkedCorrect)
791
0
    {
792
0
        for (int iParam = 0; iParam < psDefn->nParms; iParam++)
793
0
        {
794
0
            switch (psDefn->ProjParmId[iParam])
795
0
            {
796
0
                case ProjFalseEastingGeoKey:
797
0
                case ProjFalseNorthingGeoKey:
798
0
                case ProjFalseOriginEastingGeoKey:
799
0
                case ProjFalseOriginNorthingGeoKey:
800
0
                case ProjCenterEastingGeoKey:
801
0
                case ProjCenterNorthingGeoKey:
802
0
                    if (psDefn->UOMLengthInMeters != 0 &&
803
0
                        psDefn->UOMLengthInMeters != 1.0)
804
0
                    {
805
0
                        psDefn->ProjParm[iParam] /= psDefn->UOMLengthInMeters;
806
0
                        CPLDebug(
807
0
                            "GTIFF",
808
0
                            "Converting geokey to accommodate old broken file "
809
0
                            "due to GTIFF_LINEAR_UNITS=BROKEN setting.");
810
0
                    }
811
0
                    break;
812
813
0
                default:
814
0
                    break;
815
0
            }
816
0
        }
817
0
    }
818
819
    /* -------------------------------------------------------------------- */
820
    /*      If this is a projected SRS we set the PROJCS keyword first      */
821
    /*      to ensure that the GEOGCS will be a child.                      */
822
    /* -------------------------------------------------------------------- */
823
0
    OGRBoolean linearUnitIsSet = FALSE;
824
0
    if (psDefn->Model == ModelTypeProjected)
825
0
    {
826
0
        char szCTString[512] = {'\0'};
827
0
        if (psDefn->PCS != KvUserDefined)
828
0
        {
829
0
            char *pszPCSName = nullptr;
830
831
0
#if LIBGEOTIFF_VERSION >= 1600
832
0
            GTIFGetPCSInfoEx(projContext,
833
#else
834
            GTIFGetPCSInfo(
835
#endif
836
0
                             psDefn->PCS, &pszPCSName, nullptr, nullptr,
837
0
                             nullptr);
838
839
0
            oSRS.SetProjCS(pszPCSName ? pszPCSName : "unnamed");
840
0
            if (pszPCSName)
841
0
                GTIFFreeMemory(pszPCSName);
842
843
0
            oSRS.SetLinearUnits("unknown", 1.0);
844
0
        }
845
0
        else
846
0
        {
847
0
            bool bTryGTCitationGeoKey = true;
848
0
            if (GDALGTIFKeyGetASCII(hGTIF, PCSCitationGeoKey, szCTString,
849
0
                                    sizeof(szCTString)))
850
0
            {
851
0
                bTryGTCitationGeoKey = false;
852
0
                if (!SetCitationToSRS(hGTIF, szCTString, sizeof(szCTString),
853
0
                                      PCSCitationGeoKey, &oSRS,
854
0
                                      &linearUnitIsSet))
855
0
                {
856
0
                    if (!STARTS_WITH_CI(szCTString, "LUnits = "))
857
0
                    {
858
0
                        oSRS.SetProjCS(szCTString);
859
0
                        oSRS.SetLinearUnits("unknown", 1.0);
860
0
                    }
861
0
                    else
862
0
                    {
863
0
                        bTryGTCitationGeoKey = true;
864
0
                    }
865
0
                }
866
0
            }
867
868
0
            if (bTryGTCitationGeoKey)
869
0
            {
870
0
                if (GDALGTIFKeyGetASCII(hGTIF, GTCitationGeoKey, szCTString,
871
0
                                        sizeof(szCTString)) &&
872
0
                    !SetCitationToSRS(hGTIF, szCTString, sizeof(szCTString),
873
0
                                      GTCitationGeoKey, &oSRS,
874
0
                                      &linearUnitIsSet))
875
0
                {
876
0
                    oSRS.SetNode("PROJCS", szCTString);
877
0
                    oSRS.SetLinearUnits("unknown", 1.0);
878
0
                }
879
0
                else
880
0
                {
881
0
                    oSRS.SetNode("PROJCS", "unnamed");
882
0
                    oSRS.SetLinearUnits("unknown", 1.0);
883
0
                }
884
0
            }
885
0
        }
886
887
        /* Handle ESRI/Erdas style state plane and UTM in citation key */
888
0
        if (CheckCitationKeyForStatePlaneUTM(hGTIF, psDefn, &oSRS,
889
0
                                             &linearUnitIsSet))
890
0
        {
891
0
            return OGRSpatialReference::ToHandle(oSRS.Clone());
892
0
        }
893
894
        /* Handle ESRI PE string in citation */
895
0
        szCTString[0] = '\0';
896
0
        if (GDALGTIFKeyGetASCII(hGTIF, GTCitationGeoKey, szCTString,
897
0
                                sizeof(szCTString)))
898
0
            SetCitationToSRS(hGTIF, szCTString, sizeof(szCTString),
899
0
                             GTCitationGeoKey, &oSRS, &linearUnitIsSet);
900
0
    }
901
902
    /* ==================================================================== */
903
    /*      Setup the GeogCS                                                */
904
    /* ==================================================================== */
905
0
    char *pszGeogName = nullptr;
906
0
    char *pszDatumName = nullptr;
907
0
    char *pszPMName = nullptr;
908
0
    char *pszSpheroidName = nullptr;
909
0
    char *pszAngularUnits = nullptr;
910
0
    char szGCSName[512] = {'\0'};
911
912
0
    if (!
913
0
#if LIBGEOTIFF_VERSION >= 1600
914
0
        GTIFGetGCSInfoEx(projContext,
915
#else
916
        GTIFGetGCSInfo(
917
#endif
918
0
                         psDefn->GCS, &pszGeogName, nullptr, nullptr,
919
0
                         nullptr) &&
920
0
        GDALGTIFKeyGetASCII(hGTIF, GeogCitationGeoKey, szGCSName,
921
0
                            sizeof(szGCSName)))
922
0
    {
923
0
        GetGeogCSFromCitation(szGCSName, sizeof(szGCSName), GeogCitationGeoKey,
924
0
                              &pszGeogName, &pszDatumName, &pszPMName,
925
0
                              &pszSpheroidName, &pszAngularUnits);
926
0
    }
927
0
    else
928
0
    {
929
0
        GTIFToCPLRecycleString(&pszGeogName);
930
0
    }
931
932
0
    if (pszGeogName && STARTS_WITH(pszGeogName, "GCS_"))
933
0
    {
934
        // Morph from potential ESRI name
935
0
        char *pszOfficialName = GTIFGetEPSGOfficialName(
936
0
            hGTIF, PJ_TYPE_GEOGRAPHIC_2D_CRS, pszGeogName);
937
0
        if (pszOfficialName)
938
0
        {
939
0
            CPLFree(pszGeogName);
940
0
            pszGeogName = pszOfficialName;
941
0
        }
942
0
    }
943
944
0
    if (pszDatumName && strchr(pszDatumName, '_'))
945
0
    {
946
        // Morph from potential ESRI name
947
0
        char *pszOfficialName = GTIFGetEPSGOfficialName(
948
0
            hGTIF, PJ_TYPE_GEODETIC_REFERENCE_FRAME, pszDatumName);
949
0
        if (pszOfficialName)
950
0
        {
951
0
            CPLFree(pszDatumName);
952
0
            pszDatumName = pszOfficialName;
953
0
        }
954
0
    }
955
956
0
    if (pszSpheroidName && strchr(pszSpheroidName, '_'))
957
0
    {
958
        // Morph from potential ESRI name
959
0
        char *pszOfficialName =
960
0
            GTIFGetEPSGOfficialName(hGTIF, PJ_TYPE_ELLIPSOID, pszSpheroidName);
961
0
        if (pszOfficialName)
962
0
        {
963
0
            CPLFree(pszSpheroidName);
964
0
            pszSpheroidName = pszOfficialName;
965
0
        }
966
0
    }
967
968
0
    if (!pszDatumName)
969
0
    {
970
0
#if LIBGEOTIFF_VERSION >= 1600
971
0
        GTIFGetDatumInfoEx(projContext,
972
#else
973
        GTIFGetDatumInfo(
974
#endif
975
0
                           psDefn->Datum, &pszDatumName, nullptr);
976
0
        GTIFToCPLRecycleString(&pszDatumName);
977
0
    }
978
979
0
    double dfSemiMajor = 0.0;
980
0
    double dfInvFlattening = 0.0;
981
0
    if (!pszSpheroidName)
982
0
    {
983
0
#if LIBGEOTIFF_VERSION >= 1600
984
0
        GTIFGetEllipsoidInfoEx(projContext,
985
#else
986
        GTIFGetEllipsoidInfo(
987
#endif
988
0
                               psDefn->Ellipsoid, &pszSpheroidName, nullptr,
989
0
                               nullptr);
990
0
        GTIFToCPLRecycleString(&pszSpheroidName);
991
0
    }
992
0
    else
993
0
    {
994
0
        CPL_IGNORE_RET_VAL(GDALGTIFKeyGetDOUBLE(hGTIF, GeogSemiMajorAxisGeoKey,
995
0
                                                &(psDefn->SemiMajor), 0, 1));
996
0
        CPL_IGNORE_RET_VAL(GDALGTIFKeyGetDOUBLE(hGTIF, GeogInvFlatteningGeoKey,
997
0
                                                &dfInvFlattening, 0, 1));
998
0
        if (std::isinf(dfInvFlattening))
999
0
        {
1000
            // Deal with the non-nominal case of
1001
            // https://github.com/OSGeo/PROJ/issues/2317
1002
0
            dfInvFlattening = 0;
1003
0
        }
1004
0
    }
1005
0
    if (!pszPMName)
1006
0
    {
1007
0
#if LIBGEOTIFF_VERSION >= 1600
1008
0
        GTIFGetPMInfoEx(projContext,
1009
#else
1010
        GTIFGetPMInfo(
1011
#endif
1012
0
                        psDefn->PM, &pszPMName, nullptr);
1013
0
        GTIFToCPLRecycleString(&pszPMName);
1014
0
    }
1015
0
    else
1016
0
    {
1017
0
        CPL_IGNORE_RET_VAL(
1018
0
            GDALGTIFKeyGetDOUBLE(hGTIF, GeogPrimeMeridianLongGeoKey,
1019
0
                                 &(psDefn->PMLongToGreenwich), 0, 1));
1020
0
    }
1021
1022
0
    if (!pszAngularUnits)
1023
0
    {
1024
0
#if LIBGEOTIFF_VERSION >= 1600
1025
0
        GTIFGetUOMAngleInfoEx(projContext,
1026
#else
1027
        GTIFGetUOMAngleInfo(
1028
#endif
1029
0
                              psDefn->UOMAngle, &pszAngularUnits,
1030
0
                              &psDefn->UOMAngleInDegrees);
1031
0
        if (pszAngularUnits == nullptr)
1032
0
            pszAngularUnits = CPLStrdup("unknown");
1033
0
        else
1034
0
            GTIFToCPLRecycleString(&pszAngularUnits);
1035
0
    }
1036
0
    else
1037
0
    {
1038
0
        double dfRadians = 0.0;
1039
0
        if (GDALGTIFKeyGetDOUBLE(hGTIF, GeogAngularUnitSizeGeoKey, &dfRadians,
1040
0
                                 0, 1))
1041
0
        {
1042
0
            psDefn->UOMAngleInDegrees = dfRadians / CPLAtof(SRS_UA_DEGREE_CONV);
1043
0
        }
1044
0
    }
1045
1046
    // Avoid later division by zero.
1047
0
    if (psDefn->UOMAngleInDegrees == 0)
1048
0
    {
1049
0
        CPLError(CE_Warning, CPLE_AppDefined,
1050
0
                 "Invalid value for GeogAngularUnitSizeGeoKey.");
1051
0
        psDefn->UOMAngleInDegrees = 1;
1052
0
    }
1053
1054
0
    dfSemiMajor = psDefn->SemiMajor;
1055
0
    if (dfSemiMajor == 0.0)
1056
0
    {
1057
0
        CPLFree(pszSpheroidName);
1058
0
        pszSpheroidName = CPLStrdup("unretrievable - using WGS84");
1059
0
        dfSemiMajor = SRS_WGS84_SEMIMAJOR;
1060
0
        dfInvFlattening = SRS_WGS84_INVFLATTENING;
1061
0
    }
1062
0
    else if (dfInvFlattening == 0.0 &&
1063
0
             ((psDefn->SemiMinor / psDefn->SemiMajor) < 0.99999999999999999 ||
1064
0
              (psDefn->SemiMinor / psDefn->SemiMajor) > 1.00000000000000001))
1065
0
    {
1066
0
        dfInvFlattening =
1067
0
            OSRCalcInvFlattening(psDefn->SemiMajor, psDefn->SemiMinor);
1068
1069
        /* Take official inverse flattening definition in the WGS84 case */
1070
0
        if (fabs(dfSemiMajor - SRS_WGS84_SEMIMAJOR) < 1e-10 &&
1071
0
            fabs(dfInvFlattening - SRS_WGS84_INVFLATTENING) < 1e-10)
1072
0
            dfInvFlattening = SRS_WGS84_INVFLATTENING;
1073
0
    }
1074
0
    if (!pszGeogName || strlen(pszGeogName) == 0)
1075
0
    {
1076
0
        CPLFree(pszGeogName);
1077
0
        pszGeogName = CPLStrdup(pszDatumName ? pszDatumName : "unknown");
1078
0
    }
1079
1080
0
    oSRS.SetGeogCS(pszGeogName, pszDatumName, pszSpheroidName, dfSemiMajor,
1081
0
                   dfInvFlattening, pszPMName, psDefn->PMLongToGreenwich,
1082
0
                   pszAngularUnits,
1083
0
                   psDefn->UOMAngleInDegrees * CPLAtof(SRS_UA_DEGREE_CONV));
1084
1085
0
    bool bGeog3DCRS = false;
1086
0
    bool bSetDatumEllipsoidCode = true;
1087
0
    bool bHasWarnedInconsistentGeogCRSEPSG = false;
1088
0
    {
1089
0
        const int nGCS = psDefn->GCS;
1090
0
        if (nGCS != KvUserDefined && nGCS > 0 &&
1091
0
            psDefn->Model != ModelTypeGeocentric)
1092
0
        {
1093
0
            OGRSpatialReference oSRSGeog;
1094
0
            const bool bGCSCodeValid =
1095
0
                oSRSGeog.importFromEPSG(nGCS) == OGRERR_NONE;
1096
1097
0
            const std::string osGTiffSRSSource =
1098
0
                CPLGetConfigOption("GTIFF_SRS_SOURCE", "");
1099
1100
            // GeoTIFF 1.0 might put a Geographic 3D code in GeodeticCRSGeoKey
1101
0
            bool bTryCompareToEPSG = oSRSGeog.GetAxesCount() == 2;
1102
1103
0
            if (psDefn->Datum != KvUserDefined)
1104
0
            {
1105
0
                char szCode[12];
1106
0
                snprintf(szCode, sizeof(szCode), "%d", psDefn->Datum);
1107
0
                auto ctx = static_cast<PJ_CONTEXT *>(
1108
0
                    GTIFGetPROJContext(hGTIF, true, nullptr));
1109
0
                auto datum = proj_create_from_database(
1110
0
                    ctx, "EPSG", szCode, PJ_CATEGORY_DATUM, 0, nullptr);
1111
0
                if (datum)
1112
0
                {
1113
0
                    if (proj_get_type(datum) ==
1114
0
                        PJ_TYPE_DYNAMIC_GEODETIC_REFERENCE_FRAME)
1115
0
                    {
1116
                        // Current PROJ versions will not manage to
1117
                        // consider a CRS with a regular datum and another one
1118
                        // with a dynamic datum as being equivalent.
1119
0
                        bTryCompareToEPSG = false;
1120
0
                    }
1121
0
                    proj_destroy(datum);
1122
0
                }
1123
0
            }
1124
1125
0
            if (bTryCompareToEPSG && !oSRSGeog.IsSameGeogCS(&oSRS) &&
1126
0
                osGTiffSRSSource.empty())
1127
0
            {
1128
                // See https://github.com/OSGeo/gdal/issues/5399
1129
                // where a file has inconsistent GeogSemiMinorAxisGeoKey /
1130
                // GeogInvFlatteningGeoKey values, which cause its datum to be
1131
                // considered as non-equivalent to the EPSG one.
1132
0
                CPLError(
1133
0
                    CE_Warning, CPLE_AppDefined,
1134
0
                    "The definition of geographic CRS EPSG:%d got from GeoTIFF "
1135
0
                    "keys "
1136
0
                    "is not the same as the one from the EPSG registry, "
1137
0
                    "which may cause issues during reprojection operations. "
1138
0
                    "Set GTIFF_SRS_SOURCE configuration option to EPSG to "
1139
0
                    "use official parameters (overriding the ones from GeoTIFF "
1140
0
                    "keys), "
1141
0
                    "or to GEOKEYS to use custom values from GeoTIFF keys "
1142
0
                    "and drop the EPSG code.",
1143
0
                    nGCS);
1144
0
                bHasWarnedInconsistentGeogCRSEPSG = true;
1145
0
            }
1146
0
            if (EQUAL(osGTiffSRSSource.c_str(), "EPSG"))
1147
0
            {
1148
0
                oSRS.CopyGeogCSFrom(&oSRSGeog);
1149
0
            }
1150
0
            else if (osGTiffSRSSource.empty() && oSRSGeog.IsDynamic() &&
1151
0
                     psDefn->Model == ModelTypeGeographic)
1152
0
            {
1153
                // We should perhaps be more careful and detect overrides
1154
                // of geokeys...
1155
0
                oSRS = oSRSGeog;
1156
0
                bSetDatumEllipsoidCode = false;
1157
0
            }
1158
0
            else if (bGCSCodeValid && osGTiffSRSSource.empty())
1159
0
            {
1160
0
                oSRS.SetAuthority("GEOGCS", "EPSG", nGCS);
1161
0
            }
1162
0
            else
1163
0
            {
1164
0
                bSetDatumEllipsoidCode = false;
1165
0
            }
1166
1167
0
            int nVertSRSCode = verticalCSType;
1168
0
            if (verticalDatum == 6030 && nGCS == 4326)  // DatumE_WGS84
1169
0
            {
1170
0
                nVertSRSCode = 4979;
1171
0
            }
1172
1173
            // Try to reconstruct a Geographic3D CRS from the
1174
            // GeodeticCRSGeoKey and the VerticalGeoKey, when they are
1175
            // consistent
1176
0
            if (nVertSRSCode > 0 && nVertSRSCode != KvUserDefined)
1177
0
            {
1178
0
                OGRSpatialReference oTmpVertSRS;
1179
0
                if (oSRSGeog.IsGeographic() && oSRSGeog.GetAxesCount() == 2 &&
1180
0
                    oTmpVertSRS.importFromEPSG(nVertSRSCode) == OGRERR_NONE &&
1181
0
                    oTmpVertSRS.IsGeographic() &&
1182
0
                    oTmpVertSRS.GetAxesCount() == 3)
1183
0
                {
1184
0
                    const char *pszTmpCode =
1185
0
                        oSRSGeog.GetAuthorityCode("GEOGCS|DATUM");
1186
0
                    const char *pszTmpVertCode =
1187
0
                        oTmpVertSRS.GetAuthorityCode("GEOGCS|DATUM");
1188
0
                    if (pszTmpCode && pszTmpVertCode &&
1189
0
                        atoi(pszTmpCode) == atoi(pszTmpVertCode))
1190
0
                    {
1191
0
                        verticalCSType = 0;
1192
0
                        verticalDatum = 0;
1193
0
                        verticalUnits = 0;
1194
0
                        oSRS.CopyGeogCSFrom(&oTmpVertSRS);
1195
0
                        bSetDatumEllipsoidCode = false;
1196
0
                        bGeog3DCRS = true;
1197
0
                    }
1198
0
                }
1199
0
            }
1200
0
        }
1201
0
    }
1202
0
    if (bSetDatumEllipsoidCode)
1203
0
    {
1204
0
        if (psDefn->Datum != KvUserDefined)
1205
0
            oSRS.SetAuthority("DATUM", "EPSG", psDefn->Datum);
1206
1207
0
        if (psDefn->Ellipsoid != KvUserDefined)
1208
0
            oSRS.SetAuthority("SPHEROID", "EPSG", psDefn->Ellipsoid);
1209
0
    }
1210
1211
0
    CPLFree(pszGeogName);
1212
0
    CPLFree(pszDatumName);
1213
0
    CPLFree(pszSpheroidName);
1214
0
    CPLFree(pszPMName);
1215
0
    CPLFree(pszAngularUnits);
1216
1217
    /* -------------------------------------------------------------------- */
1218
    /*      Set projection units if not yet done                            */
1219
    /* -------------------------------------------------------------------- */
1220
0
    if (psDefn->Model == ModelTypeProjected && !linearUnitIsSet)
1221
0
    {
1222
0
        char *pszUnitsName = nullptr;
1223
1224
0
        if (psDefn->UOMLength != KvUserDefined)
1225
0
        {
1226
0
#if LIBGEOTIFF_VERSION >= 1600
1227
0
            GTIFGetUOMLengthInfoEx(projContext,
1228
#else
1229
            GTIFGetUOMLengthInfo(
1230
#endif
1231
0
                                   psDefn->UOMLength, &pszUnitsName, nullptr);
1232
0
        }
1233
1234
0
        if (pszUnitsName != nullptr)
1235
0
        {
1236
0
            char szUOMLength[12];
1237
0
            snprintf(szUOMLength, sizeof(szUOMLength), "%d", psDefn->UOMLength);
1238
0
            oSRS.SetTargetLinearUnits(nullptr, pszUnitsName,
1239
0
                                      psDefn->UOMLengthInMeters, "EPSG",
1240
0
                                      szUOMLength);
1241
0
        }
1242
0
        else
1243
0
            oSRS.SetLinearUnits("unknown", psDefn->UOMLengthInMeters);
1244
1245
0
        GTIFFreeMemory(pszUnitsName);
1246
0
    }
1247
1248
    /* ==================================================================== */
1249
    /*      Try to import PROJCS from ProjectedCSTypeGeoKey if we           */
1250
    /*      have essentially only it. We could relax a bit the constraints  */
1251
    /*      but that should do for now. This may mask shortcomings in the   */
1252
    /*      libgeotiff GTIFGetDefn() function.                              */
1253
    /* ==================================================================== */
1254
0
    unsigned short tmp = 0;
1255
0
    bool bGotFromEPSG = false;
1256
0
    if (psDefn->Model == ModelTypeProjected && psDefn->PCS != KvUserDefined &&
1257
0
        GDALGTIFKeyGetSHORT(hGTIF, ProjectionGeoKey, &tmp, 0, 1) == 0 &&
1258
0
        GDALGTIFKeyGetSHORT(hGTIF, ProjCoordTransGeoKey, &tmp, 0, 1) == 0 &&
1259
0
        GDALGTIFKeyGetSHORT(hGTIF, GeographicTypeGeoKey, &tmp, 0, 1) == 0 &&
1260
0
        GDALGTIFKeyGetSHORT(hGTIF, GeogGeodeticDatumGeoKey, &tmp, 0, 1) == 0 &&
1261
0
        GDALGTIFKeyGetSHORT(hGTIF, GeogEllipsoidGeoKey, &tmp, 0, 1) == 0 &&
1262
0
        CPLTestBool(CPLGetConfigOption("GTIFF_IMPORT_FROM_EPSG", "YES")))
1263
0
    {
1264
        // Save error state as importFromEPSGA() will call CPLReset()
1265
0
        CPLErrorNum errNo = CPLGetLastErrorNo();
1266
0
        CPLErr eErr = CPLGetLastErrorType();
1267
0
        const char *pszTmp = CPLGetLastErrorMsg();
1268
0
        char *pszLastErrorMsg = CPLStrdup(pszTmp ? pszTmp : "");
1269
0
        CPLPushErrorHandler(CPLQuietErrorHandler);
1270
0
        OGRSpatialReference oSRSTmp;
1271
0
        OGRErr eImportErr = oSRSTmp.importFromEPSG(psDefn->PCS);
1272
0
        CPLPopErrorHandler();
1273
        // Restore error state
1274
0
        CPLErrorSetState(eErr, errNo, pszLastErrorMsg);
1275
0
        CPLFree(pszLastErrorMsg);
1276
0
        bGotFromEPSG = eImportErr == OGRERR_NONE;
1277
1278
0
        if (bGotFromEPSG)
1279
0
        {
1280
            // See #6210. In case there's an overridden linear units, take it
1281
            // into account
1282
0
            const char *pszUnitsName = nullptr;
1283
0
            double dfUOMLengthInMeters = oSRS.GetLinearUnits(&pszUnitsName);
1284
            // Non exact comparison, as there's a slight difference between
1285
            // the evaluation of US Survey foot hardcoded in geo_normalize.c to
1286
            // 12.0 / 39.37, and the corresponding value returned by
1287
            // PROJ >= 6.0.0 and <= 7.0.0 for EPSG:9003
1288
0
            if (fabs(dfUOMLengthInMeters - oSRSTmp.GetLinearUnits(nullptr)) >
1289
0
                1e-15 * dfUOMLengthInMeters)
1290
0
            {
1291
0
                CPLDebug("GTiff", "Modify EPSG:%d to have %s linear units...",
1292
0
                         psDefn->PCS, pszUnitsName ? pszUnitsName : "unknown");
1293
1294
0
                const char *pszUnitAuthorityCode =
1295
0
                    oSRS.GetAuthorityCode("PROJCS|UNIT");
1296
0
                const char *pszUnitAuthorityName =
1297
0
                    oSRS.GetAuthorityName("PROJCS|UNIT");
1298
1299
0
                if (pszUnitsName)
1300
0
                    oSRSTmp.SetLinearUnitsAndUpdateParameters(
1301
0
                        pszUnitsName, dfUOMLengthInMeters, pszUnitAuthorityCode,
1302
0
                        pszUnitAuthorityName);
1303
0
            }
1304
1305
0
            if (bGeog3DCRS)
1306
0
            {
1307
0
                oSRSTmp.CopyGeogCSFrom(&oSRS);
1308
0
                oSRSTmp.UpdateCoordinateSystemFromGeogCRS();
1309
0
            }
1310
0
            oSRS = std::move(oSRSTmp);
1311
0
        }
1312
0
    }
1313
1314
0
#if !defined(GEO_NORMALIZE_DISABLE_TOWGS84)
1315
0
    if (psDefn->TOWGS84Count > 0 && bGotFromEPSG &&
1316
0
        CPLTestBool(CPLGetConfigOption("OSR_STRIP_TOWGS84", "YES")))
1317
0
    {
1318
0
        CPLDebug("OSR", "TOWGS84 information has been removed. "
1319
0
                        "It can be kept by setting the OSR_STRIP_TOWGS84 "
1320
0
                        "configuration option to NO");
1321
0
    }
1322
0
    else if (psDefn->TOWGS84Count > 0 &&
1323
0
             (!bGotFromEPSG ||
1324
0
              !CPLTestBool(CPLGetConfigOption("OSR_STRIP_TOWGS84", "YES"))))
1325
0
    {
1326
0
        if (bGotFromEPSG)
1327
0
        {
1328
0
            double adfTOWGS84[7] = {0.0};
1329
0
            CPL_IGNORE_RET_VAL(oSRS.GetTOWGS84(adfTOWGS84));
1330
0
            bool bSame = true;
1331
0
            for (int i = 0; i < 7; i++)
1332
0
            {
1333
0
                if (fabs(adfTOWGS84[i] - psDefn->TOWGS84[i]) > 1e-5)
1334
0
                {
1335
0
                    bSame = false;
1336
0
                    break;
1337
0
                }
1338
0
            }
1339
0
            if (!bSame)
1340
0
            {
1341
0
                CPLDebug("GTiff",
1342
0
                         "Modify EPSG:%d to have "
1343
0
                         "TOWGS84=%f,%f,%f,%f,%f,%f,%f "
1344
0
                         "coming from GeogTOWGS84GeoKey, instead of "
1345
0
                         "%f,%f,%f,%f,%f,%f,%f coming from EPSG",
1346
0
                         psDefn->PCS, psDefn->TOWGS84[0], psDefn->TOWGS84[1],
1347
0
                         psDefn->TOWGS84[2], psDefn->TOWGS84[3],
1348
0
                         psDefn->TOWGS84[4], psDefn->TOWGS84[5],
1349
0
                         psDefn->TOWGS84[6], adfTOWGS84[0], adfTOWGS84[1],
1350
0
                         adfTOWGS84[2], adfTOWGS84[3], adfTOWGS84[4],
1351
0
                         adfTOWGS84[5], adfTOWGS84[6]);
1352
0
            }
1353
0
        }
1354
1355
0
        oSRS.SetTOWGS84(psDefn->TOWGS84[0], psDefn->TOWGS84[1],
1356
0
                        psDefn->TOWGS84[2], psDefn->TOWGS84[3],
1357
0
                        psDefn->TOWGS84[4], psDefn->TOWGS84[5],
1358
0
                        psDefn->TOWGS84[6]);
1359
0
    }
1360
0
#endif
1361
1362
    /* ==================================================================== */
1363
    /*      Handle projection parameters.                                   */
1364
    /* ==================================================================== */
1365
0
    if (psDefn->Model == ModelTypeProjected && !bGotFromEPSG)
1366
0
    {
1367
        /* --------------------------------------------------------------------
1368
         */
1369
        /*      Make a local copy of params, and convert back into the */
1370
        /*      angular units of the GEOGCS and the linear units of the */
1371
        /*      projection. */
1372
        /* --------------------------------------------------------------------
1373
         */
1374
0
        double adfParam[10] = {0.0};
1375
0
        int i = 0;  // Used after for.
1376
1377
0
        for (; i < std::min(10, psDefn->nParms); i++)
1378
0
            adfParam[i] = psDefn->ProjParm[i];
1379
1380
0
        for (; i < 10; i++)
1381
0
            adfParam[i] = 0.0;
1382
1383
#if LIBGEOTIFF_VERSION <= 1730
1384
        // libgeotiff <= 1.7.3 is unfortunately inconsistent. When it synthetizes the
1385
        // projection parameters from the EPSG ProjectedCRS code, it returns
1386
        // them normalized in degrees. But when it gets them from
1387
        // ProjCoordTransGeoKey and other Proj....GeoKey's it return them in
1388
        // a raw way, that is in the units of GeogAngularUnitSizeGeoKey
1389
        // The below oSRS.SetXXXXX() methods assume the angular projection
1390
        // parameters to be in degrees, so convert them to degrees in that later case.
1391
        // From GDAL 3.0 to 3.9.0, we didn't do that conversion...
1392
        // And all versions of GDAL <= 3.9.0 when writing those geokeys, wrote
1393
        // them as degrees, hence this GTIFF_READ_ANGULAR_PARAMS_IN_DEGREE
1394
        // config option that can be set to YES to avoid that conversion and
1395
        // assume that the angular parameters have been written as degree.
1396
        if (GDALGTIFKeyGetSHORT(hGTIF, ProjCoordTransGeoKey, &tmp, 0, 1) &&
1397
            !CPLTestBool(CPLGetConfigOption(
1398
                "GTIFF_READ_ANGULAR_PARAMS_IN_DEGREE", "NO")))
1399
        {
1400
            adfParam[0] *= psDefn->UOMAngleInDegrees;
1401
            adfParam[1] *= psDefn->UOMAngleInDegrees;
1402
            adfParam[2] *= psDefn->UOMAngleInDegrees;
1403
            adfParam[3] *= psDefn->UOMAngleInDegrees;
1404
        }
1405
#else
1406
        // If GTIFF_READ_ANGULAR_PARAMS_IN_DEGREE=YES (non-nominal case), undo
1407
        // the conversion to degrees, that has been done by libgeotiff > 1.7.3
1408
0
        if (GDALGTIFKeyGetSHORT(hGTIF, ProjCoordTransGeoKey, &tmp, 0, 1) &&
1409
0
            psDefn->UOMAngleInDegrees != 0 && psDefn->UOMAngleInDegrees != 1 &&
1410
0
            CPLTestBool(CPLGetConfigOption(
1411
0
                "GTIFF_READ_ANGULAR_PARAMS_IN_DEGREE", "NO")))
1412
0
        {
1413
0
            adfParam[0] /= psDefn->UOMAngleInDegrees;
1414
0
            adfParam[1] /= psDefn->UOMAngleInDegrees;
1415
0
            adfParam[2] /= psDefn->UOMAngleInDegrees;
1416
0
            adfParam[3] /= psDefn->UOMAngleInDegrees;
1417
0
        }
1418
0
#endif
1419
1420
        /* --------------------------------------------------------------------
1421
         */
1422
        /*      Translation the fundamental projection. */
1423
        /* --------------------------------------------------------------------
1424
         */
1425
0
        switch (psDefn->CTProjection)
1426
0
        {
1427
0
            case CT_TransverseMercator:
1428
0
                oSRS.SetTM(adfParam[0], adfParam[1], adfParam[4], adfParam[5],
1429
0
                           adfParam[6]);
1430
0
                break;
1431
1432
0
            case CT_TransvMercator_SouthOriented:
1433
0
                oSRS.SetTMSO(adfParam[0], adfParam[1], adfParam[4], adfParam[5],
1434
0
                             adfParam[6]);
1435
0
                break;
1436
1437
0
            case CT_Mercator:
1438
                // If a lat_ts was specified use 2SP, otherwise use 1SP.
1439
0
                if (psDefn->ProjParmId[2] == ProjStdParallel1GeoKey)
1440
0
                {
1441
0
                    if (psDefn->ProjParmId[4] == ProjScaleAtNatOriginGeoKey)
1442
0
                        CPLError(CE_Warning, CPLE_AppDefined,
1443
0
                                 "Mercator projection should not define "
1444
0
                                 "both StdParallel1 and ScaleAtNatOrigin.  "
1445
0
                                 "Using StdParallel1 and ignoring "
1446
0
                                 "ScaleAtNatOrigin.");
1447
0
                    oSRS.SetMercator2SP(adfParam[2], adfParam[0], adfParam[1],
1448
0
                                        adfParam[5], adfParam[6]);
1449
0
                }
1450
0
                else
1451
0
                    oSRS.SetMercator(adfParam[0], adfParam[1], adfParam[4],
1452
0
                                     adfParam[5], adfParam[6]);
1453
1454
                // Override hack for google mercator.
1455
0
                if (psDefn->Projection == 1024 || psDefn->Projection == 9841)
1456
0
                {
1457
0
                    oSRS.SetExtension(
1458
0
                        "PROJCS", "PROJ4",
1459
0
                        "+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 "
1460
0
                        "+lon_0=0.0 "
1461
0
                        "+x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null "
1462
0
                        "+wktext  +no_defs");  // TODO(schwehr): Why 2 spaces?
1463
0
                }
1464
0
                break;
1465
1466
0
            case CT_ObliqueStereographic:
1467
0
                oSRS.SetOS(adfParam[0], adfParam[1], adfParam[4], adfParam[5],
1468
0
                           adfParam[6]);
1469
0
                break;
1470
1471
0
            case CT_Stereographic:
1472
0
                oSRS.SetStereographic(adfParam[0], adfParam[1], adfParam[4],
1473
0
                                      adfParam[5], adfParam[6]);
1474
0
                break;
1475
1476
0
            case CT_ObliqueMercator:  // Hotine.
1477
0
                oSRS.SetHOM(adfParam[0], adfParam[1], adfParam[2], adfParam[3],
1478
0
                            adfParam[4], adfParam[5], adfParam[6]);
1479
0
                break;
1480
1481
0
            case CT_HotineObliqueMercatorAzimuthCenter:
1482
0
                oSRS.SetHOMAC(adfParam[0], adfParam[1], adfParam[2],
1483
0
                              adfParam[3], adfParam[4], adfParam[5],
1484
0
                              adfParam[6]);
1485
0
                break;
1486
1487
0
            case CT_ObliqueMercator_Laborde:
1488
0
                oSRS.SetLOM(adfParam[0], adfParam[1], adfParam[2], adfParam[4],
1489
0
                            adfParam[5], adfParam[6]);
1490
0
                break;
1491
1492
0
            case CT_EquidistantConic:
1493
0
                oSRS.SetEC(adfParam[0], adfParam[1], adfParam[2], adfParam[3],
1494
0
                           adfParam[5], adfParam[6]);
1495
0
                break;
1496
1497
0
            case CT_CassiniSoldner:
1498
0
                oSRS.SetCS(adfParam[0], adfParam[1], adfParam[5], adfParam[6]);
1499
0
                break;
1500
1501
0
            case CT_Polyconic:
1502
0
                oSRS.SetPolyconic(adfParam[0], adfParam[1], adfParam[5],
1503
0
                                  adfParam[6]);
1504
0
                break;
1505
1506
0
            case CT_AzimuthalEquidistant:
1507
0
                oSRS.SetAE(adfParam[0], adfParam[1], adfParam[5], adfParam[6]);
1508
0
                break;
1509
1510
0
            case CT_MillerCylindrical:
1511
0
                oSRS.SetMC(adfParam[0], adfParam[1], adfParam[5], adfParam[6]);
1512
0
                break;
1513
1514
0
            case CT_Equirectangular:
1515
0
                oSRS.SetEquirectangular2(adfParam[0], adfParam[1], adfParam[2],
1516
0
                                         adfParam[5], adfParam[6]);
1517
0
                break;
1518
1519
0
            case CT_Gnomonic:
1520
0
                oSRS.SetGnomonic(adfParam[0], adfParam[1], adfParam[5],
1521
0
                                 adfParam[6]);
1522
0
                break;
1523
1524
0
            case CT_LambertAzimEqualArea:
1525
0
                oSRS.SetLAEA(adfParam[0], adfParam[1], adfParam[5],
1526
0
                             adfParam[6]);
1527
0
                break;
1528
1529
0
            case CT_Orthographic:
1530
0
                oSRS.SetOrthographic(adfParam[0], adfParam[1], adfParam[5],
1531
0
                                     adfParam[6]);
1532
0
                break;
1533
1534
0
            case CT_Robinson:
1535
0
                oSRS.SetRobinson(adfParam[1], adfParam[5], adfParam[6]);
1536
0
                break;
1537
1538
0
            case CT_Sinusoidal:
1539
0
                oSRS.SetSinusoidal(adfParam[1], adfParam[5], adfParam[6]);
1540
0
                break;
1541
1542
0
            case CT_VanDerGrinten:
1543
0
                oSRS.SetVDG(adfParam[1], adfParam[5], adfParam[6]);
1544
0
                break;
1545
1546
0
            case CT_PolarStereographic:
1547
0
                oSRS.SetPS(adfParam[0], adfParam[1], adfParam[4], adfParam[5],
1548
0
                           adfParam[6]);
1549
0
                break;
1550
1551
0
            case CT_LambertConfConic_2SP:
1552
0
                oSRS.SetLCC(adfParam[2], adfParam[3], adfParam[0], adfParam[1],
1553
0
                            adfParam[5], adfParam[6]);
1554
0
                break;
1555
1556
0
            case CT_LambertConfConic_1SP:
1557
0
                oSRS.SetLCC1SP(adfParam[0], adfParam[1], adfParam[4],
1558
0
                               adfParam[5], adfParam[6]);
1559
0
                break;
1560
1561
0
            case CT_AlbersEqualArea:
1562
0
                oSRS.SetACEA(adfParam[0], adfParam[1], adfParam[2], adfParam[3],
1563
0
                             adfParam[5], adfParam[6]);
1564
0
                break;
1565
1566
0
            case CT_NewZealandMapGrid:
1567
0
                oSRS.SetNZMG(adfParam[0], adfParam[1], adfParam[5],
1568
0
                             adfParam[6]);
1569
0
                break;
1570
1571
0
            case CT_CylindricalEqualArea:
1572
0
                oSRS.SetCEA(adfParam[0], adfParam[1], adfParam[5], adfParam[6]);
1573
0
                break;
1574
0
            default:
1575
0
                if (oSRS.IsProjected())
1576
0
                {
1577
0
                    const char *pszName = oSRS.GetName();
1578
0
                    std::string osName(pszName ? pszName : "unnamed");
1579
0
                    oSRS.Clear();
1580
0
                    oSRS.SetLocalCS(osName.c_str());
1581
0
                }
1582
0
                break;
1583
0
        }
1584
0
    }
1585
1586
0
    if (psDefn->Model == ModelTypeProjected && psDefn->PCS != KvUserDefined &&
1587
0
        !bGotFromEPSG)
1588
0
    {
1589
0
        OGRSpatialReference oSRSTest(oSRS);
1590
0
        OGRSpatialReference oSRSTmp;
1591
1592
0
        const bool bPCSCodeValid =
1593
0
            oSRSTmp.importFromEPSG(psDefn->PCS) == OGRERR_NONE;
1594
0
        oSRSTmp.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
1595
1596
        // Force axis to avoid issues with SRS with northing, easting order
1597
0
        oSRSTest.SetAxes(nullptr, "X", OAO_East, "Y", OAO_North);
1598
0
        oSRSTmp.SetAxes(nullptr, "X", OAO_East, "Y", OAO_North);
1599
1600
0
        const std::string osGTiffSRSSource =
1601
0
            CPLGetConfigOption("GTIFF_SRS_SOURCE", "");
1602
0
        const char *const apszOptions[] = {
1603
0
            "IGNORE_DATA_AXIS_TO_SRS_AXIS_MAPPING=YES", nullptr};
1604
0
        if (!bHasWarnedInconsistentGeogCRSEPSG &&
1605
0
            !oSRSTmp.IsSame(&oSRS, apszOptions) &&
1606
0
            EQUAL(osGTiffSRSSource.c_str(), ""))
1607
0
        {
1608
            // See https://github.com/OSGeo/gdal/issues/5399
1609
            // where a file has inconsistent GeogSemiMinorAxisGeoKey /
1610
            // GeogInvFlatteningGeoKey values, which cause its datum to be
1611
            // considered as non-equivalent to the EPSG one.
1612
0
            CPLError(
1613
0
                CE_Warning, CPLE_AppDefined,
1614
0
                "The definition of projected CRS EPSG:%d got from GeoTIFF keys "
1615
0
                "is not the same as the one from the EPSG registry, "
1616
0
                "which may cause issues during reprojection operations. "
1617
0
                "Set GTIFF_SRS_SOURCE configuration option to EPSG to "
1618
0
                "use official parameters (overriding the ones from GeoTIFF "
1619
0
                "keys), "
1620
0
                "or to GEOKEYS to use custom values from GeoTIFF keys "
1621
0
                "and drop the EPSG code.",
1622
0
                psDefn->PCS);
1623
0
        }
1624
0
        if (EQUAL(osGTiffSRSSource.c_str(), "EPSG"))
1625
0
        {
1626
0
            oSRS = std::move(oSRSTmp);
1627
0
        }
1628
0
        else if (bPCSCodeValid && EQUAL(osGTiffSRSSource.c_str(), ""))
1629
0
        {
1630
0
            oSRS.SetAuthority(nullptr, "EPSG", psDefn->PCS);
1631
0
        }
1632
0
    }
1633
1634
0
    if (oSRS.IsProjected() && oSRS.GetAxesCount() == 2)
1635
0
    {
1636
0
        const char *pszProjCRSName = oSRS.GetAttrValue("PROJCS");
1637
0
        if (pszProjCRSName)
1638
0
        {
1639
            // Hack to be able to read properly what we have written for
1640
            // ESRI:102113 (ESRI ancient WebMercator).
1641
0
            if (EQUAL(pszProjCRSName, "WGS_1984_Web_Mercator"))
1642
0
                oSRS.SetFromUserInput("ESRI:102113");
1643
            // And for EPSG:900913
1644
0
            else if (EQUAL(pszProjCRSName, "Google Maps Global Mercator"))
1645
0
                oSRS.importFromEPSG(900913);
1646
0
            else if (strchr(pszProjCRSName, '_'))
1647
0
            {
1648
                // Morph from potential ESRI name
1649
0
                char *pszOfficialName = GTIFGetEPSGOfficialName(
1650
0
                    hGTIF, PJ_TYPE_PROJECTED_CRS, pszProjCRSName);
1651
0
                if (pszOfficialName)
1652
0
                {
1653
0
                    oSRS.SetProjCS(pszOfficialName);
1654
0
                    CPLFree(pszOfficialName);
1655
0
                }
1656
0
            }
1657
0
        }
1658
0
    }
1659
1660
    /* ==================================================================== */
1661
    /*      Handle vertical coordinate system information if we have it.    */
1662
    /* ==================================================================== */
1663
0
    bool bNeedManualVertCS = false;
1664
0
    char citation[2048] = {'\0'};
1665
1666
    // See https://github.com/OSGeo/gdal/pull/4197
1667
0
    if (verticalCSType > KvUserDefined || verticalDatum > KvUserDefined ||
1668
0
        verticalUnits > KvUserDefined)
1669
0
    {
1670
0
        CPLError(CE_Warning, CPLE_AppDefined,
1671
0
                 "At least one of VerticalCSTypeGeoKey, VerticalDatumGeoKey or "
1672
0
                 "VerticalUnitsGeoKey has a value in the private user range. "
1673
0
                 "Ignoring vertical information.");
1674
0
        verticalCSType = 0;
1675
0
        verticalDatum = 0;
1676
0
        verticalUnits = 0;
1677
0
    }
1678
1679
0
    if ((verticalCSType != 0 || verticalDatum != 0 || verticalUnits != 0) &&
1680
0
        (oSRS.IsGeographic() || oSRS.IsProjected() || oSRS.IsLocal()))
1681
0
    {
1682
0
        std::string osVertCRSName;
1683
0
        if (GDALGTIFKeyGetASCII(hGTIF, VerticalCitationGeoKey, citation,
1684
0
                                sizeof(citation)))
1685
0
        {
1686
0
            if (STARTS_WITH_CI(citation, "VCS Name = "))
1687
0
            {
1688
0
                memmove(citation, citation + strlen("VCS Name = "),
1689
0
                        strlen(citation + strlen("VCS Name = ")) + 1);
1690
0
                char *pszPipeChar = strchr(citation, '|');
1691
0
                if (pszPipeChar)
1692
0
                    *pszPipeChar = '\0';
1693
0
                osVertCRSName = citation;
1694
0
            }
1695
0
        }
1696
1697
0
        OGRSpatialReference oVertSRS;
1698
0
        bool bCanBuildCompoundCRS = oSRS.GetRoot() != nullptr;
1699
0
        if (verticalCSType != KvUserDefined && verticalCSType > 0)
1700
0
        {
1701
0
            if (!(oVertSRS.importFromEPSG(verticalCSType) == OGRERR_NONE &&
1702
0
                  oVertSRS.IsVertical()))
1703
0
            {
1704
0
                bCanBuildCompoundCRS = false;
1705
0
            }
1706
0
            else
1707
0
            {
1708
0
                osVertCRSName = oVertSRS.GetName();
1709
0
            }
1710
0
        }
1711
0
        if (osVertCRSName.empty())
1712
0
            osVertCRSName = "unknown";
1713
1714
0
        if (bCanBuildCompoundCRS)
1715
0
        {
1716
0
            const bool bHorizontalHasCode =
1717
0
                oSRS.GetAuthorityCode(nullptr) != nullptr;
1718
0
            const char *pszHorizontalName = oSRS.GetName();
1719
0
            const std::string osHorizontalName(
1720
0
                pszHorizontalName ? pszHorizontalName : "unnamed");
1721
            /* --------------------------------------------------------------------
1722
             */
1723
            /*      Promote to being a compound coordinate system. */
1724
            /* --------------------------------------------------------------------
1725
             */
1726
0
            OGR_SRSNode *poOldRoot = oSRS.GetRoot()->Clone();
1727
1728
0
            oSRS.Clear();
1729
1730
            /* --------------------------------------------------------------------
1731
             */
1732
            /*      Set COMPD_CS name. */
1733
            /* --------------------------------------------------------------------
1734
             */
1735
0
            char szCTString[512];
1736
0
            szCTString[0] = '\0';
1737
0
            if (GDALGTIFKeyGetASCII(hGTIF, GTCitationGeoKey, szCTString,
1738
0
                                    sizeof(szCTString)) &&
1739
0
                strstr(szCTString, " = ") == nullptr)
1740
0
            {
1741
0
                oSRS.SetNode("COMPD_CS", szCTString);
1742
0
            }
1743
0
            else
1744
0
            {
1745
0
                oSRS.SetNode(
1746
0
                    "COMPD_CS",
1747
0
                    (osHorizontalName + " + " + osVertCRSName).c_str());
1748
0
            }
1749
1750
0
            oSRS.GetRoot()->AddChild(poOldRoot);
1751
1752
            /* --------------------------------------------------------------------
1753
             */
1754
            /*      If we have the vertical cs, try to look it up, and use the
1755
             */
1756
            /*      definition provided by that. */
1757
            /* --------------------------------------------------------------------
1758
             */
1759
0
            bNeedManualVertCS = true;
1760
1761
0
            if (!oVertSRS.IsEmpty())
1762
0
            {
1763
0
                oSRS.GetRoot()->AddChild(oVertSRS.GetRoot()->Clone());
1764
0
                bNeedManualVertCS = false;
1765
1766
                // GeoTIFF doesn't store EPSG code of CompoundCRS, so
1767
                // if we have an EPSG code for both the horizontal and vertical
1768
                // parts, check if there's a known CompoundCRS associating
1769
                // both
1770
0
                if (bHorizontalHasCode && verticalCSType != KvUserDefined &&
1771
0
                    verticalCSType > 0)
1772
0
                {
1773
0
                    const auto *poSRSMatch = oSRS.FindBestMatch(100);
1774
0
                    if (poSRSMatch)
1775
0
                        oSRS = *poSRSMatch;
1776
0
                    delete poSRSMatch;
1777
0
                }
1778
0
            }
1779
0
        }
1780
0
    }
1781
1782
    /* -------------------------------------------------------------------- */
1783
    /*      Collect some information from the VerticalCS if not provided    */
1784
    /*      via geokeys.                                                    */
1785
    /* -------------------------------------------------------------------- */
1786
0
    if (bNeedManualVertCS)
1787
0
    {
1788
0
        FillCompoundCRSWithManualVertCS(hGTIF, oSRS, citation, verticalDatum,
1789
0
                                        verticalUnits);
1790
0
    }
1791
1792
    // Hack for tiff_read.py:test_tiff_grads so as to normalize angular
1793
    // parameters to grad
1794
0
    if (psDefn->UOMAngleInDegrees != 1.0)
1795
0
    {
1796
0
        char *pszWKT = nullptr;
1797
0
        const char *const apszOptions[] = {
1798
0
            "FORMAT=WKT1", "ADD_TOWGS84_ON_EXPORT_TO_WKT1=NO", nullptr};
1799
0
        if (oSRS.exportToWkt(&pszWKT, apszOptions) == OGRERR_NONE)
1800
0
        {
1801
0
            const char *const apszOptionsImport[] = {
1802
0
#if PROJ_AT_LEAST_VERSION(9, 1, 0)
1803
0
                "UNSET_IDENTIFIERS_IF_INCOMPATIBLE_DEF=NO",
1804
0
#endif
1805
0
                nullptr
1806
0
            };
1807
0
            oSRS.importFromWkt(pszWKT, apszOptionsImport);
1808
0
        }
1809
0
        CPLFree(pszWKT);
1810
0
    }
1811
1812
0
    oSRS.StripTOWGS84IfKnownDatumAndAllowed();
1813
1814
0
    double dfCoordinateEpoch = 0.0;
1815
0
    if (GDALGTIFKeyGetDOUBLE(hGTIF, CoordinateEpochGeoKey, &dfCoordinateEpoch,
1816
0
                             0, 1))
1817
0
    {
1818
0
        oSRS.SetCoordinateEpoch(dfCoordinateEpoch);
1819
0
    }
1820
1821
0
    return OGRSpatialReference::ToHandle(oSRS.Clone());
1822
0
}
1823
1824
/************************************************************************/
1825
/*                          GTIFGetOGISDefn()                           */
1826
/************************************************************************/
1827
1828
char *GTIFGetOGISDefn(GTIF *hGTIF, GTIFDefn *psDefn)
1829
0
{
1830
0
    OGRSpatialReferenceH hSRS = GTIFGetOGISDefnAsOSR(hGTIF, psDefn);
1831
1832
0
    char *pszWKT = nullptr;
1833
0
    if (hSRS && OGRSpatialReference::FromHandle(hSRS)->exportToWkt(&pszWKT) ==
1834
0
                    OGRERR_NONE)
1835
0
    {
1836
0
        OSRDestroySpatialReference(hSRS);
1837
0
        return pszWKT;
1838
0
    }
1839
0
    CPLFree(pszWKT);
1840
0
    OSRDestroySpatialReference(hSRS);
1841
1842
0
    return nullptr;
1843
0
}
1844
1845
/************************************************************************/
1846
/*                     OGCDatumName2EPSGDatumCode()                     */
1847
/************************************************************************/
1848
1849
static int OGCDatumName2EPSGDatumCode(GTIF *psGTIF, const char *pszOGCName)
1850
1851
0
{
1852
0
    int nReturn = KvUserDefined;
1853
1854
    /* -------------------------------------------------------------------- */
1855
    /*      Do we know it as a built in?                                    */
1856
    /* -------------------------------------------------------------------- */
1857
0
    if (EQUAL(pszOGCName, "NAD27") ||
1858
0
        EQUAL(pszOGCName, "North_American_Datum_1927"))
1859
0
        return Datum_North_American_Datum_1927;
1860
0
    else if (EQUAL(pszOGCName, "NAD83") ||
1861
0
             EQUAL(pszOGCName, "North_American_Datum_1983"))
1862
0
        return Datum_North_American_Datum_1983;
1863
0
    else if (EQUAL(pszOGCName, "WGS84") || EQUAL(pszOGCName, "WGS_1984") ||
1864
0
             EQUAL(pszOGCName, "WGS 84"))
1865
0
        return Datum_WGS84;
1866
0
    else if (EQUAL(pszOGCName, "WGS72") || EQUAL(pszOGCName, "WGS_1972"))
1867
0
        return Datum_WGS72;
1868
1869
    /* Search in database */
1870
0
    auto ctx =
1871
0
        static_cast<PJ_CONTEXT *>(GTIFGetPROJContext(psGTIF, true, nullptr));
1872
0
    const PJ_TYPE searchType = PJ_TYPE_GEODETIC_REFERENCE_FRAME;
1873
0
    auto list = proj_create_from_name(ctx, "EPSG", pszOGCName, &searchType, 1,
1874
0
                                      true, /* approximate match */
1875
0
                                      10, nullptr);
1876
0
    if (list)
1877
0
    {
1878
0
        const auto listSize = proj_list_get_count(list);
1879
0
        for (int i = 0; nReturn == KvUserDefined && i < listSize; i++)
1880
0
        {
1881
0
            auto datum = proj_list_get(ctx, list, i);
1882
0
            if (datum)
1883
0
            {
1884
0
                const char *pszDatumName = proj_get_name(datum);
1885
0
                if (pszDatumName)
1886
0
                {
1887
0
                    char *pszTmp = CPLStrdup(pszDatumName);
1888
0
                    WKTMassageDatum(&pszTmp);
1889
0
                    if (EQUAL(pszTmp, pszOGCName))
1890
0
                    {
1891
0
                        const char *pszCode = proj_get_id_code(datum, 0);
1892
0
                        if (pszCode)
1893
0
                        {
1894
0
                            nReturn = atoi(pszCode);
1895
0
                        }
1896
0
                    }
1897
0
                    CPLFree(pszTmp);
1898
0
                }
1899
0
            }
1900
0
            proj_destroy(datum);
1901
0
        }
1902
0
        proj_list_destroy(list);
1903
0
    }
1904
1905
0
    return nReturn;
1906
0
}
1907
1908
/************************************************************************/
1909
/*                        GTIFSetFromOGISDefn()                         */
1910
/*                                                                      */
1911
/*      Write GeoTIFF projection tags from an OGC WKT definition.       */
1912
/************************************************************************/
1913
1914
int GTIFSetFromOGISDefn(GTIF *psGTIF, const char *pszOGCWKT)
1915
1916
0
{
1917
    /* -------------------------------------------------------------------- */
1918
    /*      Create an OGRSpatialReference object corresponding to the       */
1919
    /*      string.                                                         */
1920
    /* -------------------------------------------------------------------- */
1921
1922
0
    OGRSpatialReference oSRS;
1923
0
    if (oSRS.importFromWkt(pszOGCWKT) != OGRERR_NONE)
1924
0
    {
1925
0
        return FALSE;
1926
0
    }
1927
0
    return GTIFSetFromOGISDefnEx(psGTIF, OGRSpatialReference::ToHandle(&oSRS),
1928
0
                                 GEOTIFF_KEYS_STANDARD, GEOTIFF_VERSION_1_0);
1929
0
}
1930
1931
int GTIFSetFromOGISDefnEx(GTIF *psGTIF, OGRSpatialReferenceH hSRS,
1932
                          GTIFFKeysFlavorEnum eFlavor,
1933
                          GeoTIFFVersionEnum eVersion)
1934
0
{
1935
0
    std::map<geokey_t, std::string> oMapAsciiKeys;
1936
1937
0
    GTIFKeySet(psGTIF, GTRasterTypeGeoKey, TYPE_SHORT, 1, RasterPixelIsArea);
1938
1939
0
    const OGRSpatialReference *poSRS = OGRSpatialReference::FromHandle(hSRS);
1940
1941
    /* -------------------------------------------------------------------- */
1942
    /*      Set version number.                                             */
1943
    /* -------------------------------------------------------------------- */
1944
0
    if (eVersion == GEOTIFF_VERSION_AUTO)
1945
0
    {
1946
0
        if (poSRS->IsCompound() ||
1947
0
            (poSRS->IsGeographic() && poSRS->GetAxesCount() == 3))
1948
0
        {
1949
0
            eVersion = GEOTIFF_VERSION_1_1;
1950
0
        }
1951
0
        else
1952
0
        {
1953
0
            eVersion = GEOTIFF_VERSION_1_0;
1954
0
        }
1955
0
    }
1956
0
    CPLAssert(eVersion == GEOTIFF_VERSION_1_0 ||
1957
0
              eVersion == GEOTIFF_VERSION_1_1);
1958
0
    if (eVersion >= GEOTIFF_VERSION_1_1)
1959
0
    {
1960
0
#if LIBGEOTIFF_VERSION >= 1600
1961
0
        GTIFSetVersionNumbers(psGTIF, GEOTIFF_SPEC_1_1_VERSION,
1962
0
                              GEOTIFF_SPEC_1_1_KEY_REVISION,
1963
0
                              GEOTIFF_SPEC_1_1_MINOR_REVISION);
1964
#else
1965
        CPLError(CE_Warning, CPLE_AppDefined,
1966
                 "Setting GeoTIFF 1.1 requires libgeotiff >= 1.6. Key values "
1967
                 "will be written as GeoTIFF 1.1, but the version number "
1968
                 "will be seen as 1.0, which might confuse GeoTIFF readers");
1969
#endif
1970
0
    }
1971
1972
    /* -------------------------------------------------------------------- */
1973
    /*      Get the ellipsoid definition.                                   */
1974
    /* -------------------------------------------------------------------- */
1975
0
    short nSpheroid = KvUserDefined;
1976
1977
0
    if (poSRS->GetAuthorityName("PROJCS|GEOGCS|DATUM|SPHEROID") != nullptr &&
1978
0
        EQUAL(poSRS->GetAuthorityName("PROJCS|GEOGCS|DATUM|SPHEROID"), "EPSG"))
1979
0
    {
1980
0
        nSpheroid = static_cast<short>(
1981
0
            atoi(poSRS->GetAuthorityCode("PROJCS|GEOGCS|DATUM|SPHEROID")));
1982
0
    }
1983
0
    else if (poSRS->GetAuthorityName("GEOGCS|DATUM|SPHEROID") != nullptr &&
1984
0
             EQUAL(poSRS->GetAuthorityName("GEOGCS|DATUM|SPHEROID"), "EPSG"))
1985
0
    {
1986
0
        nSpheroid = static_cast<short>(
1987
0
            atoi(poSRS->GetAuthorityCode("GEOGCS|DATUM|SPHEROID")));
1988
0
    }
1989
1990
0
    OGRErr eErr = OGRERR_NONE;
1991
0
    double dfSemiMajor = 0;
1992
0
    double dfInvFlattening = 0;
1993
0
    bool bHasEllipsoid = false;
1994
0
    if (!poSRS->IsLocal())
1995
0
    {
1996
0
        bHasEllipsoid = true;
1997
0
        if (poSRS->IsCompound())
1998
0
        {
1999
0
            OGRSpatialReference oSRSTmp(*poSRS);
2000
0
            oSRSTmp.StripVertical();
2001
0
            bHasEllipsoid = CPL_TO_BOOL(!oSRSTmp.IsLocal());
2002
0
        }
2003
0
        if (bHasEllipsoid)
2004
0
        {
2005
0
            dfSemiMajor = poSRS->GetSemiMajor(&eErr);
2006
0
            dfInvFlattening = poSRS->GetInvFlattening(&eErr);
2007
0
            if (eErr != OGRERR_NONE)
2008
0
            {
2009
0
                dfSemiMajor = 0.0;
2010
0
                dfInvFlattening = 0.0;
2011
0
            }
2012
0
        }
2013
0
    }
2014
2015
    /* -------------------------------------------------------------------- */
2016
    /*      Get the Datum so we can special case a few PCS codes.           */
2017
    /* -------------------------------------------------------------------- */
2018
0
    int nDatum = KvUserDefined;
2019
2020
0
    if (poSRS->GetAuthorityName("PROJCS|GEOGCS|DATUM") != nullptr &&
2021
0
        EQUAL(poSRS->GetAuthorityName("PROJCS|GEOGCS|DATUM"), "EPSG"))
2022
0
        nDatum = atoi(poSRS->GetAuthorityCode("PROJCS|GEOGCS|DATUM"));
2023
0
    else if (poSRS->GetAuthorityName("GEOGCS|DATUM") != nullptr &&
2024
0
             EQUAL(poSRS->GetAuthorityName("GEOGCS|DATUM"), "EPSG"))
2025
0
        nDatum = atoi(poSRS->GetAuthorityCode("GEOGCS|DATUM"));
2026
0
    else if (poSRS->GetAttrValue("DATUM") != nullptr)
2027
0
        nDatum =
2028
0
            OGCDatumName2EPSGDatumCode(psGTIF, poSRS->GetAttrValue("DATUM"));
2029
2030
    /* -------------------------------------------------------------------- */
2031
    /*      Get the GCS if possible.                                        */
2032
    /* -------------------------------------------------------------------- */
2033
0
    int nGCS = KvUserDefined;
2034
2035
0
    if (poSRS->GetAuthorityName("PROJCS|GEOGCS") != nullptr &&
2036
0
        EQUAL(poSRS->GetAuthorityName("PROJCS|GEOGCS"), "EPSG"))
2037
0
        nGCS = atoi(poSRS->GetAuthorityCode("PROJCS|GEOGCS"));
2038
0
    else if (poSRS->GetAuthorityName("GEOGCS") != nullptr &&
2039
0
             EQUAL(poSRS->GetAuthorityName("GEOGCS"), "EPSG"))
2040
0
        nGCS = atoi(poSRS->GetAuthorityCode("GEOGCS"));
2041
2042
0
    int nVerticalCSKeyValue = 0;
2043
0
    bool hasEllipsoidHeight = !poSRS->IsCompound() && poSRS->IsGeographic() &&
2044
0
                              poSRS->GetAxesCount() == 3;
2045
0
    if (nGCS == 4937 && eVersion >= GEOTIFF_VERSION_1_1)
2046
0
    {
2047
        // Workaround a bug of PROJ 6.3.0
2048
        // See https://github.com/OSGeo/PROJ/pull/1880
2049
        // EPSG:4937 = ETRS89 3D
2050
0
        hasEllipsoidHeight = true;
2051
0
        nVerticalCSKeyValue = nGCS;
2052
0
        nGCS = 4258;  // ETRS89 2D
2053
0
    }
2054
0
    else if (nGCS != KvUserDefined)
2055
0
    {
2056
0
        OGRSpatialReference oGeogCRS;
2057
0
        if (oGeogCRS.importFromEPSG(nGCS) == OGRERR_NONE &&
2058
0
            oGeogCRS.IsGeographic() && oGeogCRS.GetAxesCount() == 3)
2059
0
        {
2060
0
            hasEllipsoidHeight = true;
2061
0
            if (eVersion >= GEOTIFF_VERSION_1_1)
2062
0
            {
2063
0
                const auto candidate_nVerticalCSKeyValue = nGCS;
2064
0
                nGCS = KvUserDefined;
2065
2066
                // In case of a geographic 3D CRS, find the corresponding
2067
                // geographic 2D CRS
2068
0
                auto ctx = static_cast<PJ_CONTEXT *>(
2069
0
                    GTIFGetPROJContext(psGTIF, true, nullptr));
2070
0
                const auto type = PJ_TYPE_GEOGRAPHIC_2D_CRS;
2071
0
                auto list = proj_create_from_name(ctx, "EPSG",
2072
0
                                                  oGeogCRS.GetName(), &type, 1,
2073
0
                                                  false,  // exact match
2074
0
                                                  1,  // result set limit size,
2075
0
                                                  nullptr);
2076
0
                if (list && proj_list_get_count(list) == 1)
2077
0
                {
2078
0
                    auto crs2D = proj_list_get(ctx, list, 0);
2079
0
                    if (crs2D)
2080
0
                    {
2081
0
                        const char *pszCode = proj_get_id_code(crs2D, 0);
2082
0
                        if (pszCode)
2083
0
                        {
2084
0
                            nVerticalCSKeyValue = candidate_nVerticalCSKeyValue;
2085
0
                            nGCS = atoi(pszCode);
2086
0
                        }
2087
0
                        proj_destroy(crs2D);
2088
0
                    }
2089
0
                }
2090
0
                proj_list_destroy(list);
2091
0
            }
2092
0
        }
2093
0
    }
2094
2095
    // Deprecated way of encoding ellipsoidal height
2096
0
    if (hasEllipsoidHeight && nVerticalCSKeyValue == 0)
2097
0
    {
2098
0
        if (nGCS == 4979 || nDatum == 6326 || nSpheroid == 7030)
2099
0
        {
2100
0
            nVerticalCSKeyValue = 5030;  // WGS_84_ellipsoid
2101
0
            if (nGCS == 4979 || nDatum == 6326)
2102
0
            {
2103
0
                nGCS = 4326;
2104
0
            }
2105
0
        }
2106
0
        else if (nDatum >= 6001 && nDatum <= 6033)
2107
0
        {
2108
0
            nVerticalCSKeyValue = nDatum - 1000;
2109
0
        }
2110
0
        else if (nSpheroid >= 7001 && nSpheroid <= 7033)
2111
0
        {
2112
0
            nVerticalCSKeyValue = nSpheroid - 2000;
2113
0
        }
2114
0
    }
2115
2116
0
    if (nGCS > 32767)
2117
0
        nGCS = KvUserDefined;
2118
2119
    /* -------------------------------------------------------------------- */
2120
    /*      Get the linear units.                                           */
2121
    /* -------------------------------------------------------------------- */
2122
0
    const char *pszLinearUOMNameTmp = nullptr;
2123
0
    const double dfLinearUOM = poSRS->GetLinearUnits(&pszLinearUOMNameTmp);
2124
0
    const std::string osLinearUOMName(pszLinearUOMNameTmp ? pszLinearUOMNameTmp
2125
0
                                                          : "");
2126
0
    int nUOMLengthCode = 9001;  // Meters.
2127
2128
0
    if (poSRS->GetAuthorityName("PROJCS|UNIT") != nullptr &&
2129
0
        EQUAL(poSRS->GetAuthorityName("PROJCS|UNIT"), "EPSG") &&
2130
0
        poSRS->GetAttrNode("PROJCS|UNIT") != poSRS->GetAttrNode("GEOGCS|UNIT"))
2131
0
        nUOMLengthCode = atoi(poSRS->GetAuthorityCode("PROJCS|UNIT"));
2132
0
    else if (EQUAL(osLinearUOMName.c_str(), SRS_UL_FOOT) ||
2133
0
             fabs(dfLinearUOM - GTIFAtof(SRS_UL_FOOT_CONV)) < 0.0000001)
2134
0
        nUOMLengthCode = 9002;  // International foot.
2135
0
    else if (EQUAL(osLinearUOMName.c_str(), SRS_UL_US_FOOT) ||
2136
0
             std::abs(dfLinearUOM - GTIFAtof(SRS_UL_US_FOOT_CONV)) < 0.0000001)
2137
0
        nUOMLengthCode = 9003;  // US survey foot.
2138
0
    else if (fabs(dfLinearUOM - 1.0) > 0.00000001)
2139
0
        nUOMLengthCode = KvUserDefined;
2140
2141
    /* -------------------------------------------------------------------- */
2142
    /*      Get some authority values.                                      */
2143
    /* -------------------------------------------------------------------- */
2144
0
    int nPCS = KvUserDefined;
2145
2146
0
    if (poSRS->GetAuthorityName("PROJCS") != nullptr &&
2147
0
        EQUAL(poSRS->GetAuthorityName("PROJCS"), "EPSG"))
2148
0
    {
2149
0
        nPCS = atoi(poSRS->GetAuthorityCode("PROJCS"));
2150
0
        if (nPCS > 32767)
2151
0
            nPCS = KvUserDefined;
2152
0
    }
2153
2154
    /* -------------------------------------------------------------------- */
2155
    /*      Handle the projection transformation.                           */
2156
    /* -------------------------------------------------------------------- */
2157
0
    bool bWritePEString = false;
2158
0
    bool bUnknownProjection = false;
2159
2160
0
    const OGRSpatialReference *poSRSCompatibleOfWKT1 = poSRS;
2161
0
#if PROJ_AT_LEAST_VERSION(6, 3, 0)
2162
0
    OGRSpatialReference oSRSCompatibleOfWKT1;
2163
0
    if (poSRS->IsProjected() && !poSRS->IsCompound() &&
2164
0
        poSRS->GetAxesCount() == 3)
2165
0
    {
2166
0
        oSRSCompatibleOfWKT1 = *poSRS;
2167
0
        oSRSCompatibleOfWKT1.DemoteTo2D(nullptr);
2168
0
        poSRSCompatibleOfWKT1 = &oSRSCompatibleOfWKT1;
2169
0
        bWritePEString = true;
2170
0
    }
2171
0
#endif
2172
2173
0
    double dfAngUnitValue = 0;
2174
0
    std::string osAngUnitName;
2175
0
    if (bHasEllipsoid)
2176
0
    {
2177
0
        const char *angUnitNameTmp = "";
2178
0
        dfAngUnitValue = poSRS->GetAngularUnits(&angUnitNameTmp);
2179
0
        osAngUnitName = angUnitNameTmp;
2180
0
    }
2181
2182
    // Convert angular projection parameters from its normalized value in degree
2183
    // to the units of GeogAngularUnitsGeoKey.
2184
    // Note: for GDAL <= 3.9.0, we always have written them in degrees !
2185
    // We can set GTIFF_WRITE_ANGULAR_PARAMS_IN_DEGREE=YES to get that
2186
    // non-conformant behavior...
2187
0
    const auto ConvertAngularParam = [dfAngUnitValue](double dfValInDeg)
2188
0
    {
2189
0
        constexpr double DEG_TO_RAD = M_PI / 180.0;
2190
0
        return dfAngUnitValue != 0 &&
2191
0
                       std::fabs(dfAngUnitValue - DEG_TO_RAD) > 1e-10 &&
2192
0
                       !CPLTestBool(CPLGetConfigOption(
2193
0
                           "GTIFF_WRITE_ANGULAR_PARAMS_IN_DEGREE", "NO"))
2194
0
                   ? dfValInDeg * DEG_TO_RAD / dfAngUnitValue
2195
0
                   : dfValInDeg;
2196
0
    };
2197
2198
0
    const char *pszProjection =
2199
0
        poSRSCompatibleOfWKT1->GetAttrValue("PROJECTION");
2200
0
    if (nPCS != KvUserDefined)
2201
0
    {
2202
        // If ESRI_PE flavor is explicitly required, then for EPSG:3857
2203
        // we will have to write a completely non-standard definition
2204
        // that requires not setting GTModelTypeGeoKey to ProjectedCSTypeGeoKey.
2205
0
        if (eFlavor == GEOTIFF_KEYS_ESRI_PE && nPCS == 3857)
2206
0
        {
2207
0
            bWritePEString = true;
2208
0
        }
2209
0
        else
2210
0
        {
2211
0
            GTIFKeySet(psGTIF, GTModelTypeGeoKey, TYPE_SHORT, 1,
2212
0
                       ModelTypeProjected);
2213
0
            GTIFKeySet(psGTIF, ProjectedCSTypeGeoKey, TYPE_SHORT, 1, nPCS);
2214
0
        }
2215
0
    }
2216
0
    else if (poSRSCompatibleOfWKT1->IsGeocentric())
2217
0
    {
2218
0
        GTIFKeySet(psGTIF, GTModelTypeGeoKey, TYPE_SHORT, 1,
2219
0
                   ModelTypeGeocentric);
2220
0
    }
2221
0
    else if (pszProjection == nullptr)
2222
0
    {
2223
0
        if (poSRSCompatibleOfWKT1->IsGeographic())
2224
0
            GTIFKeySet(psGTIF, GTModelTypeGeoKey, TYPE_SHORT, 1,
2225
0
                       ModelTypeGeographic);
2226
        // Otherwise, presumably something like LOCAL_CS.
2227
0
    }
2228
0
    else if (EQUAL(pszProjection, SRS_PT_ALBERS_CONIC_EQUAL_AREA))
2229
0
    {
2230
0
        GTIFKeySet(psGTIF, GTModelTypeGeoKey, TYPE_SHORT, 1,
2231
0
                   ModelTypeProjected);
2232
0
        GTIFKeySet(psGTIF, ProjectedCSTypeGeoKey, TYPE_SHORT, 1, KvUserDefined);
2233
0
        GTIFKeySet(psGTIF, ProjectionGeoKey, TYPE_SHORT, 1, KvUserDefined);
2234
2235
0
        GTIFKeySet(psGTIF, ProjCoordTransGeoKey, TYPE_SHORT, 1,
2236
0
                   CT_AlbersEqualArea);
2237
2238
0
        GTIFKeySet(psGTIF, ProjStdParallelGeoKey, TYPE_DOUBLE, 1,
2239
0
                   ConvertAngularParam(poSRSCompatibleOfWKT1->GetNormProjParm(
2240
0
                       SRS_PP_STANDARD_PARALLEL_1, 0.0)));
2241
2242
0
        GTIFKeySet(psGTIF, ProjStdParallel2GeoKey, TYPE_DOUBLE, 1,
2243
0
                   ConvertAngularParam(poSRSCompatibleOfWKT1->GetNormProjParm(
2244
0
                       SRS_PP_STANDARD_PARALLEL_2, 0.0)));
2245
2246
0
        GTIFKeySet(psGTIF, ProjNatOriginLatGeoKey, TYPE_DOUBLE, 1,
2247
0
                   ConvertAngularParam(poSRSCompatibleOfWKT1->GetNormProjParm(
2248
0
                       SRS_PP_LATITUDE_OF_CENTER, 0.0)));
2249
2250
0
        GTIFKeySet(psGTIF, ProjNatOriginLongGeoKey, TYPE_DOUBLE, 1,
2251
0
                   ConvertAngularParam(poSRSCompatibleOfWKT1->GetNormProjParm(
2252
0
                       SRS_PP_LONGITUDE_OF_CENTER, 0.0)));
2253
2254
0
        GTIFKeySet(
2255
0
            psGTIF, ProjFalseEastingGeoKey, TYPE_DOUBLE, 1,
2256
0
            poSRSCompatibleOfWKT1->GetProjParm(SRS_PP_FALSE_EASTING, 0.0));
2257
2258
0
        GTIFKeySet(
2259
0
            psGTIF, ProjFalseNorthingGeoKey, TYPE_DOUBLE, 1,
2260
0
            poSRSCompatibleOfWKT1->GetProjParm(SRS_PP_FALSE_NORTHING, 0.0));
2261
0
    }
2262
2263
0
    else if (poSRSCompatibleOfWKT1->GetUTMZone() != 0)
2264
0
    {
2265
0
        GTIFKeySet(psGTIF, GTModelTypeGeoKey, TYPE_SHORT, 1,
2266
0
                   ModelTypeProjected);
2267
2268
0
        int bNorth = 0;
2269
0
        const int nZone = poSRSCompatibleOfWKT1->GetUTMZone(&bNorth);
2270
2271
0
        if (nDatum == Datum_North_American_Datum_1983 && nZone >= 3 &&
2272
0
            nZone <= 22 && bNorth && nUOMLengthCode == 9001)
2273
0
        {
2274
0
            nPCS = 26900 + nZone;
2275
2276
0
            GTIFKeySet(psGTIF, ProjectedCSTypeGeoKey, TYPE_SHORT, 1, nPCS);
2277
0
        }
2278
0
        else if (nDatum == Datum_North_American_Datum_1927 && nZone >= 3 &&
2279
0
                 nZone <= 22 && bNorth && nUOMLengthCode == 9001)
2280
0
        {
2281
0
            nPCS = 26700 + nZone;
2282
2283
0
            GTIFKeySet(psGTIF, ProjectedCSTypeGeoKey, TYPE_SHORT, 1, nPCS);
2284
0
        }
2285
0
        else if (nDatum == Datum_WGS84 && nUOMLengthCode == 9001)
2286
0
        {
2287
0
            if (bNorth)
2288
0
                nPCS = 32600 + nZone;
2289
0
            else
2290
0
                nPCS = 32700 + nZone;
2291
2292
0
            GTIFKeySet(psGTIF, ProjectedCSTypeGeoKey, TYPE_SHORT, 1, nPCS);
2293
0
        }
2294
0
        else
2295
0
        {
2296
0
            const int nProjection = nZone + (bNorth ? 16000 : 16100);
2297
0
            GTIFKeySet(psGTIF, ProjectedCSTypeGeoKey, TYPE_SHORT, 1,
2298
0
                       KvUserDefined);
2299
2300
0
            GTIFKeySet(psGTIF, ProjectionGeoKey, TYPE_SHORT, 1, nProjection);
2301
0
        }
2302
0
    }
2303
2304
0
    else if (EQUAL(pszProjection, SRS_PT_TRANSVERSE_MERCATOR))
2305
0
    {
2306
0
        GTIFKeySet(psGTIF, GTModelTypeGeoKey, TYPE_SHORT, 1,
2307
0
                   ModelTypeProjected);
2308
0
        GTIFKeySet(psGTIF, ProjectedCSTypeGeoKey, TYPE_SHORT, 1, KvUserDefined);
2309
0
        GTIFKeySet(psGTIF, ProjectionGeoKey, TYPE_SHORT, 1, KvUserDefined);
2310
2311
0
        GTIFKeySet(psGTIF, ProjCoordTransGeoKey, TYPE_SHORT, 1,
2312
0
                   CT_TransverseMercator);
2313
2314
0
        GTIFKeySet(psGTIF, ProjNatOriginLatGeoKey, TYPE_DOUBLE, 1,
2315
0
                   ConvertAngularParam(poSRSCompatibleOfWKT1->GetNormProjParm(
2316
0
                       SRS_PP_LATITUDE_OF_ORIGIN, 0.0)));
2317
2318
0
        GTIFKeySet(psGTIF, ProjNatOriginLongGeoKey, TYPE_DOUBLE, 1,
2319
0
                   ConvertAngularParam(poSRSCompatibleOfWKT1->GetNormProjParm(
2320
0
                       SRS_PP_CENTRAL_MERIDIAN, 0.0)));
2321
2322
0
        GTIFKeySet(
2323
0
            psGTIF, ProjScaleAtNatOriginGeoKey, TYPE_DOUBLE, 1,
2324
0
            poSRSCompatibleOfWKT1->GetNormProjParm(SRS_PP_SCALE_FACTOR, 1.0));
2325
2326
0
        GTIFKeySet(
2327
0
            psGTIF, ProjFalseEastingGeoKey, TYPE_DOUBLE, 1,
2328
0
            poSRSCompatibleOfWKT1->GetProjParm(SRS_PP_FALSE_EASTING, 0.0));
2329
2330
0
        GTIFKeySet(
2331
0
            psGTIF, ProjFalseNorthingGeoKey, TYPE_DOUBLE, 1,
2332
0
            poSRSCompatibleOfWKT1->GetProjParm(SRS_PP_FALSE_NORTHING, 0.0));
2333
0
    }
2334
2335
0
    else if (EQUAL(pszProjection, SRS_PT_TRANSVERSE_MERCATOR_SOUTH_ORIENTED))
2336
0
    {
2337
0
        GTIFKeySet(psGTIF, GTModelTypeGeoKey, TYPE_SHORT, 1,
2338
0
                   ModelTypeProjected);
2339
0
        GTIFKeySet(psGTIF, ProjectedCSTypeGeoKey, TYPE_SHORT, 1, KvUserDefined);
2340
0
        GTIFKeySet(psGTIF, ProjectionGeoKey, TYPE_SHORT, 1, KvUserDefined);
2341
2342
0
        GTIFKeySet(psGTIF, ProjCoordTransGeoKey, TYPE_SHORT, 1,
2343
0
                   CT_TransvMercator_SouthOriented);
2344
2345
0
        GTIFKeySet(psGTIF, ProjNatOriginLatGeoKey, TYPE_DOUBLE, 1,
2346
0
                   ConvertAngularParam(poSRSCompatibleOfWKT1->GetNormProjParm(
2347
0
                       SRS_PP_LATITUDE_OF_ORIGIN, 0.0)));
2348
2349
0
        GTIFKeySet(psGTIF, ProjNatOriginLongGeoKey, TYPE_DOUBLE, 1,
2350
0
                   ConvertAngularParam(poSRSCompatibleOfWKT1->GetNormProjParm(
2351
0
                       SRS_PP_CENTRAL_MERIDIAN, 0.0)));
2352
2353
0
        GTIFKeySet(
2354
0
            psGTIF, ProjScaleAtNatOriginGeoKey, TYPE_DOUBLE, 1,
2355
0
            poSRSCompatibleOfWKT1->GetNormProjParm(SRS_PP_SCALE_FACTOR, 1.0));
2356
2357
0
        GTIFKeySet(
2358
0
            psGTIF, ProjFalseEastingGeoKey, TYPE_DOUBLE, 1,
2359
0
            poSRSCompatibleOfWKT1->GetProjParm(SRS_PP_FALSE_EASTING, 0.0));
2360
2361
0
        GTIFKeySet(
2362
0
            psGTIF, ProjFalseNorthingGeoKey, TYPE_DOUBLE, 1,
2363
0
            poSRSCompatibleOfWKT1->GetProjParm(SRS_PP_FALSE_NORTHING, 0.0));
2364
0
    }
2365
2366
0
    else if (EQUAL(pszProjection, SRS_PT_MERCATOR_2SP) ||
2367
0
             EQUAL(pszProjection, SRS_PT_MERCATOR_1SP))
2368
2369
0
    {
2370
0
        GTIFKeySet(psGTIF, GTModelTypeGeoKey, TYPE_SHORT, 1,
2371
0
                   ModelTypeProjected);
2372
0
        GTIFKeySet(psGTIF, ProjectedCSTypeGeoKey, TYPE_SHORT, 1, KvUserDefined);
2373
0
        GTIFKeySet(psGTIF, ProjectionGeoKey, TYPE_SHORT, 1, KvUserDefined);
2374
2375
0
        GTIFKeySet(psGTIF, ProjCoordTransGeoKey, TYPE_SHORT, 1, CT_Mercator);
2376
2377
0
        GTIFKeySet(psGTIF, ProjNatOriginLatGeoKey, TYPE_DOUBLE, 1,
2378
0
                   ConvertAngularParam(poSRSCompatibleOfWKT1->GetNormProjParm(
2379
0
                       SRS_PP_LATITUDE_OF_ORIGIN, 0.0)));
2380
2381
0
        GTIFKeySet(psGTIF, ProjNatOriginLongGeoKey, TYPE_DOUBLE, 1,
2382
0
                   ConvertAngularParam(poSRSCompatibleOfWKT1->GetNormProjParm(
2383
0
                       SRS_PP_CENTRAL_MERIDIAN, 0.0)));
2384
2385
0
        if (EQUAL(pszProjection, SRS_PT_MERCATOR_2SP))
2386
0
            GTIFKeySet(
2387
0
                psGTIF, ProjStdParallel1GeoKey, TYPE_DOUBLE, 1,
2388
0
                ConvertAngularParam(poSRSCompatibleOfWKT1->GetNormProjParm(
2389
0
                    SRS_PP_STANDARD_PARALLEL_1, 0.0)));
2390
0
        else
2391
0
            GTIFKeySet(psGTIF, ProjScaleAtNatOriginGeoKey, TYPE_DOUBLE, 1,
2392
0
                       poSRSCompatibleOfWKT1->GetNormProjParm(
2393
0
                           SRS_PP_SCALE_FACTOR, 1.0));
2394
2395
0
        GTIFKeySet(
2396
0
            psGTIF, ProjFalseEastingGeoKey, TYPE_DOUBLE, 1,
2397
0
            poSRSCompatibleOfWKT1->GetProjParm(SRS_PP_FALSE_EASTING, 0.0));
2398
2399
0
        GTIFKeySet(
2400
0
            psGTIF, ProjFalseNorthingGeoKey, TYPE_DOUBLE, 1,
2401
0
            poSRSCompatibleOfWKT1->GetProjParm(SRS_PP_FALSE_NORTHING, 0.0));
2402
0
    }
2403
2404
0
    else if (EQUAL(pszProjection, SRS_PT_OBLIQUE_STEREOGRAPHIC))
2405
0
    {
2406
0
        GTIFKeySet(psGTIF, GTModelTypeGeoKey, TYPE_SHORT, 1,
2407
0
                   ModelTypeProjected);
2408
0
        GTIFKeySet(psGTIF, ProjectedCSTypeGeoKey, TYPE_SHORT, 1, KvUserDefined);
2409
0
        GTIFKeySet(psGTIF, ProjectionGeoKey, TYPE_SHORT, 1, KvUserDefined);
2410
2411
0
        GTIFKeySet(psGTIF, ProjCoordTransGeoKey, TYPE_SHORT, 1,
2412
0
                   CT_ObliqueStereographic);
2413
2414
0
        GTIFKeySet(psGTIF, ProjNatOriginLatGeoKey, TYPE_DOUBLE, 1,
2415
0
                   ConvertAngularParam(poSRSCompatibleOfWKT1->GetNormProjParm(
2416
0
                       SRS_PP_LATITUDE_OF_ORIGIN, 0.0)));
2417
2418
0
        GTIFKeySet(psGTIF, ProjNatOriginLongGeoKey, TYPE_DOUBLE, 1,
2419
0
                   ConvertAngularParam(poSRSCompatibleOfWKT1->GetNormProjParm(
2420
0
                       SRS_PP_CENTRAL_MERIDIAN, 0.0)));
2421
2422
0
        GTIFKeySet(
2423
0
            psGTIF, ProjScaleAtNatOriginGeoKey, TYPE_DOUBLE, 1,
2424
0
            poSRSCompatibleOfWKT1->GetNormProjParm(SRS_PP_SCALE_FACTOR, 1.0));
2425
2426
0
        GTIFKeySet(
2427
0
            psGTIF, ProjFalseEastingGeoKey, TYPE_DOUBLE, 1,
2428
0
            poSRSCompatibleOfWKT1->GetProjParm(SRS_PP_FALSE_EASTING, 0.0));
2429
2430
0
        GTIFKeySet(
2431
0
            psGTIF, ProjFalseNorthingGeoKey, TYPE_DOUBLE, 1,
2432
0
            poSRSCompatibleOfWKT1->GetProjParm(SRS_PP_FALSE_NORTHING, 0.0));
2433
0
    }
2434
2435
0
    else if (EQUAL(pszProjection, SRS_PT_STEREOGRAPHIC))
2436
0
    {
2437
0
        GTIFKeySet(psGTIF, GTModelTypeGeoKey, TYPE_SHORT, 1,
2438
0
                   ModelTypeProjected);
2439
0
        GTIFKeySet(psGTIF, ProjectedCSTypeGeoKey, TYPE_SHORT, 1, KvUserDefined);
2440
0
        GTIFKeySet(psGTIF, ProjectionGeoKey, TYPE_SHORT, 1, KvUserDefined);
2441
2442
0
        GTIFKeySet(psGTIF, ProjCoordTransGeoKey, TYPE_SHORT, 1,
2443
0
                   CT_Stereographic);
2444
2445
0
        GTIFKeySet(psGTIF, ProjCenterLatGeoKey, TYPE_DOUBLE, 1,
2446
0
                   ConvertAngularParam(poSRSCompatibleOfWKT1->GetNormProjParm(
2447
0
                       SRS_PP_LATITUDE_OF_ORIGIN, 0.0)));
2448
2449
0
        GTIFKeySet(psGTIF, ProjCenterLongGeoKey, TYPE_DOUBLE, 1,
2450
0
                   ConvertAngularParam(poSRSCompatibleOfWKT1->GetNormProjParm(
2451
0
                       SRS_PP_CENTRAL_MERIDIAN, 0.0)));
2452
2453
0
        GTIFKeySet(
2454
0
            psGTIF, ProjScaleAtNatOriginGeoKey, TYPE_DOUBLE, 1,
2455
0
            poSRSCompatibleOfWKT1->GetNormProjParm(SRS_PP_SCALE_FACTOR, 1.0));
2456
2457
0
        GTIFKeySet(
2458
0
            psGTIF, ProjFalseEastingGeoKey, TYPE_DOUBLE, 1,
2459
0
            poSRSCompatibleOfWKT1->GetProjParm(SRS_PP_FALSE_EASTING, 0.0));
2460
2461
0
        GTIFKeySet(
2462
0
            psGTIF, ProjFalseNorthingGeoKey, TYPE_DOUBLE, 1,
2463
0
            poSRSCompatibleOfWKT1->GetProjParm(SRS_PP_FALSE_NORTHING, 0.0));
2464
0
    }
2465
2466
0
    else if (EQUAL(pszProjection, SRS_PT_POLAR_STEREOGRAPHIC))
2467
0
    {
2468
0
        GTIFKeySet(psGTIF, GTModelTypeGeoKey, TYPE_SHORT, 1,
2469
0
                   ModelTypeProjected);
2470
0
        GTIFKeySet(psGTIF, ProjectedCSTypeGeoKey, TYPE_SHORT, 1, KvUserDefined);
2471
0
        GTIFKeySet(psGTIF, ProjectionGeoKey, TYPE_SHORT, 1, KvUserDefined);
2472
2473
0
        GTIFKeySet(psGTIF, ProjCoordTransGeoKey, TYPE_SHORT, 1,
2474
0
                   CT_PolarStereographic);
2475
2476
0
        GTIFKeySet(psGTIF, ProjNatOriginLatGeoKey, TYPE_DOUBLE, 1,
2477
0
                   ConvertAngularParam(poSRSCompatibleOfWKT1->GetNormProjParm(
2478
0
                       SRS_PP_LATITUDE_OF_ORIGIN, 0.0)));
2479
2480
0
        GTIFKeySet(psGTIF, ProjStraightVertPoleLongGeoKey, TYPE_DOUBLE, 1,
2481
0
                   ConvertAngularParam(poSRSCompatibleOfWKT1->GetNormProjParm(
2482
0
                       SRS_PP_CENTRAL_MERIDIAN, 0.0)));
2483
2484
0
        GTIFKeySet(
2485
0
            psGTIF, ProjScaleAtNatOriginGeoKey, TYPE_DOUBLE, 1,
2486
0
            poSRSCompatibleOfWKT1->GetNormProjParm(SRS_PP_SCALE_FACTOR, 1.0));
2487
2488
0
        GTIFKeySet(
2489
0
            psGTIF, ProjFalseEastingGeoKey, TYPE_DOUBLE, 1,
2490
0
            poSRSCompatibleOfWKT1->GetProjParm(SRS_PP_FALSE_EASTING, 0.0));
2491
2492
0
        GTIFKeySet(
2493
0
            psGTIF, ProjFalseNorthingGeoKey, TYPE_DOUBLE, 1,
2494
0
            poSRSCompatibleOfWKT1->GetProjParm(SRS_PP_FALSE_NORTHING, 0.0));
2495
0
    }
2496
2497
0
    else if (EQUAL(pszProjection, SRS_PT_HOTINE_OBLIQUE_MERCATOR))
2498
0
    {
2499
0
        GTIFKeySet(psGTIF, GTModelTypeGeoKey, TYPE_SHORT, 1,
2500
0
                   ModelTypeProjected);
2501
0
        GTIFKeySet(psGTIF, ProjectedCSTypeGeoKey, TYPE_SHORT, 1, KvUserDefined);
2502
0
        GTIFKeySet(psGTIF, ProjectionGeoKey, TYPE_SHORT, 1, KvUserDefined);
2503
2504
0
        GTIFKeySet(psGTIF, ProjCoordTransGeoKey, TYPE_SHORT, 1,
2505
0
                   CT_ObliqueMercator);
2506
2507
0
        GTIFKeySet(psGTIF, ProjCenterLatGeoKey, TYPE_DOUBLE, 1,
2508
0
                   ConvertAngularParam(poSRSCompatibleOfWKT1->GetNormProjParm(
2509
0
                       SRS_PP_LATITUDE_OF_CENTER, 0.0)));
2510
2511
0
        GTIFKeySet(psGTIF, ProjCenterLongGeoKey, TYPE_DOUBLE, 1,
2512
0
                   ConvertAngularParam(poSRSCompatibleOfWKT1->GetNormProjParm(
2513
0
                       SRS_PP_LONGITUDE_OF_CENTER, 0.0)));
2514
2515
0
        GTIFKeySet(psGTIF, ProjAzimuthAngleGeoKey, TYPE_DOUBLE, 1,
2516
0
                   poSRSCompatibleOfWKT1->GetNormProjParm(SRS_PP_AZIMUTH, 0.0));
2517
2518
0
        GTIFKeySet(psGTIF, ProjRectifiedGridAngleGeoKey, TYPE_DOUBLE, 1,
2519
0
                   ConvertAngularParam(poSRSCompatibleOfWKT1->GetNormProjParm(
2520
0
                       SRS_PP_RECTIFIED_GRID_ANGLE, 0.0)));
2521
2522
0
        GTIFKeySet(
2523
0
            psGTIF, ProjScaleAtCenterGeoKey, TYPE_DOUBLE, 1,
2524
0
            poSRSCompatibleOfWKT1->GetNormProjParm(SRS_PP_SCALE_FACTOR, 1.0));
2525
2526
0
        GTIFKeySet(
2527
0
            psGTIF, ProjFalseEastingGeoKey, TYPE_DOUBLE, 1,
2528
0
            poSRSCompatibleOfWKT1->GetProjParm(SRS_PP_FALSE_EASTING, 0.0));
2529
2530
0
        GTIFKeySet(
2531
0
            psGTIF, ProjFalseNorthingGeoKey, TYPE_DOUBLE, 1,
2532
0
            poSRSCompatibleOfWKT1->GetProjParm(SRS_PP_FALSE_NORTHING, 0.0));
2533
0
    }
2534
2535
0
    else if (EQUAL(pszProjection,
2536
0
                   SRS_PT_HOTINE_OBLIQUE_MERCATOR_AZIMUTH_CENTER))
2537
0
    {
2538
0
        GTIFKeySet(psGTIF, GTModelTypeGeoKey, TYPE_SHORT, 1,
2539
0
                   ModelTypeProjected);
2540
0
        GTIFKeySet(psGTIF, ProjectedCSTypeGeoKey, TYPE_SHORT, 1, KvUserDefined);
2541
0
        GTIFKeySet(psGTIF, ProjectionGeoKey, TYPE_SHORT, 1, KvUserDefined);
2542
2543
0
        GTIFKeySet(psGTIF, ProjCoordTransGeoKey, TYPE_SHORT, 1,
2544
0
                   CT_HotineObliqueMercatorAzimuthCenter);
2545
2546
0
        GTIFKeySet(psGTIF, ProjCenterLatGeoKey, TYPE_DOUBLE, 1,
2547
0
                   ConvertAngularParam(poSRSCompatibleOfWKT1->GetNormProjParm(
2548
0
                       SRS_PP_LATITUDE_OF_CENTER, 0.0)));
2549
2550
0
        GTIFKeySet(psGTIF, ProjCenterLongGeoKey, TYPE_DOUBLE, 1,
2551
0
                   ConvertAngularParam(poSRSCompatibleOfWKT1->GetNormProjParm(
2552
0
                       SRS_PP_LONGITUDE_OF_CENTER, 0.0)));
2553
2554
0
        GTIFKeySet(psGTIF, ProjAzimuthAngleGeoKey, TYPE_DOUBLE, 1,
2555
0
                   poSRSCompatibleOfWKT1->GetNormProjParm(SRS_PP_AZIMUTH, 0.0));
2556
2557
0
        GTIFKeySet(psGTIF, ProjRectifiedGridAngleGeoKey, TYPE_DOUBLE, 1,
2558
0
                   ConvertAngularParam(poSRSCompatibleOfWKT1->GetNormProjParm(
2559
0
                       SRS_PP_RECTIFIED_GRID_ANGLE, 0.0)));
2560
2561
0
        GTIFKeySet(
2562
0
            psGTIF, ProjScaleAtCenterGeoKey, TYPE_DOUBLE, 1,
2563
0
            poSRSCompatibleOfWKT1->GetNormProjParm(SRS_PP_SCALE_FACTOR, 1.0));
2564
2565
0
        GTIFKeySet(
2566
0
            psGTIF, ProjFalseEastingGeoKey, TYPE_DOUBLE, 1,
2567
0
            poSRSCompatibleOfWKT1->GetProjParm(SRS_PP_FALSE_EASTING, 0.0));
2568
2569
0
        GTIFKeySet(
2570
0
            psGTIF, ProjFalseNorthingGeoKey, TYPE_DOUBLE, 1,
2571
0
            poSRSCompatibleOfWKT1->GetProjParm(SRS_PP_FALSE_NORTHING, 0.0));
2572
0
    }
2573
2574
0
    else if (EQUAL(pszProjection, "Laborde_Oblique_Mercator"))
2575
0
    {
2576
0
        GTIFKeySet(psGTIF, GTModelTypeGeoKey, TYPE_SHORT, 1,
2577
0
                   ModelTypeProjected);
2578
0
        GTIFKeySet(psGTIF, ProjectedCSTypeGeoKey, TYPE_SHORT, 1, KvUserDefined);
2579
0
        GTIFKeySet(psGTIF, ProjectionGeoKey, TYPE_SHORT, 1, KvUserDefined);
2580
2581
0
        GTIFKeySet(psGTIF, ProjCoordTransGeoKey, TYPE_SHORT, 1,
2582
0
                   CT_ObliqueMercator_Laborde);
2583
2584
0
        GTIFKeySet(psGTIF, ProjCenterLatGeoKey, TYPE_DOUBLE, 1,
2585
0
                   ConvertAngularParam(poSRSCompatibleOfWKT1->GetNormProjParm(
2586
0
                       SRS_PP_LATITUDE_OF_CENTER, 0.0)));
2587
2588
0
        GTIFKeySet(psGTIF, ProjCenterLongGeoKey, TYPE_DOUBLE, 1,
2589
0
                   ConvertAngularParam(poSRSCompatibleOfWKT1->GetNormProjParm(
2590
0
                       SRS_PP_LONGITUDE_OF_CENTER, 0.0)));
2591
2592
0
        GTIFKeySet(psGTIF, ProjAzimuthAngleGeoKey, TYPE_DOUBLE, 1,
2593
0
                   poSRSCompatibleOfWKT1->GetNormProjParm(SRS_PP_AZIMUTH, 0.0));
2594
2595
0
        GTIFKeySet(
2596
0
            psGTIF, ProjScaleAtCenterGeoKey, TYPE_DOUBLE, 1,
2597
0
            poSRSCompatibleOfWKT1->GetNormProjParm(SRS_PP_SCALE_FACTOR, 1.0));
2598
2599
0
        GTIFKeySet(
2600
0
            psGTIF, ProjFalseEastingGeoKey, TYPE_DOUBLE, 1,
2601
0
            poSRSCompatibleOfWKT1->GetProjParm(SRS_PP_FALSE_EASTING, 0.0));
2602
2603
0
        GTIFKeySet(
2604
0
            psGTIF, ProjFalseNorthingGeoKey, TYPE_DOUBLE, 1,
2605
0
            poSRSCompatibleOfWKT1->GetProjParm(SRS_PP_FALSE_NORTHING, 0.0));
2606
0
    }
2607
2608
0
    else if (EQUAL(pszProjection, SRS_PT_CASSINI_SOLDNER))
2609
0
    {
2610
0
        GTIFKeySet(psGTIF, GTModelTypeGeoKey, TYPE_SHORT, 1,
2611
0
                   ModelTypeProjected);
2612
0
        GTIFKeySet(psGTIF, ProjectedCSTypeGeoKey, TYPE_SHORT, 1, KvUserDefined);
2613
0
        GTIFKeySet(psGTIF, ProjectionGeoKey, TYPE_SHORT, 1, KvUserDefined);
2614
2615
0
        GTIFKeySet(psGTIF, ProjCoordTransGeoKey, TYPE_SHORT, 1,
2616
0
                   CT_CassiniSoldner);
2617
2618
0
        GTIFKeySet(psGTIF, ProjNatOriginLatGeoKey, TYPE_DOUBLE, 1,
2619
0
                   ConvertAngularParam(poSRSCompatibleOfWKT1->GetNormProjParm(
2620
0
                       SRS_PP_LATITUDE_OF_ORIGIN, 0.0)));
2621
2622
0
        GTIFKeySet(psGTIF, ProjNatOriginLongGeoKey, TYPE_DOUBLE, 1,
2623
0
                   ConvertAngularParam(poSRSCompatibleOfWKT1->GetNormProjParm(
2624
0
                       SRS_PP_CENTRAL_MERIDIAN, 0.0)));
2625
2626
0
        GTIFKeySet(
2627
0
            psGTIF, ProjFalseEastingGeoKey, TYPE_DOUBLE, 1,
2628
0
            poSRSCompatibleOfWKT1->GetProjParm(SRS_PP_FALSE_EASTING, 0.0));
2629
2630
0
        GTIFKeySet(
2631
0
            psGTIF, ProjFalseNorthingGeoKey, TYPE_DOUBLE, 1,
2632
0
            poSRSCompatibleOfWKT1->GetProjParm(SRS_PP_FALSE_NORTHING, 0.0));
2633
0
    }
2634
2635
0
    else if (EQUAL(pszProjection, SRS_PT_EQUIDISTANT_CONIC))
2636
0
    {
2637
0
        GTIFKeySet(psGTIF, GTModelTypeGeoKey, TYPE_SHORT, 1,
2638
0
                   ModelTypeProjected);
2639
0
        GTIFKeySet(psGTIF, ProjectedCSTypeGeoKey, TYPE_SHORT, 1, KvUserDefined);
2640
0
        GTIFKeySet(psGTIF, ProjectionGeoKey, TYPE_SHORT, 1, KvUserDefined);
2641
2642
0
        GTIFKeySet(psGTIF, ProjCoordTransGeoKey, TYPE_SHORT, 1,
2643
0
                   CT_EquidistantConic);
2644
2645
0
        GTIFKeySet(psGTIF, ProjStdParallel1GeoKey, TYPE_DOUBLE, 1,
2646
0
                   ConvertAngularParam(poSRSCompatibleOfWKT1->GetNormProjParm(
2647
0
                       SRS_PP_STANDARD_PARALLEL_1, 0.0)));
2648
2649
0
        GTIFKeySet(psGTIF, ProjStdParallel2GeoKey, TYPE_DOUBLE, 1,
2650
0
                   ConvertAngularParam(poSRSCompatibleOfWKT1->GetNormProjParm(
2651
0
                       SRS_PP_STANDARD_PARALLEL_2, 0.0)));
2652
2653
0
        GTIFKeySet(psGTIF, ProjNatOriginLatGeoKey, TYPE_DOUBLE, 1,
2654
0
                   ConvertAngularParam(poSRSCompatibleOfWKT1->GetNormProjParm(
2655
0
                       SRS_PP_LATITUDE_OF_CENTER, 0.0)));
2656
2657
0
        GTIFKeySet(psGTIF, ProjNatOriginLongGeoKey, TYPE_DOUBLE, 1,
2658
0
                   ConvertAngularParam(poSRSCompatibleOfWKT1->GetNormProjParm(
2659
0
                       SRS_PP_LONGITUDE_OF_CENTER, 0.0)));
2660
2661
0
        GTIFKeySet(
2662
0
            psGTIF, ProjFalseEastingGeoKey, TYPE_DOUBLE, 1,
2663
0
            poSRSCompatibleOfWKT1->GetProjParm(SRS_PP_FALSE_EASTING, 0.0));
2664
2665
0
        GTIFKeySet(
2666
0
            psGTIF, ProjFalseNorthingGeoKey, TYPE_DOUBLE, 1,
2667
0
            poSRSCompatibleOfWKT1->GetProjParm(SRS_PP_FALSE_NORTHING, 0.0));
2668
0
    }
2669
2670
0
    else if (EQUAL(pszProjection, SRS_PT_POLYCONIC))
2671
0
    {
2672
0
        GTIFKeySet(psGTIF, GTModelTypeGeoKey, TYPE_SHORT, 1,
2673
0
                   ModelTypeProjected);
2674
0
        GTIFKeySet(psGTIF, ProjectedCSTypeGeoKey, TYPE_SHORT, 1, KvUserDefined);
2675
0
        GTIFKeySet(psGTIF, ProjectionGeoKey, TYPE_SHORT, 1, KvUserDefined);
2676
2677
0
        GTIFKeySet(psGTIF, ProjCoordTransGeoKey, TYPE_SHORT, 1, CT_Polyconic);
2678
2679
0
        GTIFKeySet(psGTIF, ProjNatOriginLatGeoKey, TYPE_DOUBLE, 1,
2680
0
                   poSRSCompatibleOfWKT1->GetNormProjParm(
2681
0
                       SRS_PP_LATITUDE_OF_ORIGIN, 0.0));
2682
2683
0
        GTIFKeySet(psGTIF, ProjNatOriginLongGeoKey, TYPE_DOUBLE, 1,
2684
0
                   poSRSCompatibleOfWKT1->GetNormProjParm(
2685
0
                       SRS_PP_CENTRAL_MERIDIAN, 0.0));
2686
2687
0
        GTIFKeySet(
2688
0
            psGTIF, ProjScaleAtNatOriginGeoKey, TYPE_DOUBLE, 1,
2689
0
            poSRSCompatibleOfWKT1->GetNormProjParm(SRS_PP_SCALE_FACTOR, 1.0));
2690
2691
0
        GTIFKeySet(
2692
0
            psGTIF, ProjFalseEastingGeoKey, TYPE_DOUBLE, 1,
2693
0
            poSRSCompatibleOfWKT1->GetProjParm(SRS_PP_FALSE_EASTING, 0.0));
2694
2695
0
        GTIFKeySet(
2696
0
            psGTIF, ProjFalseNorthingGeoKey, TYPE_DOUBLE, 1,
2697
0
            poSRSCompatibleOfWKT1->GetProjParm(SRS_PP_FALSE_NORTHING, 0.0));
2698
0
    }
2699
2700
0
    else if (EQUAL(pszProjection, SRS_PT_AZIMUTHAL_EQUIDISTANT))
2701
0
    {
2702
0
        GTIFKeySet(psGTIF, GTModelTypeGeoKey, TYPE_SHORT, 1,
2703
0
                   ModelTypeProjected);
2704
0
        GTIFKeySet(psGTIF, ProjectedCSTypeGeoKey, TYPE_SHORT, 1, KvUserDefined);
2705
0
        GTIFKeySet(psGTIF, ProjectionGeoKey, TYPE_SHORT, 1, KvUserDefined);
2706
2707
0
        GTIFKeySet(psGTIF, ProjCoordTransGeoKey, TYPE_SHORT, 1,
2708
0
                   CT_AzimuthalEquidistant);
2709
2710
0
        GTIFKeySet(psGTIF, ProjCenterLatGeoKey, TYPE_DOUBLE, 1,
2711
0
                   ConvertAngularParam(poSRSCompatibleOfWKT1->GetNormProjParm(
2712
0
                       SRS_PP_LATITUDE_OF_CENTER, 0.0)));
2713
2714
0
        GTIFKeySet(psGTIF, ProjCenterLongGeoKey, TYPE_DOUBLE, 1,
2715
0
                   ConvertAngularParam(poSRSCompatibleOfWKT1->GetNormProjParm(
2716
0
                       SRS_PP_LONGITUDE_OF_CENTER, 0.0)));
2717
2718
0
        GTIFKeySet(
2719
0
            psGTIF, ProjFalseEastingGeoKey, TYPE_DOUBLE, 1,
2720
0
            poSRSCompatibleOfWKT1->GetProjParm(SRS_PP_FALSE_EASTING, 0.0));
2721
2722
0
        GTIFKeySet(
2723
0
            psGTIF, ProjFalseNorthingGeoKey, TYPE_DOUBLE, 1,
2724
0
            poSRSCompatibleOfWKT1->GetProjParm(SRS_PP_FALSE_NORTHING, 0.0));
2725
0
    }
2726
2727
0
    else if (EQUAL(pszProjection, SRS_PT_MILLER_CYLINDRICAL))
2728
0
    {
2729
0
        GTIFKeySet(psGTIF, GTModelTypeGeoKey, TYPE_SHORT, 1,
2730
0
                   ModelTypeProjected);
2731
0
        GTIFKeySet(psGTIF, ProjectedCSTypeGeoKey, TYPE_SHORT, 1, KvUserDefined);
2732
0
        GTIFKeySet(psGTIF, ProjectionGeoKey, TYPE_SHORT, 1, KvUserDefined);
2733
2734
0
        GTIFKeySet(psGTIF, ProjCoordTransGeoKey, TYPE_SHORT, 1,
2735
0
                   CT_MillerCylindrical);
2736
2737
0
        GTIFKeySet(psGTIF, ProjCenterLatGeoKey, TYPE_DOUBLE, 1,
2738
0
                   ConvertAngularParam(poSRSCompatibleOfWKT1->GetNormProjParm(
2739
0
                       SRS_PP_LATITUDE_OF_CENTER, 0.0)));
2740
2741
0
        GTIFKeySet(psGTIF, ProjCenterLongGeoKey, TYPE_DOUBLE, 1,
2742
0
                   ConvertAngularParam(poSRSCompatibleOfWKT1->GetNormProjParm(
2743
0
                       SRS_PP_LONGITUDE_OF_CENTER, 0.0)));
2744
2745
0
        GTIFKeySet(
2746
0
            psGTIF, ProjFalseEastingGeoKey, TYPE_DOUBLE, 1,
2747
0
            poSRSCompatibleOfWKT1->GetProjParm(SRS_PP_FALSE_EASTING, 0.0));
2748
2749
0
        GTIFKeySet(
2750
0
            psGTIF, ProjFalseNorthingGeoKey, TYPE_DOUBLE, 1,
2751
0
            poSRSCompatibleOfWKT1->GetProjParm(SRS_PP_FALSE_NORTHING, 0.0));
2752
0
    }
2753
2754
0
    else if (EQUAL(pszProjection, SRS_PT_EQUIRECTANGULAR))
2755
0
    {
2756
0
        GTIFKeySet(psGTIF, GTModelTypeGeoKey, TYPE_SHORT, 1,
2757
0
                   ModelTypeProjected);
2758
0
        GTIFKeySet(psGTIF, ProjectedCSTypeGeoKey, TYPE_SHORT, 1, KvUserDefined);
2759
0
        GTIFKeySet(psGTIF, ProjectionGeoKey, TYPE_SHORT, 1, KvUserDefined);
2760
2761
0
        GTIFKeySet(psGTIF, ProjCoordTransGeoKey, TYPE_SHORT, 1,
2762
0
                   CT_Equirectangular);
2763
2764
0
        GTIFKeySet(psGTIF, ProjCenterLatGeoKey, TYPE_DOUBLE, 1,
2765
0
                   ConvertAngularParam(poSRSCompatibleOfWKT1->GetNormProjParm(
2766
0
                       SRS_PP_LATITUDE_OF_ORIGIN, 0.0)));
2767
2768
0
        GTIFKeySet(psGTIF, ProjCenterLongGeoKey, TYPE_DOUBLE, 1,
2769
0
                   ConvertAngularParam(poSRSCompatibleOfWKT1->GetNormProjParm(
2770
0
                       SRS_PP_CENTRAL_MERIDIAN, 0.0)));
2771
2772
0
        GTIFKeySet(psGTIF, ProjStdParallel1GeoKey, TYPE_DOUBLE, 1,
2773
0
                   ConvertAngularParam(poSRSCompatibleOfWKT1->GetNormProjParm(
2774
0
                       SRS_PP_STANDARD_PARALLEL_1, 0.0)));
2775
2776
0
        GTIFKeySet(
2777
0
            psGTIF, ProjFalseEastingGeoKey, TYPE_DOUBLE, 1,
2778
0
            poSRSCompatibleOfWKT1->GetProjParm(SRS_PP_FALSE_EASTING, 0.0));
2779
2780
0
        GTIFKeySet(
2781
0
            psGTIF, ProjFalseNorthingGeoKey, TYPE_DOUBLE, 1,
2782
0
            poSRSCompatibleOfWKT1->GetProjParm(SRS_PP_FALSE_NORTHING, 0.0));
2783
0
    }
2784
2785
0
    else if (EQUAL(pszProjection, SRS_PT_GNOMONIC))
2786
0
    {
2787
0
        GTIFKeySet(psGTIF, GTModelTypeGeoKey, TYPE_SHORT, 1,
2788
0
                   ModelTypeProjected);
2789
0
        GTIFKeySet(psGTIF, ProjectedCSTypeGeoKey, TYPE_SHORT, 1, KvUserDefined);
2790
0
        GTIFKeySet(psGTIF, ProjectionGeoKey, TYPE_SHORT, 1, KvUserDefined);
2791
2792
0
        GTIFKeySet(psGTIF, ProjCoordTransGeoKey, TYPE_SHORT, 1, CT_Gnomonic);
2793
2794
0
        GTIFKeySet(psGTIF, ProjCenterLatGeoKey, TYPE_DOUBLE, 1,
2795
0
                   ConvertAngularParam(poSRSCompatibleOfWKT1->GetNormProjParm(
2796
0
                       SRS_PP_LATITUDE_OF_ORIGIN, 0.0)));
2797
2798
0
        GTIFKeySet(psGTIF, ProjCenterLongGeoKey, TYPE_DOUBLE, 1,
2799
0
                   ConvertAngularParam(poSRSCompatibleOfWKT1->GetNormProjParm(
2800
0
                       SRS_PP_CENTRAL_MERIDIAN, 0.0)));
2801
2802
0
        GTIFKeySet(
2803
0
            psGTIF, ProjFalseEastingGeoKey, TYPE_DOUBLE, 1,
2804
0
            poSRSCompatibleOfWKT1->GetProjParm(SRS_PP_FALSE_EASTING, 0.0));
2805
2806
0
        GTIFKeySet(
2807
0
            psGTIF, ProjFalseNorthingGeoKey, TYPE_DOUBLE, 1,
2808
0
            poSRSCompatibleOfWKT1->GetProjParm(SRS_PP_FALSE_NORTHING, 0.0));
2809
0
    }
2810
2811
0
    else if (EQUAL(pszProjection, SRS_PT_LAMBERT_AZIMUTHAL_EQUAL_AREA))
2812
0
    {
2813
0
        GTIFKeySet(psGTIF, GTModelTypeGeoKey, TYPE_SHORT, 1,
2814
0
                   ModelTypeProjected);
2815
0
        GTIFKeySet(psGTIF, ProjectedCSTypeGeoKey, TYPE_SHORT, 1, KvUserDefined);
2816
0
        GTIFKeySet(psGTIF, ProjectionGeoKey, TYPE_SHORT, 1, KvUserDefined);
2817
2818
0
        GTIFKeySet(psGTIF, ProjCoordTransGeoKey, TYPE_SHORT, 1,
2819
0
                   CT_LambertAzimEqualArea);
2820
2821
0
        GTIFKeySet(psGTIF, ProjCenterLatGeoKey, TYPE_DOUBLE, 1,
2822
0
                   ConvertAngularParam(poSRSCompatibleOfWKT1->GetNormProjParm(
2823
0
                       SRS_PP_LATITUDE_OF_CENTER, 0.0)));
2824
2825
0
        GTIFKeySet(psGTIF, ProjCenterLongGeoKey, TYPE_DOUBLE, 1,
2826
0
                   ConvertAngularParam(poSRSCompatibleOfWKT1->GetNormProjParm(
2827
0
                       SRS_PP_LONGITUDE_OF_CENTER, 0.0)));
2828
2829
0
        GTIFKeySet(
2830
0
            psGTIF, ProjFalseEastingGeoKey, TYPE_DOUBLE, 1,
2831
0
            poSRSCompatibleOfWKT1->GetProjParm(SRS_PP_FALSE_EASTING, 0.0));
2832
2833
0
        GTIFKeySet(
2834
0
            psGTIF, ProjFalseNorthingGeoKey, TYPE_DOUBLE, 1,
2835
0
            poSRSCompatibleOfWKT1->GetProjParm(SRS_PP_FALSE_NORTHING, 0.0));
2836
0
    }
2837
2838
0
    else if (EQUAL(pszProjection, SRS_PT_ORTHOGRAPHIC))
2839
0
    {
2840
0
        GTIFKeySet(psGTIF, GTModelTypeGeoKey, TYPE_SHORT, 1,
2841
0
                   ModelTypeProjected);
2842
0
        GTIFKeySet(psGTIF, ProjectedCSTypeGeoKey, TYPE_SHORT, 1, KvUserDefined);
2843
0
        GTIFKeySet(psGTIF, ProjectionGeoKey, TYPE_SHORT, 1, KvUserDefined);
2844
2845
0
        GTIFKeySet(psGTIF, ProjCoordTransGeoKey, TYPE_SHORT, 1,
2846
0
                   CT_Orthographic);
2847
2848
0
        GTIFKeySet(psGTIF, ProjCenterLatGeoKey, TYPE_DOUBLE, 1,
2849
0
                   ConvertAngularParam(poSRSCompatibleOfWKT1->GetNormProjParm(
2850
0
                       SRS_PP_LATITUDE_OF_ORIGIN, 0.0)));
2851
2852
0
        GTIFKeySet(psGTIF, ProjCenterLongGeoKey, TYPE_DOUBLE, 1,
2853
0
                   ConvertAngularParam(poSRSCompatibleOfWKT1->GetNormProjParm(
2854
0
                       SRS_PP_CENTRAL_MERIDIAN, 0.0)));
2855
2856
0
        GTIFKeySet(
2857
0
            psGTIF, ProjFalseEastingGeoKey, TYPE_DOUBLE, 1,
2858
0
            poSRSCompatibleOfWKT1->GetProjParm(SRS_PP_FALSE_EASTING, 0.0));
2859
2860
0
        GTIFKeySet(
2861
0
            psGTIF, ProjFalseNorthingGeoKey, TYPE_DOUBLE, 1,
2862
0
            poSRSCompatibleOfWKT1->GetProjParm(SRS_PP_FALSE_NORTHING, 0.0));
2863
0
    }
2864
2865
0
    else if (EQUAL(pszProjection, SRS_PT_NEW_ZEALAND_MAP_GRID))
2866
0
    {
2867
0
        GTIFKeySet(psGTIF, GTModelTypeGeoKey, TYPE_SHORT, 1,
2868
0
                   ModelTypeProjected);
2869
0
        GTIFKeySet(psGTIF, ProjectedCSTypeGeoKey, TYPE_SHORT, 1, KvUserDefined);
2870
0
        GTIFKeySet(psGTIF, ProjectionGeoKey, TYPE_SHORT, 1, KvUserDefined);
2871
2872
0
        GTIFKeySet(psGTIF, ProjCoordTransGeoKey, TYPE_SHORT, 1,
2873
0
                   CT_NewZealandMapGrid);
2874
2875
0
        GTIFKeySet(psGTIF, ProjNatOriginLatGeoKey, TYPE_DOUBLE, 1,
2876
0
                   ConvertAngularParam(poSRSCompatibleOfWKT1->GetNormProjParm(
2877
0
                       SRS_PP_LATITUDE_OF_ORIGIN, 0.0)));
2878
2879
0
        GTIFKeySet(psGTIF, ProjNatOriginLongGeoKey, TYPE_DOUBLE, 1,
2880
0
                   ConvertAngularParam(poSRSCompatibleOfWKT1->GetNormProjParm(
2881
0
                       SRS_PP_CENTRAL_MERIDIAN, 0.0)));
2882
2883
0
        GTIFKeySet(
2884
0
            psGTIF, ProjFalseEastingGeoKey, TYPE_DOUBLE, 1,
2885
0
            poSRSCompatibleOfWKT1->GetProjParm(SRS_PP_FALSE_EASTING, 0.0));
2886
2887
0
        GTIFKeySet(
2888
0
            psGTIF, ProjFalseNorthingGeoKey, TYPE_DOUBLE, 1,
2889
0
            poSRSCompatibleOfWKT1->GetProjParm(SRS_PP_FALSE_NORTHING, 0.0));
2890
0
    }
2891
2892
0
    else if (EQUAL(pszProjection, SRS_PT_ROBINSON))
2893
0
    {
2894
0
        GTIFKeySet(psGTIF, GTModelTypeGeoKey, TYPE_SHORT, 1,
2895
0
                   ModelTypeProjected);
2896
0
        GTIFKeySet(psGTIF, ProjectedCSTypeGeoKey, TYPE_SHORT, 1, KvUserDefined);
2897
0
        GTIFKeySet(psGTIF, ProjectionGeoKey, TYPE_SHORT, 1, KvUserDefined);
2898
2899
0
        GTIFKeySet(psGTIF, ProjCoordTransGeoKey, TYPE_SHORT, 1, CT_Robinson);
2900
2901
0
        GTIFKeySet(psGTIF, ProjCenterLongGeoKey, TYPE_DOUBLE, 1,
2902
0
                   ConvertAngularParam(poSRSCompatibleOfWKT1->GetNormProjParm(
2903
0
                       SRS_PP_LONGITUDE_OF_CENTER, 0.0)));
2904
2905
0
        GTIFKeySet(
2906
0
            psGTIF, ProjFalseEastingGeoKey, TYPE_DOUBLE, 1,
2907
0
            poSRSCompatibleOfWKT1->GetProjParm(SRS_PP_FALSE_EASTING, 0.0));
2908
2909
0
        GTIFKeySet(
2910
0
            psGTIF, ProjFalseNorthingGeoKey, TYPE_DOUBLE, 1,
2911
0
            poSRSCompatibleOfWKT1->GetProjParm(SRS_PP_FALSE_NORTHING, 0.0));
2912
0
    }
2913
2914
0
    else if (EQUAL(pszProjection, SRS_PT_SINUSOIDAL))
2915
0
    {
2916
0
        GTIFKeySet(psGTIF, GTModelTypeGeoKey, TYPE_SHORT, 1,
2917
0
                   ModelTypeProjected);
2918
0
        GTIFKeySet(psGTIF, ProjectedCSTypeGeoKey, TYPE_SHORT, 1, KvUserDefined);
2919
0
        GTIFKeySet(psGTIF, ProjectionGeoKey, TYPE_SHORT, 1, KvUserDefined);
2920
2921
0
        GTIFKeySet(psGTIF, ProjCoordTransGeoKey, TYPE_SHORT, 1, CT_Sinusoidal);
2922
2923
0
        GTIFKeySet(psGTIF, ProjCenterLongGeoKey, TYPE_DOUBLE, 1,
2924
0
                   ConvertAngularParam(poSRSCompatibleOfWKT1->GetNormProjParm(
2925
0
                       SRS_PP_LONGITUDE_OF_CENTER, 0.0)));
2926
2927
0
        GTIFKeySet(
2928
0
            psGTIF, ProjFalseEastingGeoKey, TYPE_DOUBLE, 1,
2929
0
            poSRSCompatibleOfWKT1->GetProjParm(SRS_PP_FALSE_EASTING, 0.0));
2930
2931
0
        GTIFKeySet(
2932
0
            psGTIF, ProjFalseNorthingGeoKey, TYPE_DOUBLE, 1,
2933
0
            poSRSCompatibleOfWKT1->GetProjParm(SRS_PP_FALSE_NORTHING, 0.0));
2934
0
    }
2935
2936
0
    else if (EQUAL(pszProjection, SRS_PT_VANDERGRINTEN))
2937
0
    {
2938
0
        GTIFKeySet(psGTIF, GTModelTypeGeoKey, TYPE_SHORT, 1,
2939
0
                   ModelTypeProjected);
2940
0
        GTIFKeySet(psGTIF, ProjectedCSTypeGeoKey, TYPE_SHORT, 1, KvUserDefined);
2941
0
        GTIFKeySet(psGTIF, ProjectionGeoKey, TYPE_SHORT, 1, KvUserDefined);
2942
2943
0
        GTIFKeySet(psGTIF, ProjCoordTransGeoKey, TYPE_SHORT, 1,
2944
0
                   CT_VanDerGrinten);
2945
2946
0
        GTIFKeySet(psGTIF, ProjCenterLongGeoKey, TYPE_DOUBLE, 1,
2947
0
                   ConvertAngularParam(poSRSCompatibleOfWKT1->GetNormProjParm(
2948
0
                       SRS_PP_CENTRAL_MERIDIAN, 0.0)));
2949
2950
0
        GTIFKeySet(
2951
0
            psGTIF, ProjFalseEastingGeoKey, TYPE_DOUBLE, 1,
2952
0
            poSRSCompatibleOfWKT1->GetProjParm(SRS_PP_FALSE_EASTING, 0.0));
2953
2954
0
        GTIFKeySet(
2955
0
            psGTIF, ProjFalseNorthingGeoKey, TYPE_DOUBLE, 1,
2956
0
            poSRSCompatibleOfWKT1->GetProjParm(SRS_PP_FALSE_NORTHING, 0.0));
2957
0
    }
2958
2959
0
    else if (EQUAL(pszProjection, SRS_PT_LAMBERT_CONFORMAL_CONIC_2SP))
2960
0
    {
2961
0
        GTIFKeySet(psGTIF, GTModelTypeGeoKey, TYPE_SHORT, 1,
2962
0
                   ModelTypeProjected);
2963
0
        GTIFKeySet(psGTIF, ProjectedCSTypeGeoKey, TYPE_SHORT, 1, KvUserDefined);
2964
0
        GTIFKeySet(psGTIF, ProjectionGeoKey, TYPE_SHORT, 1, KvUserDefined);
2965
2966
0
        GTIFKeySet(psGTIF, ProjCoordTransGeoKey, TYPE_SHORT, 1,
2967
0
                   CT_LambertConfConic_2SP);
2968
2969
0
        GTIFKeySet(psGTIF, ProjFalseOriginLatGeoKey, TYPE_DOUBLE, 1,
2970
0
                   ConvertAngularParam(poSRSCompatibleOfWKT1->GetNormProjParm(
2971
0
                       SRS_PP_LATITUDE_OF_ORIGIN, 0.0)));
2972
2973
0
        GTIFKeySet(psGTIF, ProjFalseOriginLongGeoKey, TYPE_DOUBLE, 1,
2974
0
                   ConvertAngularParam(poSRSCompatibleOfWKT1->GetNormProjParm(
2975
0
                       SRS_PP_CENTRAL_MERIDIAN, 0.0)));
2976
2977
0
        GTIFKeySet(psGTIF, ProjStdParallel1GeoKey, TYPE_DOUBLE, 1,
2978
0
                   ConvertAngularParam(poSRSCompatibleOfWKT1->GetNormProjParm(
2979
0
                       SRS_PP_STANDARD_PARALLEL_1, 0.0)));
2980
2981
0
        GTIFKeySet(psGTIF, ProjStdParallel2GeoKey, TYPE_DOUBLE, 1,
2982
0
                   ConvertAngularParam(poSRSCompatibleOfWKT1->GetNormProjParm(
2983
0
                       SRS_PP_STANDARD_PARALLEL_2, 0.0)));
2984
2985
0
        GTIFKeySet(
2986
0
            psGTIF, ProjFalseOriginEastingGeoKey, TYPE_DOUBLE, 1,
2987
0
            poSRSCompatibleOfWKT1->GetProjParm(SRS_PP_FALSE_EASTING, 0.0));
2988
2989
0
        GTIFKeySet(
2990
0
            psGTIF, ProjFalseOriginNorthingGeoKey, TYPE_DOUBLE, 1,
2991
0
            poSRSCompatibleOfWKT1->GetProjParm(SRS_PP_FALSE_NORTHING, 0.0));
2992
0
    }
2993
2994
0
    else if (EQUAL(pszProjection, SRS_PT_LAMBERT_CONFORMAL_CONIC_1SP))
2995
0
    {
2996
0
        GTIFKeySet(psGTIF, GTModelTypeGeoKey, TYPE_SHORT, 1,
2997
0
                   ModelTypeProjected);
2998
0
        GTIFKeySet(psGTIF, ProjectedCSTypeGeoKey, TYPE_SHORT, 1, KvUserDefined);
2999
0
        GTIFKeySet(psGTIF, ProjectionGeoKey, TYPE_SHORT, 1, KvUserDefined);
3000
3001
0
        GTIFKeySet(psGTIF, ProjCoordTransGeoKey, TYPE_SHORT, 1,
3002
0
                   CT_LambertConfConic_1SP);
3003
3004
0
        GTIFKeySet(psGTIF, ProjNatOriginLatGeoKey, TYPE_DOUBLE, 1,
3005
0
                   ConvertAngularParam(poSRSCompatibleOfWKT1->GetNormProjParm(
3006
0
                       SRS_PP_LATITUDE_OF_ORIGIN, 0.0)));
3007
3008
0
        GTIFKeySet(psGTIF, ProjNatOriginLongGeoKey, TYPE_DOUBLE, 1,
3009
0
                   ConvertAngularParam(poSRSCompatibleOfWKT1->GetNormProjParm(
3010
0
                       SRS_PP_CENTRAL_MERIDIAN, 0.0)));
3011
3012
0
        GTIFKeySet(
3013
0
            psGTIF, ProjScaleAtNatOriginGeoKey, TYPE_DOUBLE, 1,
3014
0
            poSRSCompatibleOfWKT1->GetNormProjParm(SRS_PP_SCALE_FACTOR, 1.0));
3015
3016
0
        GTIFKeySet(
3017
0
            psGTIF, ProjFalseEastingGeoKey, TYPE_DOUBLE, 1,
3018
0
            poSRSCompatibleOfWKT1->GetProjParm(SRS_PP_FALSE_EASTING, 0.0));
3019
3020
0
        GTIFKeySet(
3021
0
            psGTIF, ProjFalseNorthingGeoKey, TYPE_DOUBLE, 1,
3022
0
            poSRSCompatibleOfWKT1->GetProjParm(SRS_PP_FALSE_NORTHING, 0.0));
3023
0
    }
3024
3025
0
    else if (EQUAL(pszProjection, SRS_PT_CYLINDRICAL_EQUAL_AREA))
3026
0
    {
3027
0
        GTIFKeySet(psGTIF, GTModelTypeGeoKey, TYPE_SHORT, 1,
3028
0
                   ModelTypeProjected);
3029
0
        GTIFKeySet(psGTIF, ProjectedCSTypeGeoKey, TYPE_SHORT, 1, KvUserDefined);
3030
0
        GTIFKeySet(psGTIF, ProjectionGeoKey, TYPE_SHORT, 1, KvUserDefined);
3031
3032
0
        GTIFKeySet(psGTIF, ProjCoordTransGeoKey, TYPE_SHORT, 1,
3033
0
                   CT_CylindricalEqualArea);
3034
3035
0
        GTIFKeySet(psGTIF, ProjNatOriginLongGeoKey, TYPE_DOUBLE, 1,
3036
0
                   ConvertAngularParam(poSRSCompatibleOfWKT1->GetNormProjParm(
3037
0
                       SRS_PP_CENTRAL_MERIDIAN, 0.0)));
3038
3039
0
        GTIFKeySet(psGTIF, ProjStdParallel1GeoKey, TYPE_DOUBLE, 1,
3040
0
                   ConvertAngularParam(poSRSCompatibleOfWKT1->GetNormProjParm(
3041
0
                       SRS_PP_STANDARD_PARALLEL_1, 0.0)));
3042
3043
0
        GTIFKeySet(
3044
0
            psGTIF, ProjFalseEastingGeoKey, TYPE_DOUBLE, 1,
3045
0
            poSRSCompatibleOfWKT1->GetProjParm(SRS_PP_FALSE_EASTING, 0.0));
3046
3047
0
        GTIFKeySet(
3048
0
            psGTIF, ProjFalseNorthingGeoKey, TYPE_DOUBLE, 1,
3049
0
            poSRSCompatibleOfWKT1->GetProjParm(SRS_PP_FALSE_NORTHING, 0.0));
3050
0
    }
3051
0
    else
3052
0
    {
3053
0
        bWritePEString = true;
3054
0
        bUnknownProjection = true;
3055
0
    }
3056
3057
    // Note that VERTCS is an ESRI "spelling" of VERT_CS so we assume if
3058
    // we find it that we should try to treat this as a PE string.
3059
0
    if (eFlavor == GEOTIFF_KEYS_ESRI_PE ||
3060
0
        poSRS->GetAttrValue("VERTCS") != nullptr)
3061
0
    {
3062
0
        bWritePEString = true;
3063
0
    }
3064
3065
0
    if (nPCS == KvUserDefined)
3066
0
    {
3067
0
        const char *pszPROJ4Ext =
3068
0
            poSRS->GetExtension("PROJCS", "PROJ4", nullptr);
3069
0
        if (pszPROJ4Ext &&
3070
0
            strstr(pszPROJ4Ext, "+proj=merc +a=6378137 +b=6378137"))
3071
0
        {
3072
0
            bWritePEString = true;
3073
0
        }
3074
0
    }
3075
3076
0
    bWritePEString &=
3077
0
        CPLTestBool(CPLGetConfigOption("GTIFF_ESRI_CITATION", "YES"));
3078
3079
0
    bool peStrStored = false;
3080
3081
0
    if (bWritePEString)
3082
0
    {
3083
        // Anything we can't map, store as an ESRI PE string with a citation
3084
        // key.
3085
0
        char *pszPEString = nullptr;
3086
        // We cheat a bit, but if we have a custom_proj4, do not morph to ESRI
3087
        // so as to keep the EXTENSION PROJ4 node
3088
0
        const char *const apszOptionsDefault[] = {nullptr};
3089
0
        const char *const apszOptionsEsri[] = {"FORMAT=WKT1_ESRI", nullptr};
3090
0
        const char *const *papszOptions = apszOptionsDefault;
3091
0
        if (!(bUnknownProjection &&
3092
0
              poSRS->GetExtension("PROJCS", "PROJ4", nullptr) != nullptr) &&
3093
0
            !(poSRS->IsProjected() && !poSRS->IsCompound() &&
3094
0
              poSRS->GetAxesCount() == 3))
3095
0
        {
3096
0
            papszOptions = apszOptionsEsri;
3097
0
        }
3098
0
        poSRS->exportToWkt(&pszPEString, papszOptions);
3099
0
        const int peStrLen = static_cast<int>(strlen(pszPEString));
3100
0
        if (peStrLen > 0)
3101
0
        {
3102
0
            char *outPeStr = static_cast<char *>(
3103
0
                CPLMalloc(peStrLen + strlen("ESRI PE String = ") + 1));
3104
0
            strcpy(outPeStr, "ESRI PE String = ");
3105
0
            strcat(outPeStr, pszPEString);
3106
0
            oMapAsciiKeys[PCSCitationGeoKey] = outPeStr;
3107
0
            peStrStored = true;
3108
0
            CPLFree(outPeStr);
3109
0
        }
3110
0
        CPLFree(pszPEString);
3111
0
        GTIFKeySet(psGTIF, GTModelTypeGeoKey, TYPE_SHORT, 1, KvUserDefined);
3112
3113
        // Not completely sure we really need to imitate ArcGIS to that point
3114
        // but that cannot hurt.
3115
0
        if (nPCS == 3857)
3116
0
        {
3117
0
            oMapAsciiKeys[GTCitationGeoKey] =
3118
0
                "PCS Name = WGS_1984_Web_Mercator_Auxiliary_Sphere";
3119
0
            GTIFKeySet(psGTIF, GeographicTypeGeoKey, TYPE_SHORT, 1, GCS_WGS_84);
3120
0
            GTIFKeySet(psGTIF, GeogSemiMajorAxisGeoKey, TYPE_DOUBLE, 1,
3121
0
                       6378137.0);
3122
0
            GTIFKeySet(psGTIF, GeogInvFlatteningGeoKey, TYPE_DOUBLE, 1,
3123
0
                       298.257223563);
3124
0
        }
3125
0
    }
3126
3127
    /* -------------------------------------------------------------------- */
3128
    /*      Is there a false easting/northing set?  If so, write out a      */
3129
    /*      special geokey tag to indicate that GDAL has written these      */
3130
    /*      with the proper interpretation of the linear units.             */
3131
    /* -------------------------------------------------------------------- */
3132
0
    double dfFE = 0.0;
3133
0
    double dfFN = 0.0;
3134
3135
0
    if (eVersion == GEOTIFF_VERSION_1_0 &&
3136
0
        (GDALGTIFKeyGetDOUBLE(psGTIF, ProjFalseEastingGeoKey, &dfFE, 0, 1) ||
3137
0
         GDALGTIFKeyGetDOUBLE(psGTIF, ProjFalseNorthingGeoKey, &dfFN, 0, 1) ||
3138
0
         GDALGTIFKeyGetDOUBLE(psGTIF, ProjFalseOriginEastingGeoKey, &dfFE, 0,
3139
0
                              1) ||
3140
0
         GDALGTIFKeyGetDOUBLE(psGTIF, ProjFalseOriginNorthingGeoKey, &dfFN, 0,
3141
0
                              1)) &&
3142
0
        (dfFE != 0.0 || dfFN != 0.0) && nUOMLengthCode != 9001)
3143
0
    {
3144
0
        GTIFKeySet(psGTIF, ProjLinearUnitsInterpCorrectGeoKey, TYPE_SHORT, 1,
3145
0
                   static_cast<short>(1));
3146
0
    }
3147
3148
    /* -------------------------------------------------------------------- */
3149
    /*      Write linear units information.                                 */
3150
    /* -------------------------------------------------------------------- */
3151
0
    if (poSRS->IsGeocentric())
3152
0
    {
3153
0
        GTIFKeySet(psGTIF, GeogLinearUnitsGeoKey, TYPE_SHORT, 1,
3154
0
                   nUOMLengthCode);
3155
0
        if (nUOMLengthCode == KvUserDefined)
3156
0
            GTIFKeySet(psGTIF, GeogLinearUnitSizeGeoKey, TYPE_DOUBLE, 1,
3157
0
                       dfLinearUOM);
3158
0
    }
3159
0
    else if (!poSRS->IsGeographic() &&
3160
0
             (nPCS == KvUserDefined || eVersion == GEOTIFF_VERSION_1_0))
3161
0
    {
3162
0
        GTIFKeySet(psGTIF, ProjLinearUnitsGeoKey, TYPE_SHORT, 1,
3163
0
                   nUOMLengthCode);
3164
0
        if (nUOMLengthCode == KvUserDefined)
3165
0
            GTIFKeySet(psGTIF, ProjLinearUnitSizeGeoKey, TYPE_DOUBLE, 1,
3166
0
                       dfLinearUOM);
3167
3168
        // If linear units name is available and user defined, store it as
3169
        // citation.
3170
0
        if (!peStrStored && nUOMLengthCode == KvUserDefined &&
3171
0
            !osLinearUOMName.empty() &&
3172
0
            CPLTestBool(CPLGetConfigOption("GTIFF_ESRI_CITATION", "YES")))
3173
0
        {
3174
0
            SetLinearUnitCitation(oMapAsciiKeys, osLinearUOMName.c_str());
3175
0
        }
3176
0
    }
3177
3178
    /* -------------------------------------------------------------------- */
3179
    /*      Write angular units.                                            */
3180
    /* -------------------------------------------------------------------- */
3181
3182
0
    if (bHasEllipsoid &&
3183
0
        (nGCS == KvUserDefined || eVersion == GEOTIFF_VERSION_1_0))
3184
0
    {
3185
0
        if (EQUAL(osAngUnitName.c_str(), "Degree"))
3186
0
            GTIFKeySet(psGTIF, GeogAngularUnitsGeoKey, TYPE_SHORT, 1,
3187
0
                       Angular_Degree);
3188
0
        else if (EQUAL(osAngUnitName.c_str(), "arc-second"))
3189
0
            GTIFKeySet(psGTIF, GeogAngularUnitsGeoKey, TYPE_SHORT, 1,
3190
0
                       Angular_Arc_Second);
3191
0
        else if (EQUAL(osAngUnitName.c_str(), "arc-minute"))
3192
0
            GTIFKeySet(psGTIF, GeogAngularUnitsGeoKey, TYPE_SHORT, 1,
3193
0
                       Angular_Arc_Minute);
3194
0
        else if (EQUAL(osAngUnitName.c_str(), "grad"))
3195
0
            GTIFKeySet(psGTIF, GeogAngularUnitsGeoKey, TYPE_SHORT, 1,
3196
0
                       Angular_Grad);
3197
0
        else if (EQUAL(osAngUnitName.c_str(), "gon"))
3198
0
            GTIFKeySet(psGTIF, GeogAngularUnitsGeoKey, TYPE_SHORT, 1,
3199
0
                       Angular_Gon);
3200
0
        else if (EQUAL(osAngUnitName.c_str(), "radian"))
3201
0
            GTIFKeySet(psGTIF, GeogAngularUnitsGeoKey, TYPE_SHORT, 1,
3202
0
                       Angular_Radian);
3203
        // else if (EQUAL(osAngUnitName.c_str(), "microradian"))
3204
        //    GTIFKeySet(psGTIF, GeogAngularUnitsGeoKey, TYPE_SHORT, 1,
3205
        //               9109);
3206
0
        else
3207
0
        {
3208
            // GeogCitationGeoKey may be rewritten if the gcs is user defined.
3209
0
            oMapAsciiKeys[GeogCitationGeoKey] = osAngUnitName;
3210
0
            GTIFKeySet(psGTIF, GeogAngularUnitSizeGeoKey, TYPE_DOUBLE, 1,
3211
0
                       dfAngUnitValue);
3212
0
        }
3213
0
    }
3214
3215
    /* -------------------------------------------------------------------- */
3216
    /*      Try to write a citation from the main coordinate system         */
3217
    /*      name.                                                           */
3218
    /* -------------------------------------------------------------------- */
3219
0
    if (poSRS->GetName() != nullptr &&
3220
0
        ((poSRS->IsProjected() &&
3221
0
          (nPCS == KvUserDefined || eVersion == GEOTIFF_VERSION_1_0)) ||
3222
0
         poSRS->IsCompound() || poSRS->IsLocal() ||
3223
0
         (poSRS->IsGeocentric() &&
3224
0
          (nGCS == KvUserDefined || eVersion == GEOTIFF_VERSION_1_0))))
3225
0
    {
3226
0
        if (!(bWritePEString && nPCS == 3857))
3227
0
        {
3228
0
            oMapAsciiKeys[GTCitationGeoKey] = poSRS->GetName();
3229
0
        }
3230
0
    }
3231
3232
    /* -------------------------------------------------------------------- */
3233
    /*      Try to write a GCS citation.                                    */
3234
    /* -------------------------------------------------------------------- */
3235
0
    if (nGCS == KvUserDefined || eVersion == GEOTIFF_VERSION_1_0)
3236
0
    {
3237
0
        const OGR_SRSNode *poGCS = poSRS->GetAttrNode("GEOGCS");
3238
3239
0
        if (poGCS != nullptr && poGCS->GetChild(0) != nullptr)
3240
0
        {
3241
0
            oMapAsciiKeys[GeogCitationGeoKey] = poGCS->GetChild(0)->GetValue();
3242
0
        }
3243
0
    }
3244
3245
    /* -------------------------------------------------------------------- */
3246
    /*      Try to identify the GCS/datum, scanning the EPSG datum file for */
3247
    /*      a match.                                                        */
3248
    /* -------------------------------------------------------------------- */
3249
0
    if (nPCS == KvUserDefined)
3250
0
    {
3251
0
        if (nGCS == KvUserDefined && poSRS->IsGeographic() &&
3252
0
            std::fabs(poSRS->GetAngularUnits() - CPLAtof(SRS_UA_DEGREE_CONV)) <
3253
0
                1e-9)
3254
0
        {
3255
0
            if (nDatum == Datum_North_American_Datum_1927)
3256
0
                nGCS = GCS_NAD27;
3257
0
            else if (nDatum == Datum_North_American_Datum_1983)
3258
0
                nGCS = GCS_NAD83;
3259
0
            else if (nDatum == Datum_WGS84 || nDatum == DatumE_WGS84)
3260
0
                nGCS = GCS_WGS_84;
3261
0
        }
3262
3263
0
        if (nGCS != KvUserDefined)
3264
0
        {
3265
0
            GTIFKeySet(psGTIF, GeographicTypeGeoKey, TYPE_SHORT, 1, nGCS);
3266
0
        }
3267
0
        else if (nDatum != KvUserDefined)
3268
0
        {
3269
0
            GTIFKeySet(psGTIF, GeographicTypeGeoKey, TYPE_SHORT, 1,
3270
0
                       KvUserDefined);
3271
0
            GTIFKeySet(psGTIF, GeogGeodeticDatumGeoKey, TYPE_SHORT, 1, nDatum);
3272
0
        }
3273
0
        else if (nSpheroid != KvUserDefined)
3274
0
        {
3275
0
            GTIFKeySet(psGTIF, GeographicTypeGeoKey, TYPE_SHORT, 1,
3276
0
                       KvUserDefined);
3277
0
            GTIFKeySet(psGTIF, GeogGeodeticDatumGeoKey, TYPE_SHORT, 1,
3278
0
                       KvUserDefined);
3279
0
            GTIFKeySet(psGTIF, GeogEllipsoidGeoKey, TYPE_SHORT, 1, nSpheroid);
3280
0
        }
3281
0
        else if (dfSemiMajor != 0.0)
3282
0
        {
3283
0
            GTIFKeySet(psGTIF, GeographicTypeGeoKey, TYPE_SHORT, 1,
3284
0
                       KvUserDefined);
3285
0
            GTIFKeySet(psGTIF, GeogGeodeticDatumGeoKey, TYPE_SHORT, 1,
3286
0
                       KvUserDefined);
3287
0
            GTIFKeySet(psGTIF, GeogEllipsoidGeoKey, TYPE_SHORT, 1,
3288
0
                       KvUserDefined);
3289
0
            GTIFKeySet(psGTIF, GeogSemiMajorAxisGeoKey, TYPE_DOUBLE, 1,
3290
0
                       dfSemiMajor);
3291
0
            if (dfInvFlattening == 0.0)
3292
0
                GTIFKeySet(psGTIF, GeogSemiMinorAxisGeoKey, TYPE_DOUBLE, 1,
3293
0
                           dfSemiMajor);
3294
0
            else
3295
0
                GTIFKeySet(psGTIF, GeogInvFlatteningGeoKey, TYPE_DOUBLE, 1,
3296
0
                           dfInvFlattening);
3297
0
        }
3298
0
        else if (poSRS->GetAttrValue("DATUM") != nullptr &&
3299
0
                 strstr(poSRS->GetAttrValue("DATUM"), "unknown") == nullptr &&
3300
0
                 strstr(poSRS->GetAttrValue("DATUM"), "unnamed") == nullptr)
3301
3302
0
        {
3303
0
            CPLError(CE_Warning, CPLE_AppDefined,
3304
0
                     "Couldn't translate `%s' to a GeoTIFF datum.",
3305
0
                     poSRS->GetAttrValue("DATUM"));
3306
0
        }
3307
3308
0
        if (nGCS == KvUserDefined || eVersion == GEOTIFF_VERSION_1_0)
3309
0
        {
3310
            // Always set InvFlattening if it is available.
3311
            // So that it doesn't need to calculate from SemiMinor.
3312
0
            if (dfInvFlattening != 0.0)
3313
0
                GTIFKeySet(psGTIF, GeogInvFlatteningGeoKey, TYPE_DOUBLE, 1,
3314
0
                           dfInvFlattening);
3315
            // Always set SemiMajor to keep the precision and in case of
3316
            // editing.
3317
0
            if (dfSemiMajor != 0.0)
3318
0
                GTIFKeySet(psGTIF, GeogSemiMajorAxisGeoKey, TYPE_DOUBLE, 1,
3319
0
                           dfSemiMajor);
3320
3321
0
            if (nGCS == KvUserDefined &&
3322
0
                CPLTestBool(CPLGetConfigOption("GTIFF_ESRI_CITATION", "YES")))
3323
0
            {
3324
0
                SetGeogCSCitation(psGTIF, oMapAsciiKeys, poSRS,
3325
0
                                  osAngUnitName.c_str(), nDatum, nSpheroid);
3326
0
            }
3327
0
        }
3328
0
    }
3329
3330
/* -------------------------------------------------------------------- */
3331
/*      Do we have TOWGS84 parameters?                                  */
3332
/* -------------------------------------------------------------------- */
3333
0
#if !defined(GEO_NORMALIZE_DISABLE_TOWGS84)
3334
0
    double adfTOWGS84[7] = {0.0};
3335
3336
0
    if ((nGCS == KvUserDefined || eVersion == GEOTIFF_VERSION_1_0) &&
3337
0
        poSRS->GetTOWGS84(adfTOWGS84) == OGRERR_NONE)
3338
0
    {
3339
        // If we are writing a SRS with a EPSG code, and that the EPSG code
3340
        // of the current SRS object and the one coming from the EPSG code
3341
        // are the same, then by default, do not write them.
3342
0
        bool bUseReferenceTOWGS84 = false;
3343
0
        const char *pszAuthName = poSRS->GetAuthorityName(nullptr);
3344
0
        const char *pszAuthCode = poSRS->GetAuthorityCode(nullptr);
3345
0
        if (pszAuthName && EQUAL(pszAuthName, "EPSG") && pszAuthCode)
3346
0
        {
3347
0
            CPLErrorHandlerPusher oErrorHandler(CPLQuietErrorHandler);
3348
0
            OGRSpatialReference oRefSRS;
3349
0
            double adfRefTOWGS84[7] = {0.0};
3350
0
            if (oRefSRS.importFromEPSG(atoi(pszAuthCode)) == OGRERR_NONE)
3351
0
            {
3352
0
                oRefSRS.AddGuessedTOWGS84();
3353
0
                if (oRefSRS.GetTOWGS84(adfRefTOWGS84) == OGRERR_NONE &&
3354
0
                    memcmp(adfRefTOWGS84, adfTOWGS84, sizeof(adfTOWGS84)) == 0)
3355
0
                {
3356
0
                    bUseReferenceTOWGS84 = true;
3357
0
                }
3358
0
            }
3359
0
        }
3360
0
        const char *pszWriteTOWGS84 =
3361
0
            CPLGetConfigOption("GTIFF_WRITE_TOWGS84", "AUTO");
3362
0
        if ((EQUAL(pszWriteTOWGS84, "YES") || EQUAL(pszWriteTOWGS84, "TRUE") ||
3363
0
             EQUAL(pszWriteTOWGS84, "ON")) ||
3364
0
            (!bUseReferenceTOWGS84 && EQUAL(pszWriteTOWGS84, "AUTO")))
3365
0
        {
3366
0
            if (adfTOWGS84[3] == 0.0 && adfTOWGS84[4] == 0.0 &&
3367
0
                adfTOWGS84[5] == 0.0 && adfTOWGS84[6] == 0.0)
3368
0
            {
3369
0
                if (nGCS == GCS_WGS_84 && adfTOWGS84[0] == 0.0 &&
3370
0
                    adfTOWGS84[1] == 0.0 && adfTOWGS84[2] == 0.0)
3371
0
                {
3372
0
                    ;  // Do nothing.
3373
0
                }
3374
0
                else
3375
0
                    GTIFKeySet(psGTIF, GeogTOWGS84GeoKey, TYPE_DOUBLE, 3,
3376
0
                               adfTOWGS84);
3377
0
            }
3378
0
            else
3379
0
                GTIFKeySet(psGTIF, GeogTOWGS84GeoKey, TYPE_DOUBLE, 7,
3380
0
                           adfTOWGS84);
3381
0
        }
3382
0
    }
3383
0
#endif
3384
3385
    /* -------------------------------------------------------------------- */
3386
    /*      Do we have vertical information to set?                         */
3387
    /* -------------------------------------------------------------------- */
3388
0
    if (poSRS->GetAttrValue("COMPD_CS|VERT_CS") != nullptr)
3389
0
    {
3390
0
        bool bGotVertCSCode = false;
3391
0
        const char *pszVertCSCode = poSRS->GetAuthorityCode("COMPD_CS|VERT_CS");
3392
0
        const char *pszVertCSAuthName =
3393
0
            poSRS->GetAuthorityName("COMPD_CS|VERT_CS");
3394
0
        if (pszVertCSCode && pszVertCSAuthName && atoi(pszVertCSCode) &&
3395
0
            EQUAL(pszVertCSAuthName, "EPSG"))
3396
0
        {
3397
0
            bGotVertCSCode = true;
3398
0
            GTIFKeySet(psGTIF, VerticalCSTypeGeoKey, TYPE_SHORT, 1,
3399
0
                       atoi(pszVertCSCode));
3400
0
        }
3401
0
        else if (eVersion >= GEOTIFF_VERSION_1_1)
3402
0
        {
3403
0
            GTIFKeySet(psGTIF, VerticalCSTypeGeoKey, TYPE_SHORT, 1,
3404
0
                       KvUserDefined);
3405
0
        }
3406
3407
0
        if (eVersion == GEOTIFF_VERSION_1_0 || !bGotVertCSCode)
3408
0
        {
3409
0
            oMapAsciiKeys[VerticalCitationGeoKey] =
3410
0
                poSRS->GetAttrValue("COMPD_CS|VERT_CS");
3411
3412
0
            const char *pszVertDatumCode =
3413
0
                poSRS->GetAuthorityCode("COMPD_CS|VERT_CS|VERT_DATUM");
3414
0
            const char *pszVertDatumAuthName =
3415
0
                poSRS->GetAuthorityName("COMPD_CS|VERT_CS|VERT_DATUM");
3416
0
            if (pszVertDatumCode && pszVertDatumAuthName &&
3417
0
                atoi(pszVertDatumCode) && EQUAL(pszVertDatumAuthName, "EPSG"))
3418
0
            {
3419
0
                GTIFKeySet(psGTIF, VerticalDatumGeoKey, TYPE_SHORT, 1,
3420
0
                           atoi(pszVertDatumCode));
3421
0
            }
3422
0
            else if (eVersion >= GEOTIFF_VERSION_1_1)
3423
0
            {
3424
0
                GTIFKeySet(psGTIF, VerticalDatumGeoKey, TYPE_SHORT, 1,
3425
0
                           KvUserDefined);
3426
0
            }
3427
3428
0
            const char *pszVertUnitCode =
3429
0
                poSRS->GetAuthorityCode("COMPD_CS|VERT_CS|UNIT");
3430
0
            const char *pszVertUnitAuthName =
3431
0
                poSRS->GetAuthorityName("COMPD_CS|VERT_CS|UNIT");
3432
0
            if (pszVertUnitCode && pszVertUnitAuthName &&
3433
0
                atoi(pszVertUnitCode) && EQUAL(pszVertUnitAuthName, "EPSG"))
3434
0
            {
3435
0
                GTIFKeySet(psGTIF, VerticalUnitsGeoKey, TYPE_SHORT, 1,
3436
0
                           atoi(pszVertUnitCode));
3437
0
            }
3438
0
            else if (eVersion >= GEOTIFF_VERSION_1_1)
3439
0
            {
3440
0
                GTIFKeySet(psGTIF, VerticalUnitsGeoKey, TYPE_SHORT, 1,
3441
0
                           KvUserDefined);
3442
0
            }
3443
0
        }
3444
0
    }
3445
0
    else if (eVersion >= GEOTIFF_VERSION_1_1 && nVerticalCSKeyValue != 0)
3446
0
    {
3447
0
        GTIFKeySet(psGTIF, VerticalCSTypeGeoKey, TYPE_SHORT, 1,
3448
0
                   nVerticalCSKeyValue);
3449
0
    }
3450
3451
0
    const double dfCoordinateEpoch = poSRS->GetCoordinateEpoch();
3452
0
    if (dfCoordinateEpoch > 0)
3453
0
    {
3454
0
        GTIFKeySet(psGTIF, CoordinateEpochGeoKey, TYPE_DOUBLE, 1,
3455
0
                   dfCoordinateEpoch);
3456
0
    }
3457
3458
    /* -------------------------------------------------------------------- */
3459
    /*      Write all ascii keys                                            */
3460
    /* -------------------------------------------------------------------- */
3461
0
    for (const auto &oIter : oMapAsciiKeys)
3462
0
    {
3463
0
        GTIFKeySet(psGTIF, oIter.first, TYPE_ASCII, 0, oIter.second.c_str());
3464
0
    }
3465
3466
0
    return TRUE;
3467
0
}
3468
3469
/************************************************************************/
3470
/*                         GTIFWktFromMemBuf()                          */
3471
/************************************************************************/
3472
3473
CPLErr GTIFWktFromMemBuf(int nSize, unsigned char *pabyBuffer, char **ppszWKT,
3474
                         double *padfGeoTransform, int *pnGCPCount,
3475
                         GDAL_GCP **ppasGCPList)
3476
0
{
3477
0
    OGRSpatialReferenceH hSRS = nullptr;
3478
0
    if (ppszWKT)
3479
0
        *ppszWKT = nullptr;
3480
0
    CPLErr eErr =
3481
0
        GTIFWktFromMemBufEx(nSize, pabyBuffer, &hSRS, padfGeoTransform,
3482
0
                            pnGCPCount, ppasGCPList, nullptr, nullptr);
3483
0
    if (eErr == CE_None)
3484
0
    {
3485
0
        if (hSRS && ppszWKT)
3486
0
        {
3487
0
            OSRExportToWkt(hSRS, ppszWKT);
3488
0
        }
3489
0
    }
3490
0
    OSRDestroySpatialReference(hSRS);
3491
0
    return eErr;
3492
0
}
3493
3494
CPLErr GTIFWktFromMemBufEx(int nSize, unsigned char *pabyBuffer,
3495
                           OGRSpatialReferenceH *phSRS,
3496
                           double *padfGeoTransform, int *pnGCPCount,
3497
                           GDAL_GCP **ppasGCPList, int *pbPixelIsPoint,
3498
                           char ***ppapszRPCMD)
3499
3500
0
{
3501
0
    const std::string osFilename(
3502
0
        VSIMemGenerateHiddenFilename("wkt_from_mem_buf.tif"));
3503
3504
    /* -------------------------------------------------------------------- */
3505
    /*      Initialization of libtiff and libgeotiff.                       */
3506
    /* -------------------------------------------------------------------- */
3507
0
    GTiffOneTimeInit();  // For RPC tag.
3508
0
    LibgeotiffOneTimeInit();
3509
3510
    /* -------------------------------------------------------------------- */
3511
    /*      Create a memory file from the buffer.                           */
3512
    /* -------------------------------------------------------------------- */
3513
0
    VSILFILE *fp =
3514
0
        VSIFileFromMemBuffer(osFilename.c_str(), pabyBuffer, nSize, FALSE);
3515
0
    if (fp == nullptr)
3516
0
        return CE_Failure;
3517
3518
    /* -------------------------------------------------------------------- */
3519
    /*      Initialize access to the memory geotiff structure.              */
3520
    /* -------------------------------------------------------------------- */
3521
0
    TIFF *hTIFF = VSI_TIFFOpen(osFilename.c_str(), "rc", fp);
3522
3523
0
    if (hTIFF == nullptr)
3524
0
    {
3525
0
        CPLError(CE_Failure, CPLE_AppDefined,
3526
0
                 "TIFF/GeoTIFF structure is corrupt.");
3527
0
        VSIUnlink(osFilename.c_str());
3528
0
        CPL_IGNORE_RET_VAL(VSIFCloseL(fp));
3529
0
        return CE_Failure;
3530
0
    }
3531
3532
    /* -------------------------------------------------------------------- */
3533
    /*      Get the projection definition.                                  */
3534
    /* -------------------------------------------------------------------- */
3535
0
    bool bPixelIsPoint = false;
3536
0
    bool bPointGeoIgnore = false;
3537
0
    unsigned short nRasterType = 0;
3538
3539
0
    GTIF *hGTIF = GTIFNew(hTIFF);
3540
0
    if (hGTIF)
3541
0
        GTIFAttachPROJContext(hGTIF, OSRGetProjTLSContext());
3542
3543
0
    if (hGTIF != nullptr &&
3544
0
        GDALGTIFKeyGetSHORT(hGTIF, GTRasterTypeGeoKey, &nRasterType, 0, 1) ==
3545
0
            1 &&
3546
0
        nRasterType == static_cast<unsigned short>(RasterPixelIsPoint))
3547
0
    {
3548
0
        bPixelIsPoint = true;
3549
0
        bPointGeoIgnore =
3550
0
            CPLTestBool(CPLGetConfigOption("GTIFF_POINT_GEO_IGNORE", "FALSE"));
3551
0
    }
3552
0
    if (pbPixelIsPoint)
3553
0
        *pbPixelIsPoint = bPixelIsPoint;
3554
0
    if (ppapszRPCMD)
3555
0
        *ppapszRPCMD = nullptr;
3556
3557
0
    if (phSRS)
3558
0
    {
3559
0
        *phSRS = nullptr;
3560
0
        if (hGTIF != nullptr)
3561
0
        {
3562
0
            GTIFDefn *psGTIFDefn = GTIFAllocDefn();
3563
0
            if (GTIFGetDefn(hGTIF, psGTIFDefn))
3564
0
            {
3565
0
                *phSRS = GTIFGetOGISDefnAsOSR(hGTIF, psGTIFDefn);
3566
0
            }
3567
0
            GTIFFreeDefn(psGTIFDefn);
3568
0
        }
3569
0
    }
3570
0
    if (hGTIF)
3571
0
        GTIFFree(hGTIF);
3572
3573
    /* -------------------------------------------------------------------- */
3574
    /*      Get geotransform or tiepoints.                                  */
3575
    /* -------------------------------------------------------------------- */
3576
0
    double *padfTiePoints = nullptr;
3577
0
    double *padfScale = nullptr;
3578
0
    double *padfMatrix = nullptr;
3579
0
    int16_t nCount = 0;
3580
3581
0
    padfGeoTransform[0] = 0.0;
3582
0
    padfGeoTransform[1] = 1.0;
3583
0
    padfGeoTransform[2] = 0.0;
3584
0
    padfGeoTransform[3] = 0.0;
3585
0
    padfGeoTransform[4] = 0.0;
3586
0
    padfGeoTransform[5] = 1.0;
3587
3588
0
    *pnGCPCount = 0;
3589
0
    *ppasGCPList = nullptr;
3590
3591
0
    if (TIFFGetField(hTIFF, TIFFTAG_GEOPIXELSCALE, &nCount, &padfScale) &&
3592
0
        nCount >= 2)
3593
0
    {
3594
0
        padfGeoTransform[1] = padfScale[0];
3595
0
        padfGeoTransform[5] = -std::abs(padfScale[1]);
3596
3597
0
        if (TIFFGetField(hTIFF, TIFFTAG_GEOTIEPOINTS, &nCount,
3598
0
                         &padfTiePoints) &&
3599
0
            nCount >= 6)
3600
0
        {
3601
0
            padfGeoTransform[0] =
3602
0
                padfTiePoints[3] - padfTiePoints[0] * padfGeoTransform[1];
3603
0
            padfGeoTransform[3] =
3604
0
                padfTiePoints[4] - padfTiePoints[1] * padfGeoTransform[5];
3605
3606
            // Adjust for pixel is point in transform.
3607
0
            if (bPixelIsPoint && !bPointGeoIgnore)
3608
0
            {
3609
0
                padfGeoTransform[0] -=
3610
0
                    padfGeoTransform[1] * 0.5 + padfGeoTransform[2] * 0.5;
3611
0
                padfGeoTransform[3] -=
3612
0
                    padfGeoTransform[4] * 0.5 + padfGeoTransform[5] * 0.5;
3613
0
            }
3614
0
        }
3615
0
    }
3616
0
    else if (TIFFGetField(hTIFF, TIFFTAG_GEOTIEPOINTS, &nCount,
3617
0
                          &padfTiePoints) &&
3618
0
             nCount >= 6)
3619
0
    {
3620
0
        *pnGCPCount = nCount / 6;
3621
0
        *ppasGCPList =
3622
0
            static_cast<GDAL_GCP *>(CPLCalloc(sizeof(GDAL_GCP), *pnGCPCount));
3623
3624
0
        for (int iGCP = 0; iGCP < *pnGCPCount; iGCP++)
3625
0
        {
3626
0
            char szID[32] = {};
3627
0
            GDAL_GCP *psGCP = *ppasGCPList + iGCP;
3628
3629
0
            snprintf(szID, sizeof(szID), "%d", iGCP + 1);
3630
0
            psGCP->pszId = CPLStrdup(szID);
3631
0
            psGCP->pszInfo = CPLStrdup("");
3632
0
            psGCP->dfGCPPixel = padfTiePoints[iGCP * 6 + 0];
3633
0
            psGCP->dfGCPLine = padfTiePoints[iGCP * 6 + 1];
3634
0
            psGCP->dfGCPX = padfTiePoints[iGCP * 6 + 3];
3635
0
            psGCP->dfGCPY = padfTiePoints[iGCP * 6 + 4];
3636
0
            psGCP->dfGCPZ = padfTiePoints[iGCP * 6 + 5];
3637
0
        }
3638
0
    }
3639
0
    else if (TIFFGetField(hTIFF, TIFFTAG_GEOTRANSMATRIX, &nCount,
3640
0
                          &padfMatrix) &&
3641
0
             nCount == 16)
3642
0
    {
3643
0
        padfGeoTransform[0] = padfMatrix[3];
3644
0
        padfGeoTransform[1] = padfMatrix[0];
3645
0
        padfGeoTransform[2] = padfMatrix[1];
3646
0
        padfGeoTransform[3] = padfMatrix[7];
3647
0
        padfGeoTransform[4] = padfMatrix[4];
3648
0
        padfGeoTransform[5] = padfMatrix[5];
3649
0
    }
3650
3651
    /* -------------------------------------------------------------------- */
3652
    /*      Read RPC                                                        */
3653
    /* -------------------------------------------------------------------- */
3654
0
    if (ppapszRPCMD != nullptr)
3655
0
    {
3656
0
        *ppapszRPCMD = GTiffDatasetReadRPCTag(hTIFF);
3657
0
    }
3658
3659
    /* -------------------------------------------------------------------- */
3660
    /*      Cleanup.                                                        */
3661
    /* -------------------------------------------------------------------- */
3662
0
    XTIFFClose(hTIFF);
3663
0
    CPL_IGNORE_RET_VAL(VSIFCloseL(fp));
3664
3665
0
    VSIUnlink(osFilename.c_str());
3666
3667
0
    if (phSRS && *phSRS == nullptr)
3668
0
        return CE_Failure;
3669
3670
0
    return CE_None;
3671
0
}
3672
3673
/************************************************************************/
3674
/*                         GTIFMemBufFromWkt()                          */
3675
/************************************************************************/
3676
3677
CPLErr GTIFMemBufFromWkt(const char *pszWKT, const double *padfGeoTransform,
3678
                         int nGCPCount, const GDAL_GCP *pasGCPList, int *pnSize,
3679
                         unsigned char **ppabyBuffer)
3680
0
{
3681
0
    OGRSpatialReference oSRS;
3682
0
    if (pszWKT != nullptr)
3683
0
        oSRS.importFromWkt(pszWKT);
3684
0
    return GTIFMemBufFromSRS(OGRSpatialReference::ToHandle(&oSRS),
3685
0
                             padfGeoTransform, nGCPCount, pasGCPList, pnSize,
3686
0
                             ppabyBuffer, FALSE, nullptr);
3687
0
}
3688
3689
CPLErr GTIFMemBufFromSRS(OGRSpatialReferenceH hSRS,
3690
                         const double *padfGeoTransform, int nGCPCount,
3691
                         const GDAL_GCP *pasGCPList, int *pnSize,
3692
                         unsigned char **ppabyBuffer, int bPixelIsPoint,
3693
                         char **papszRPCMD)
3694
3695
0
{
3696
0
    const std::string osFilename(
3697
0
        VSIMemGenerateHiddenFilename("wkt_from_mem_buf.tif"));
3698
3699
    /* -------------------------------------------------------------------- */
3700
    /*      Initialization of libtiff and libgeotiff.                       */
3701
    /* -------------------------------------------------------------------- */
3702
0
    GTiffOneTimeInit();  // For RPC tag.
3703
0
    LibgeotiffOneTimeInit();
3704
3705
    /* -------------------------------------------------------------------- */
3706
    /*      Initialize access to the memory geotiff structure.              */
3707
    /* -------------------------------------------------------------------- */
3708
0
    VSILFILE *fpL = VSIFOpenL(osFilename.c_str(), "w");
3709
0
    if (fpL == nullptr)
3710
0
        return CE_Failure;
3711
3712
0
    TIFF *hTIFF = VSI_TIFFOpen(osFilename.c_str(), "w", fpL);
3713
3714
0
    if (hTIFF == nullptr)
3715
0
    {
3716
0
        CPLError(CE_Failure, CPLE_AppDefined,
3717
0
                 "TIFF/GeoTIFF structure is corrupt.");
3718
0
        CPL_IGNORE_RET_VAL(VSIFCloseL(fpL));
3719
0
        VSIUnlink(osFilename.c_str());
3720
0
        return CE_Failure;
3721
0
    }
3722
3723
    /* -------------------------------------------------------------------- */
3724
    /*      Write some minimal set of image parameters.                     */
3725
    /* -------------------------------------------------------------------- */
3726
0
    TIFFSetField(hTIFF, TIFFTAG_IMAGEWIDTH, 1);
3727
0
    TIFFSetField(hTIFF, TIFFTAG_IMAGELENGTH, 1);
3728
0
    TIFFSetField(hTIFF, TIFFTAG_BITSPERSAMPLE, 8);
3729
0
    TIFFSetField(hTIFF, TIFFTAG_SAMPLESPERPIXEL, 1);
3730
0
    TIFFSetField(hTIFF, TIFFTAG_ROWSPERSTRIP, 1);
3731
0
    TIFFSetField(hTIFF, TIFFTAG_PLANARCONFIG, PLANARCONFIG_CONTIG);
3732
0
    TIFFSetField(hTIFF, TIFFTAG_PHOTOMETRIC, PHOTOMETRIC_MINISBLACK);
3733
3734
    /* -------------------------------------------------------------------- */
3735
    /*      Get the projection definition.                                  */
3736
    /* -------------------------------------------------------------------- */
3737
3738
0
    bool bPointGeoIgnore = false;
3739
0
    if (bPixelIsPoint)
3740
0
    {
3741
0
        bPointGeoIgnore =
3742
0
            CPLTestBool(CPLGetConfigOption("GTIFF_POINT_GEO_IGNORE", "FALSE"));
3743
0
    }
3744
3745
0
    GTIF *hGTIF = nullptr;
3746
0
    if (hSRS != nullptr || bPixelIsPoint)
3747
0
    {
3748
0
        hGTIF = GTIFNew(hTIFF);
3749
0
        if (hGTIF)
3750
0
        {
3751
0
            GTIFAttachPROJContext(hGTIF, OSRGetProjTLSContext());
3752
3753
0
            if (hSRS != nullptr)
3754
0
                GTIFSetFromOGISDefnEx(hGTIF, hSRS, GEOTIFF_KEYS_STANDARD,
3755
0
                                      GEOTIFF_VERSION_1_0);
3756
3757
0
            if (bPixelIsPoint)
3758
0
            {
3759
0
                GTIFKeySet(hGTIF, GTRasterTypeGeoKey, TYPE_SHORT, 1,
3760
0
                           RasterPixelIsPoint);
3761
0
            }
3762
3763
0
            GTIFWriteKeys(hGTIF);
3764
0
            GTIFFree(hGTIF);
3765
0
        }
3766
0
    }
3767
3768
    /* -------------------------------------------------------------------- */
3769
    /*      Set the geotransform, or GCPs.                                  */
3770
    /* -------------------------------------------------------------------- */
3771
3772
0
    if (padfGeoTransform[0] != 0.0 || padfGeoTransform[1] != 1.0 ||
3773
0
        padfGeoTransform[2] != 0.0 || padfGeoTransform[3] != 0.0 ||
3774
0
        padfGeoTransform[4] != 0.0 || std::abs(padfGeoTransform[5]) != 1.0)
3775
0
    {
3776
3777
0
        if (padfGeoTransform[2] == 0.0 && padfGeoTransform[4] == 0.0)
3778
0
        {
3779
0
            double adfPixelScale[3] = {padfGeoTransform[1],
3780
0
                                       fabs(padfGeoTransform[5]), 0.0};
3781
3782
0
            TIFFSetField(hTIFF, TIFFTAG_GEOPIXELSCALE, 3, adfPixelScale);
3783
3784
0
            double adfTiePoints[6] = {
3785
0
                0.0, 0.0, 0.0, padfGeoTransform[0], padfGeoTransform[3], 0.0};
3786
3787
0
            if (bPixelIsPoint && !bPointGeoIgnore)
3788
0
            {
3789
0
                adfTiePoints[3] +=
3790
0
                    padfGeoTransform[1] * 0.5 + padfGeoTransform[2] * 0.5;
3791
0
                adfTiePoints[4] +=
3792
0
                    padfGeoTransform[4] * 0.5 + padfGeoTransform[5] * 0.5;
3793
0
            }
3794
3795
0
            TIFFSetField(hTIFF, TIFFTAG_GEOTIEPOINTS, 6, adfTiePoints);
3796
0
        }
3797
0
        else
3798
0
        {
3799
0
            double adfMatrix[16] = {0.0};
3800
3801
0
            adfMatrix[0] = padfGeoTransform[1];
3802
0
            adfMatrix[1] = padfGeoTransform[2];
3803
0
            adfMatrix[3] = padfGeoTransform[0];
3804
0
            adfMatrix[4] = padfGeoTransform[4];
3805
0
            adfMatrix[5] = padfGeoTransform[5];
3806
0
            adfMatrix[7] = padfGeoTransform[3];
3807
0
            adfMatrix[15] = 1.0;
3808
3809
0
            if (bPixelIsPoint && !bPointGeoIgnore)
3810
0
            {
3811
0
                adfMatrix[3] +=
3812
0
                    padfGeoTransform[1] * 0.5 + padfGeoTransform[2] * 0.5;
3813
0
                adfMatrix[7] +=
3814
0
                    padfGeoTransform[4] * 0.5 + padfGeoTransform[5] * 0.5;
3815
0
            }
3816
3817
0
            TIFFSetField(hTIFF, TIFFTAG_GEOTRANSMATRIX, 16, adfMatrix);
3818
0
        }
3819
0
    }
3820
3821
    /* -------------------------------------------------------------------- */
3822
    /*      Otherwise write tiepoints if they are available.                */
3823
    /* -------------------------------------------------------------------- */
3824
0
    else if (nGCPCount > 0)
3825
0
    {
3826
0
        double *padfTiePoints =
3827
0
            static_cast<double *>(CPLMalloc(6 * sizeof(double) * nGCPCount));
3828
3829
0
        for (int iGCP = 0; iGCP < nGCPCount; iGCP++)
3830
0
        {
3831
3832
0
            padfTiePoints[iGCP * 6 + 0] = pasGCPList[iGCP].dfGCPPixel;
3833
0
            padfTiePoints[iGCP * 6 + 1] = pasGCPList[iGCP].dfGCPLine;
3834
0
            padfTiePoints[iGCP * 6 + 2] = 0;
3835
0
            padfTiePoints[iGCP * 6 + 3] = pasGCPList[iGCP].dfGCPX;
3836
0
            padfTiePoints[iGCP * 6 + 4] = pasGCPList[iGCP].dfGCPY;
3837
0
            padfTiePoints[iGCP * 6 + 5] = pasGCPList[iGCP].dfGCPZ;
3838
0
        }
3839
3840
0
        TIFFSetField(hTIFF, TIFFTAG_GEOTIEPOINTS, 6 * nGCPCount, padfTiePoints);
3841
0
        CPLFree(padfTiePoints);
3842
0
    }
3843
3844
    /* -------------------------------------------------------------------- */
3845
    /*      Write RPC                                                       */
3846
    /* -------------------------------------------------------------------- */
3847
0
    if (papszRPCMD != nullptr)
3848
0
    {
3849
0
        GTiffDatasetWriteRPCTag(hTIFF, papszRPCMD);
3850
0
    }
3851
3852
    /* -------------------------------------------------------------------- */
3853
    /*      Cleanup and return the created memory buffer.                   */
3854
    /* -------------------------------------------------------------------- */
3855
0
    GByte bySmallImage = 0;
3856
3857
0
    TIFFWriteEncodedStrip(hTIFF, 0, reinterpret_cast<char *>(&bySmallImage), 1);
3858
0
    TIFFWriteCheck(hTIFF, TIFFIsTiled(hTIFF), "GTIFMemBufFromWkt");
3859
0
    TIFFWriteDirectory(hTIFF);
3860
3861
0
    XTIFFClose(hTIFF);
3862
0
    CPL_IGNORE_RET_VAL(VSIFCloseL(fpL));
3863
3864
    /* -------------------------------------------------------------------- */
3865
    /*      Read back from the memory buffer.  It would be preferable       */
3866
    /*      to be able to "steal" the memory buffer, but there isn't        */
3867
    /*      currently any support for this.                                 */
3868
    /* -------------------------------------------------------------------- */
3869
0
    GUIntBig nBigLength = 0;
3870
3871
0
    *ppabyBuffer = VSIGetMemFileBuffer(osFilename.c_str(), &nBigLength, TRUE);
3872
0
    *pnSize = static_cast<int>(nBigLength);
3873
3874
0
    return CE_None;
3875
0
}