• Home
  • Features
  • Pricing
  • Docs
  • Announcements
  • Sign In

OSGeo / gdal / 15797738239

21 Jun 2025 04:51PM UTC coverage: 71.076% (+0.001%) from 71.075%
15797738239

Pull #12622

github

web-flow
Merge 19559fd15 into 645a00347
Pull Request #12622: VSI Win32: implement OpenDir()

574181 of 807840 relevant lines covered (71.08%)

249831.39 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

90.58
/frmts/gtiff/gt_wkt_srs.cpp
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
#define CT_HotineObliqueMercatorAzimuthCenter 9815
54
#endif
55

56
#if !defined(GTIFAtof)
57
#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
#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()
9,775✔
105
{
106
    std::lock_guard<std::mutex> oLock(oDeleteMutex);
9,775✔
107

108
    static bool bOneTimeInitDone = false;
109

110
    if (bOneTimeInitDone)
9,775✔
111
        return;
9,166✔
112

113
    bOneTimeInitDone = true;
609✔
114

115
    // This isn't thread-safe, so better do it now
116
    XTIFFInitialize();
609✔
117
}
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)
39,582✔
127

128
{
129
    if (*ppszTarget == nullptr)
39,582✔
130
        return;
54✔
131

132
    char *pszTempString = CPLStrdup(*ppszTarget);
39,528✔
133
    GTIFFreeMemory(*ppszTarget);
39,528✔
134
    *ppszTarget = pszTempString;
39,528✔
135
}
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)
13✔
145

146
{
147
    char *pszDatum = *ppszDatum;
13✔
148
    if (!pszDatum || pszDatum[0] == '\0')
13✔
149
        return;
×
150

151
    /* -------------------------------------------------------------------- */
152
    /*      Translate non-alphanumeric values to underscores.               */
153
    /* -------------------------------------------------------------------- */
154
    for (int i = 0; pszDatum[i] != '\0'; i++)
435✔
155
    {
156
        if (pszDatum[i] != '+' && !(pszDatum[i] >= 'A' && pszDatum[i] <= 'Z') &&
422✔
157
            !(pszDatum[i] >= 'a' && pszDatum[i] <= 'z') &&
375✔
158
            !(pszDatum[i] >= '0' && pszDatum[i] <= '9'))
124✔
159
        {
160
            pszDatum[i] = '_';
62✔
161
        }
162
    }
163

164
    /* -------------------------------------------------------------------- */
165
    /*      Remove repeated and trailing underscores.                       */
166
    /* -------------------------------------------------------------------- */
167
    int j = 0;  // Used after for.
13✔
168
    for (int i = 1; pszDatum[i] != '\0'; i++)
422✔
169
    {
170
        if (pszDatum[j] == '_' && pszDatum[i] == '_')
409✔
171
            continue;
8✔
172

173
        pszDatum[++j] = pszDatum[i];
401✔
174
    }
175
    if (pszDatum[j] == '_')
13✔
176
        pszDatum[j] = '\0';
8✔
177
    else
178
        pszDatum[j + 1] = '\0';
5✔
179

180
    /* -------------------------------------------------------------------- */
181
    /*      Search for datum equivalences.  Specific massaged names get     */
182
    /*      mapped to OpenGIS specified names.                              */
183
    /* -------------------------------------------------------------------- */
184
    for (int i = 0; papszDatumEquiv[i] != nullptr; i += 2)
77✔
185
    {
186
        if (EQUAL(*ppszDatum, papszDatumEquiv[i]))
68✔
187
        {
188
            CPLFree(*ppszDatum);
4✔
189
            *ppszDatum = CPLStrdup(papszDatumEquiv[i + 1]);
4✔
190
            return;
4✔
191
        }
192
    }
193
}
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)
335✔
223

224
{
225
    if (strstr(pszCitation, "IMAGINE GeoTIFF") == nullptr)
335✔
226
        return;
331✔
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
    char *pszSkip = pszCitation + strlen(pszCitation) - 1;
4✔
233

234
    for (; pszSkip != pszCitation && *pszSkip != '$'; pszSkip--)
428✔
235
    {
236
    }
237

238
    if (*pszSkip == '$')
4✔
239
        pszSkip++;
2✔
240
    if (*pszSkip == '\n')
4✔
241
        pszSkip++;
2✔
242

243
    memmove(pszCitation, pszSkip, strlen(pszSkip) + 1);
4✔
244

245
    /* -------------------------------------------------------------------- */
246
    /*      Convert any newlines into spaces, they really gum up the        */
247
    /*      WKT.                                                            */
248
    /* -------------------------------------------------------------------- */
249
    for (int i = 0; pszCitation[i] != '\0'; i++)
428✔
250
    {
251
        if (pszCitation[i] == '\n')
424✔
252
            pszCitation[i] = ' ';
14✔
253
    }
254
}
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,
15✔
314
                                            OGRSpatialReference &oSRS,
315
                                            const char *pszVertCSName,
316
                                            int verticalDatum,
317
                                            int verticalUnits)
318
{
319
    /* -------------------------------------------------------------------- */
320
    /*      Setup VERT_CS with citation if present.                         */
321
    /* -------------------------------------------------------------------- */
322
    oSRS.SetNode("COMPD_CS|VERT_CS", pszVertCSName);
15✔
323

324
    /* -------------------------------------------------------------------- */
325
    /*      Setup the vertical datum.                                       */
326
    /* -------------------------------------------------------------------- */
327
    std::string osVDatumName = "unknown";
30✔
328
    const char *pszVDatumType = "2005";  // CS_VD_GeoidModelDerived
15✔
329
    std::string osVDatumAuthName;
30✔
330
    int nVDatumCode = 0;
15✔
331

332
    if (verticalDatum > 0 && verticalDatum != KvUserDefined)
15✔
333
    {
334
        osVDatumAuthName = "EPSG";
3✔
335
        nVDatumCode = verticalDatum;
3✔
336

337
        char szCode[12];
338
        snprintf(szCode, sizeof(szCode), "%d", verticalDatum);
3✔
339
        auto ctx =
340
            static_cast<PJ_CONTEXT *>(GTIFGetPROJContext(hGTIF, true, nullptr));
3✔
341
        auto datum = proj_create_from_database(ctx, "EPSG", szCode,
3✔
342
                                               PJ_CATEGORY_DATUM, 0, nullptr);
343
        if (datum)
3✔
344
        {
345
            const char *pszName = proj_get_name(datum);
3✔
346
            if (pszName)
3✔
347
            {
348
                osVDatumName = pszName;
3✔
349
            }
350
            proj_destroy(datum);
3✔
351
        }
3✔
352
    }
353
    else if (verticalDatum == KvUserDefined)
12✔
354
    {
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
        auto ctx =
358
            static_cast<PJ_CONTEXT *>(GTIFGetPROJContext(hGTIF, true, nullptr));
3✔
359
        const auto type = PJ_TYPE_VERTICAL_CRS;
3✔
360
        auto list = proj_create_from_name(ctx, nullptr, pszVertCSName, &type, 1,
3✔
361
                                          true,  // exact match
362
                                          -1,    // result set limit size,
363
                                          nullptr);
364
        if (list)
3✔
365
        {
366
            // If we have several matches, check they all refer to the
367
            // same datum
368
            bool bGoOn = true;
3✔
369
            int ncount = proj_list_get_count(list);
3✔
370
            for (int i = 0; bGoOn && i < ncount; ++i)
6✔
371
            {
372
                auto crs = proj_list_get(ctx, list, i);
3✔
373
                if (crs)
3✔
374
                {
375
                    auto datum = proj_crs_get_datum(ctx, crs);
3✔
376
                    if (datum)
3✔
377
                    {
378
                        osVDatumName = proj_get_name(datum);
3✔
379
                        const char *pszAuthName =
380
                            proj_get_id_auth_name(datum, 0);
3✔
381
                        const char *pszCode = proj_get_id_code(datum, 0);
3✔
382
                        if (pszCode && atoi(pszCode) && pszAuthName)
3✔
383
                        {
384
                            if (osVDatumAuthName.empty())
3✔
385
                            {
386
                                osVDatumAuthName = pszAuthName;
1✔
387
                                nVDatumCode = atoi(pszCode);
1✔
388
                            }
389
                            else if (osVDatumAuthName != pszAuthName ||
4✔
390
                                     nVDatumCode != atoi(pszCode))
2✔
391
                            {
392
                                osVDatumAuthName.clear();
×
393
                                nVDatumCode = 0;
×
394
                                bGoOn = false;
×
395
                            }
396
                        }
397
                        proj_destroy(datum);
3✔
398
                    }
399
                    proj_destroy(crs);
3✔
400
                }
401
            }
402
        }
403
        proj_list_destroy(list);
3✔
404
    }
405

406
    oSRS.SetNode("COMPD_CS|VERT_CS|VERT_DATUM", osVDatumName.c_str());
15✔
407
    oSRS.GetAttrNode("COMPD_CS|VERT_CS|VERT_DATUM")
408
        ->AddChild(new OGR_SRSNode(pszVDatumType));
15✔
409
    if (!osVDatumAuthName.empty())
15✔
410
        oSRS.SetAuthority("COMPD_CS|VERT_CS|VERT_DATUM",
4✔
411
                          osVDatumAuthName.c_str(), nVDatumCode);
412

413
    /* -------------------------------------------------------------------- */
414
    /*      Set the vertical units.                                         */
415
    /* -------------------------------------------------------------------- */
416
    if (verticalUnits > 0 && verticalUnits != KvUserDefined &&
15✔
417
        verticalUnits != 9001)
418
    {
419
        char szCode[12];
420
        snprintf(szCode, sizeof(szCode), "%d", verticalUnits);
1✔
421
        auto ctx =
422
            static_cast<PJ_CONTEXT *>(GTIFGetPROJContext(hGTIF, true, nullptr));
1✔
423
        const char *pszName = nullptr;
1✔
424
        double dfInMeters = 0.0;
1✔
425
        if (proj_uom_get_info_from_database(ctx, "EPSG", szCode, &pszName,
1✔
426
                                            &dfInMeters, nullptr))
1✔
427
        {
428
            if (pszName)
1✔
429
                oSRS.SetNode("COMPD_CS|VERT_CS|UNIT", pszName);
1✔
430

431
            char szInMeters[128] = {};
1✔
432
            CPLsnprintf(szInMeters, sizeof(szInMeters), "%.16g", dfInMeters);
1✔
433
            oSRS.GetAttrNode("COMPD_CS|VERT_CS|UNIT")
434
                ->AddChild(new OGR_SRSNode(szInMeters));
1✔
435
        }
436

437
        oSRS.SetAuthority("COMPD_CS|VERT_CS|UNIT", "EPSG", verticalUnits);
1✔
438
    }
439
    else
440
    {
441
        oSRS.SetNode("COMPD_CS|VERT_CS|UNIT", "metre");
14✔
442
        oSRS.GetAttrNode("COMPD_CS|VERT_CS|UNIT")
443
            ->AddChild(new OGR_SRSNode("1.0"));
14✔
444
        oSRS.SetAuthority("COMPD_CS|VERT_CS|UNIT", "EPSG", 9001);
14✔
445
    }
446

447
    /* -------------------------------------------------------------------- */
448
    /*      Set the axis and VERT_CS authority.                             */
449
    /* -------------------------------------------------------------------- */
450
    oSRS.SetNode("COMPD_CS|VERT_CS|AXIS", "Up");
15✔
451
    oSRS.GetAttrNode("COMPD_CS|VERT_CS|AXIS")->AddChild(new OGR_SRSNode("UP"));
15✔
452
}
15✔
453

454
/************************************************************************/
455
/*                    GTIFGetEPSGOfficialName()                         */
456
/************************************************************************/
457

458
static char *GTIFGetEPSGOfficialName(GTIF *hGTIF, PJ_TYPE searchType,
47✔
459
                                     const char *pszName)
460
{
461
    char *pszRet = nullptr;
47✔
462
    /* Search in database the corresponding EPSG 'official' name */
463
    auto ctx =
464
        static_cast<PJ_CONTEXT *>(GTIFGetPROJContext(hGTIF, true, nullptr));
47✔
465
    auto list = proj_create_from_name(ctx, "EPSG", pszName, &searchType, 1,
47✔
466
                                      false, /* approximate match */
467
                                      1, nullptr);
468
    if (list)
47✔
469
    {
470
        const auto listSize = proj_list_get_count(list);
47✔
471
        if (listSize == 1)
47✔
472
        {
473
            auto obj = proj_list_get(ctx, list, 0);
4✔
474
            if (obj)
4✔
475
            {
476
                const char *pszOfficialName = proj_get_name(obj);
4✔
477
                if (pszOfficialName)
4✔
478
                {
479
                    pszRet = CPLStrdup(pszOfficialName);
4✔
480
                }
481
            }
482
            proj_destroy(obj);
4✔
483
        }
484
        proj_list_destroy(list);
47✔
485
    }
486
    return pszRet;
47✔
487
}
488

489
/************************************************************************/
490
/*                      GTIFGetOGISDefnAsOSR()                          */
491
/************************************************************************/
492

493
OGRSpatialReferenceH GTIFGetOGISDefnAsOSR(GTIF *hGTIF, GTIFDefn *psDefn)
8,325✔
494

495
{
496
    OGRSpatialReference oSRS;
16,649✔
497

498
    LibgeotiffOneTimeInit();
8,325✔
499

500
#if LIBGEOTIFF_VERSION >= 1600
501
    void *projContext = GTIFGetPROJContext(hGTIF, FALSE, nullptr);
8,325✔
502
#endif
503

504
    /* -------------------------------------------------------------------- */
505
    /*  Handle non-standard coordinate systems where GTModelTypeGeoKey      */
506
    /*  is not defined, but ProjectedCSTypeGeoKey is defined (ticket #3019) */
507
    /* -------------------------------------------------------------------- */
508
    if (psDefn->Model == KvUserDefined && psDefn->PCS != KvUserDefined)
8,325✔
509
    {
510
        psDefn->Model = ModelTypeProjected;
5✔
511
    }
512

513
    /* ==================================================================== */
514
    /*      Read keys related to vertical component.                        */
515
    /* ==================================================================== */
516
    unsigned short verticalCSType = 0;
8,325✔
517
    unsigned short verticalDatum = 0;
8,325✔
518
    unsigned short verticalUnits = 0;
8,325✔
519

520
    GDALGTIFKeyGetSHORT(hGTIF, VerticalCSTypeGeoKey, &verticalCSType, 0, 1);
8,325✔
521
    GDALGTIFKeyGetSHORT(hGTIF, VerticalDatumGeoKey, &verticalDatum, 0, 1);
8,325✔
522
    GDALGTIFKeyGetSHORT(hGTIF, VerticalUnitsGeoKey, &verticalUnits, 0, 1);
8,325✔
523

524
    if (verticalCSType != 0 || verticalDatum != 0 || verticalUnits != 0)
8,325✔
525
    {
526
        int versions[3];
527
        GTIFDirectoryInfo(hGTIF, versions, nullptr);
42✔
528
        // GeoTIFF 1.0
529
        if (versions[0] == 1 && versions[1] == 1 && versions[2] == 0)
42✔
530
        {
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
            if ((verticalCSType >= 5101 && verticalCSType <= 5112) &&
12✔
544
                verticalDatum == 0)
×
545
            {
546
                verticalDatum = verticalCSType;
×
547
                verticalCSType = verticalDatum + 600;
×
548
            }
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
            if (verticalCSType == 5171)
12✔
558
            {
559
                verticalDatum = 5171;
×
560
                verticalCSType = 5773;
×
561
            }
562
        }
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
        if ((verticalCSType >= 5001 && verticalCSType <= 5033) &&
42✔
574
            verticalDatum == 0)
1✔
575
        {
576
            verticalDatum = verticalCSType + 1000;
1✔
577
            verticalCSType = 0;
1✔
578
        }
579
    }
580

581
    /* -------------------------------------------------------------------- */
582
    /*      Handle non-standard coordinate systems as LOCAL_CS.             */
583
    /* -------------------------------------------------------------------- */
584
    if (psDefn->Model != ModelTypeProjected &&
8,325✔
585
        psDefn->Model != ModelTypeGeographic &&
4,609✔
586
        psDefn->Model != ModelTypeGeocentric)
346✔
587
    {
588
        char szPeStr[2400] = {'\0'};
344✔
589

590
        /** check if there is a pe string citation key **/
591
        if (GDALGTIFKeyGetASCII(hGTIF, PCSCitationGeoKey, szPeStr,
344✔
592
                                sizeof(szPeStr)) &&
353✔
593
            strstr(szPeStr, "ESRI PE String = "))
9✔
594
        {
595
            const char *pszWKT = szPeStr + strlen("ESRI PE String = ");
9✔
596
            oSRS.importFromWkt(pszWKT);
9✔
597

598
            if (strstr(pszWKT,
9✔
599
                       "PROJCS[\"WGS_1984_Web_Mercator_Auxiliary_Sphere\""))
600
            {
601
                oSRS.SetExtension(
6✔
602
                    "PROJCS", "PROJ4",
603
                    "+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 "
604
                    "+x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null "
605
                    "+wktext  +no_defs");  // TODO(schwehr): Why 2 spaces?
606
            }
607

608
            return OGRSpatialReference::ToHandle(oSRS.Clone());
9✔
609
        }
610
        else
611
        {
612
            char *pszUnitsName = nullptr;
335✔
613
            char szPCSName[300] = {'\0'};
335✔
614
            int nKeyCount = 0;
335✔
615
            int anVersion[3] = {0};
335✔
616

617
            GTIFDirectoryInfo(hGTIF, anVersion, &nKeyCount);
335✔
618

619
            if (nKeyCount > 0)  // Use LOCAL_CS if we have any geokeys at all.
335✔
620
            {
621
                // Handle citation.
622
                strcpy(szPCSName, "unnamed");
335✔
623
                if (!GDALGTIFKeyGetASCII(hGTIF, GTCitationGeoKey, szPCSName,
335✔
624
                                         sizeof(szPCSName)))
625
                    GDALGTIFKeyGetASCII(hGTIF, GeogCitationGeoKey, szPCSName,
320✔
626
                                        sizeof(szPCSName));
627

628
                GTIFCleanupImagineNames(szPCSName);
335✔
629
                oSRS.SetLocalCS(szPCSName);
335✔
630

631
                // Handle units
632
                if (psDefn->UOMLength != KvUserDefined)
335✔
633
                {
634
#if LIBGEOTIFF_VERSION >= 1600
635
                    GTIFGetUOMLengthInfoEx(projContext,
323✔
636
#else
637
                    GTIFGetUOMLengthInfo(
638
#endif
639
                                           psDefn->UOMLength, &pszUnitsName,
323✔
640
                                           nullptr);
641
                }
642

643
                if (pszUnitsName != nullptr)
335✔
644
                {
645
                    char szUOMLength[12];
646
                    snprintf(szUOMLength, sizeof(szUOMLength), "%d",
323✔
647
                             psDefn->UOMLength);
323✔
648
                    oSRS.SetTargetLinearUnits(nullptr, pszUnitsName,
323✔
649
                                              psDefn->UOMLengthInMeters, "EPSG",
650
                                              szUOMLength);
651
                }
652
                else
653
                    oSRS.SetLinearUnits("unknown", psDefn->UOMLengthInMeters);
12✔
654

655
                if (verticalUnits != 0)
335✔
656
                {
657
                    char szVertCSCitation[2048] = {0};
3✔
658
                    if (GDALGTIFKeyGetASCII(hGTIF, VerticalCitationGeoKey,
3✔
659
                                            szVertCSCitation,
660
                                            sizeof(szVertCSCitation)))
3✔
661
                    {
662
                        if (STARTS_WITH_CI(szVertCSCitation, "VCS Name = "))
1✔
663
                        {
664
                            memmove(szVertCSCitation,
×
665
                                    szVertCSCitation + strlen("VCS Name = "),
666
                                    strlen(szVertCSCitation +
×
667
                                           strlen("VCS Name = ")) +
668
                                        1);
669
                            char *pszPipeChar = strchr(szVertCSCitation, '|');
×
670
                            if (pszPipeChar)
×
671
                                *pszPipeChar = '\0';
×
672
                        }
673
                    }
674
                    else
675
                    {
676
                        strcpy(szVertCSCitation, "unknown");
2✔
677
                    }
678

679
                    const char *pszHorizontalName = oSRS.GetName();
3✔
680
                    const std::string osHorizontalName(
681
                        pszHorizontalName ? pszHorizontalName : "unnamed");
6✔
682
                    /* --------------------------------------------------------------------
683
                     */
684
                    /*      Promote to being a compound coordinate system. */
685
                    /* --------------------------------------------------------------------
686
                     */
687
                    OGR_SRSNode *poOldRoot = oSRS.GetRoot()->Clone();
3✔
688

689
                    oSRS.Clear();
3✔
690

691
                    /* --------------------------------------------------------------------
692
                     */
693
                    /*      Set COMPD_CS name. */
694
                    /* --------------------------------------------------------------------
695
                     */
696
                    char szCTString[512];
697
                    szCTString[0] = '\0';
3✔
698
                    if (GDALGTIFKeyGetASCII(hGTIF, GTCitationGeoKey, szCTString,
3✔
699
                                            sizeof(szCTString)) &&
4✔
700
                        strstr(szCTString, " = ") == nullptr)
1✔
701
                    {
702
                        oSRS.SetNode("COMPD_CS", szCTString);
1✔
703
                    }
704
                    else
705
                    {
706
                        oSRS.SetNode("COMPD_CS", (osHorizontalName + " + " +
2✔
707
                                                  szVertCSCitation)
708
                                                     .c_str());
709
                    }
710

711
                    oSRS.GetRoot()->AddChild(poOldRoot);
3✔
712

713
                    FillCompoundCRSWithManualVertCS(
3✔
714
                        hGTIF, oSRS, szVertCSCitation, verticalDatum,
715
                        verticalUnits);
716
                }
717

718
                GTIFFreeMemory(pszUnitsName);
335✔
719
            }
720
            return OGRSpatialReference::ToHandle(oSRS.Clone());
335✔
721
        }
722
    }
723

724
    /* -------------------------------------------------------------------- */
725
    /*      Handle Geocentric coordinate systems.                           */
726
    /* -------------------------------------------------------------------- */
727
    if (psDefn->Model == ModelTypeGeocentric)
7,981✔
728
    {
729
        char szName[300] = {'\0'};
2✔
730

731
        strcpy(szName, "unnamed");
2✔
732
        if (!GDALGTIFKeyGetASCII(hGTIF, GTCitationGeoKey, szName,
2✔
733
                                 sizeof(szName)))
734
            GDALGTIFKeyGetASCII(hGTIF, GeogCitationGeoKey, szName,
×
735
                                sizeof(szName));
736

737
        oSRS.SetGeocCS(szName);
2✔
738

739
        char *pszUnitsName = nullptr;
2✔
740

741
        if (psDefn->UOMLength != KvUserDefined)
2✔
742
        {
743
#if LIBGEOTIFF_VERSION >= 1600
744
            GTIFGetUOMLengthInfoEx(projContext,
2✔
745
#else
746
            GTIFGetUOMLengthInfo(
747
#endif
748
                                   psDefn->UOMLength, &pszUnitsName, nullptr);
2✔
749
        }
750

751
        if (pszUnitsName != nullptr)
2✔
752
        {
753
            char szUOMLength[12];
754
            snprintf(szUOMLength, sizeof(szUOMLength), "%d", psDefn->UOMLength);
2✔
755
            oSRS.SetTargetLinearUnits(nullptr, pszUnitsName,
2✔
756
                                      psDefn->UOMLengthInMeters, "EPSG",
757
                                      szUOMLength);
758
        }
759
        else
760
            oSRS.SetLinearUnits("unknown", psDefn->UOMLengthInMeters);
×
761

762
        GTIFFreeMemory(pszUnitsName);
2✔
763
    }
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
    const char *pszLinearUnits =
774
        CPLGetConfigOption("GTIFF_LINEAR_UNITS", "DEFAULT");
7,981✔
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
    unsigned short bLinearUnitsMarkedCorrect = FALSE;
7,981✔
785

786
    GDALGTIFKeyGetSHORT(hGTIF, ProjLinearUnitsInterpCorrectGeoKey,
7,981✔
787
                        &bLinearUnitsMarkedCorrect, 0, 1);
788

789
    if (EQUAL(pszLinearUnits, "BROKEN") &&
7,981✔
790
        psDefn->Projection == KvUserDefined && !bLinearUnitsMarkedCorrect)
3✔
791
    {
792
        for (int iParam = 0; iParam < psDefn->nParms; iParam++)
16✔
793
        {
794
            switch (psDefn->ProjParmId[iParam])
14✔
795
            {
796
                case ProjFalseEastingGeoKey:
4✔
797
                case ProjFalseNorthingGeoKey:
798
                case ProjFalseOriginEastingGeoKey:
799
                case ProjFalseOriginNorthingGeoKey:
800
                case ProjCenterEastingGeoKey:
801
                case ProjCenterNorthingGeoKey:
802
                    if (psDefn->UOMLengthInMeters != 0 &&
4✔
803
                        psDefn->UOMLengthInMeters != 1.0)
4✔
804
                    {
805
                        psDefn->ProjParm[iParam] /= psDefn->UOMLengthInMeters;
4✔
806
                        CPLDebug(
4✔
807
                            "GTIFF",
808
                            "Converting geokey to accommodate old broken file "
809
                            "due to GTIFF_LINEAR_UNITS=BROKEN setting.");
810
                    }
811
                    break;
4✔
812

813
                default:
10✔
814
                    break;
10✔
815
            }
816
        }
817
    }
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
    OGRBoolean linearUnitIsSet = FALSE;
7,981✔
824
    if (psDefn->Model == ModelTypeProjected)
7,981✔
825
    {
826
        char szCTString[512] = {'\0'};
3,716✔
827
        if (psDefn->PCS != KvUserDefined)
3,716✔
828
        {
829
            char *pszPCSName = nullptr;
3,619✔
830

831
#if LIBGEOTIFF_VERSION >= 1600
832
            GTIFGetPCSInfoEx(projContext,
3,619✔
833
#else
834
            GTIFGetPCSInfo(
835
#endif
836
                             psDefn->PCS, &pszPCSName, nullptr, nullptr,
3,619✔
837
                             nullptr);
838

839
            oSRS.SetProjCS(pszPCSName ? pszPCSName : "unnamed");
3,619✔
840
            if (pszPCSName)
3,619✔
841
                GTIFFreeMemory(pszPCSName);
3,619✔
842

843
            oSRS.SetLinearUnits("unknown", 1.0);
3,619✔
844
        }
845
        else
846
        {
847
            bool bTryGTCitationGeoKey = true;
97✔
848
            if (GDALGTIFKeyGetASCII(hGTIF, PCSCitationGeoKey, szCTString,
97✔
849
                                    sizeof(szCTString)))
97✔
850
            {
851
                bTryGTCitationGeoKey = false;
6✔
852
                if (!SetCitationToSRS(hGTIF, szCTString, sizeof(szCTString),
6✔
853
                                      PCSCitationGeoKey, &oSRS,
854
                                      &linearUnitIsSet))
855
                {
856
                    if (!STARTS_WITH_CI(szCTString, "LUnits = "))
6✔
857
                    {
858
                        oSRS.SetProjCS(szCTString);
4✔
859
                        oSRS.SetLinearUnits("unknown", 1.0);
4✔
860
                    }
861
                    else
862
                    {
863
                        bTryGTCitationGeoKey = true;
2✔
864
                    }
865
                }
866
            }
867

868
            if (bTryGTCitationGeoKey)
97✔
869
            {
870
                if (GDALGTIFKeyGetASCII(hGTIF, GTCitationGeoKey, szCTString,
93✔
871
                                        sizeof(szCTString)) &&
184✔
872
                    !SetCitationToSRS(hGTIF, szCTString, sizeof(szCTString),
91✔
873
                                      GTCitationGeoKey, &oSRS,
874
                                      &linearUnitIsSet))
875
                {
876
                    oSRS.SetNode("PROJCS", szCTString);
×
877
                    oSRS.SetLinearUnits("unknown", 1.0);
×
878
                }
879
                else
880
                {
881
                    oSRS.SetNode("PROJCS", "unnamed");
93✔
882
                    oSRS.SetLinearUnits("unknown", 1.0);
93✔
883
                }
884
            }
885
        }
886

887
        /* Handle ESRI/Erdas style state plane and UTM in citation key */
888
        if (CheckCitationKeyForStatePlaneUTM(hGTIF, psDefn, &oSRS,
3,716✔
889
                                             &linearUnitIsSet))
3,716✔
890
        {
891
            return OGRSpatialReference::ToHandle(oSRS.Clone());
×
892
        }
893

894
        /* Handle ESRI PE string in citation */
895
        szCTString[0] = '\0';
3,716✔
896
        if (GDALGTIFKeyGetASCII(hGTIF, GTCitationGeoKey, szCTString,
3,716✔
897
                                sizeof(szCTString)))
3,716✔
898
            SetCitationToSRS(hGTIF, szCTString, sizeof(szCTString),
3,697✔
899
                             GTCitationGeoKey, &oSRS, &linearUnitIsSet);
900
    }
901

902
    /* ==================================================================== */
903
    /*      Setup the GeogCS                                                */
904
    /* ==================================================================== */
905
    char *pszGeogName = nullptr;
7,981✔
906
    char *pszDatumName = nullptr;
7,981✔
907
    char *pszPMName = nullptr;
7,981✔
908
    char *pszSpheroidName = nullptr;
7,981✔
909
    char *pszAngularUnits = nullptr;
7,981✔
910
    char szGCSName[512] = {'\0'};
7,981✔
911

912
    if (!
7,981✔
913
#if LIBGEOTIFF_VERSION >= 1600
914
        GTIFGetGCSInfoEx(projContext,
15,962✔
915
#else
916
        GTIFGetGCSInfo(
917
#endif
918
                         psDefn->GCS, &pszGeogName, nullptr, nullptr,
7,981✔
919
                         nullptr) &&
8,080✔
920
        GDALGTIFKeyGetASCII(hGTIF, GeogCitationGeoKey, szGCSName,
99✔
921
                            sizeof(szGCSName)))
922
    {
923
        GetGeogCSFromCitation(szGCSName, sizeof(szGCSName), GeogCitationGeoKey,
86✔
924
                              &pszGeogName, &pszDatumName, &pszPMName,
925
                              &pszSpheroidName, &pszAngularUnits);
926
    }
927
    else
928
    {
929
        GTIFToCPLRecycleString(&pszGeogName);
7,895✔
930
    }
931

932
    if (pszGeogName && STARTS_WITH(pszGeogName, "GCS_"))
7,981✔
933
    {
934
        // Morph from potential ESRI name
935
        char *pszOfficialName = GTIFGetEPSGOfficialName(
10✔
936
            hGTIF, PJ_TYPE_GEOGRAPHIC_2D_CRS, pszGeogName);
937
        if (pszOfficialName)
10✔
938
        {
939
            CPLFree(pszGeogName);
×
940
            pszGeogName = pszOfficialName;
×
941
        }
942
    }
943

944
    if (pszDatumName && strchr(pszDatumName, '_'))
7,981✔
945
    {
946
        // Morph from potential ESRI name
947
        char *pszOfficialName = GTIFGetEPSGOfficialName(
22✔
948
            hGTIF, PJ_TYPE_GEODETIC_REFERENCE_FRAME, pszDatumName);
949
        if (pszOfficialName)
22✔
950
        {
951
            CPLFree(pszDatumName);
3✔
952
            pszDatumName = pszOfficialName;
3✔
953
        }
954
    }
955

956
    if (pszSpheroidName && strchr(pszSpheroidName, '_'))
7,981✔
957
    {
958
        // Morph from potential ESRI name
959
        char *pszOfficialName =
960
            GTIFGetEPSGOfficialName(hGTIF, PJ_TYPE_ELLIPSOID, pszSpheroidName);
4✔
961
        if (pszOfficialName)
4✔
962
        {
963
            CPLFree(pszSpheroidName);
1✔
964
            pszSpheroidName = pszOfficialName;
1✔
965
        }
966
    }
967

968
    if (!pszDatumName)
7,981✔
969
    {
970
#if LIBGEOTIFF_VERSION >= 1600
971
        GTIFGetDatumInfoEx(projContext,
7,917✔
972
#else
973
        GTIFGetDatumInfo(
974
#endif
975
                           psDefn->Datum, &pszDatumName, nullptr);
7,917✔
976
        GTIFToCPLRecycleString(&pszDatumName);
7,917✔
977
    }
978

979
    double dfSemiMajor = 0.0;
7,981✔
980
    double dfInvFlattening = 0.0;
7,981✔
981
    if (!pszSpheroidName)
7,981✔
982
    {
983
#if LIBGEOTIFF_VERSION >= 1600
984
        GTIFGetEllipsoidInfoEx(projContext,
7,918✔
985
#else
986
        GTIFGetEllipsoidInfo(
987
#endif
988
                               psDefn->Ellipsoid, &pszSpheroidName, nullptr,
7,918✔
989
                               nullptr);
990
        GTIFToCPLRecycleString(&pszSpheroidName);
7,918✔
991
    }
992
    else
993
    {
994
        CPL_IGNORE_RET_VAL(GDALGTIFKeyGetDOUBLE(hGTIF, GeogSemiMajorAxisGeoKey,
63✔
995
                                                &(psDefn->SemiMajor), 0, 1));
996
        CPL_IGNORE_RET_VAL(GDALGTIFKeyGetDOUBLE(hGTIF, GeogInvFlatteningGeoKey,
63✔
997
                                                &dfInvFlattening, 0, 1));
998
        if (std::isinf(dfInvFlattening))
63✔
999
        {
1000
            // Deal with the non-nominal case of
1001
            // https://github.com/OSGeo/PROJ/issues/2317
1002
            dfInvFlattening = 0;
×
1003
        }
1004
    }
1005
    if (!pszPMName)
7,981✔
1006
    {
1007
#if LIBGEOTIFF_VERSION >= 1600
1008
        GTIFGetPMInfoEx(projContext,
7,897✔
1009
#else
1010
        GTIFGetPMInfo(
1011
#endif
1012
                        psDefn->PM, &pszPMName, nullptr);
7,897✔
1013
        GTIFToCPLRecycleString(&pszPMName);
7,897✔
1014
    }
1015
    else
1016
    {
1017
        CPL_IGNORE_RET_VAL(
84✔
1018
            GDALGTIFKeyGetDOUBLE(hGTIF, GeogPrimeMeridianLongGeoKey,
84✔
1019
                                 &(psDefn->PMLongToGreenwich), 0, 1));
1020
    }
1021

1022
    if (!pszAngularUnits)
7,981✔
1023
    {
1024
#if LIBGEOTIFF_VERSION >= 1600
1025
        GTIFGetUOMAngleInfoEx(projContext,
7,966✔
1026
#else
1027
        GTIFGetUOMAngleInfo(
1028
#endif
1029
                              psDefn->UOMAngle, &pszAngularUnits,
7,966✔
1030
                              &psDefn->UOMAngleInDegrees);
1031
        if (pszAngularUnits == nullptr)
7,966✔
1032
            pszAngularUnits = CPLStrdup("unknown");
11✔
1033
        else
1034
            GTIFToCPLRecycleString(&pszAngularUnits);
7,955✔
1035
    }
1036
    else
1037
    {
1038
        double dfRadians = 0.0;
15✔
1039
        if (GDALGTIFKeyGetDOUBLE(hGTIF, GeogAngularUnitSizeGeoKey, &dfRadians,
15✔
1040
                                 0, 1))
15✔
1041
        {
1042
            psDefn->UOMAngleInDegrees = dfRadians / CPLAtof(SRS_UA_DEGREE_CONV);
4✔
1043
        }
1044
    }
1045

1046
    // Avoid later division by zero.
1047
    if (psDefn->UOMAngleInDegrees == 0)
7,981✔
1048
    {
1049
        CPLError(CE_Warning, CPLE_AppDefined,
1✔
1050
                 "Invalid value for GeogAngularUnitSizeGeoKey.");
1051
        psDefn->UOMAngleInDegrees = 1;
1✔
1052
    }
1053

1054
    dfSemiMajor = psDefn->SemiMajor;
7,981✔
1055
    if (dfSemiMajor == 0.0)
7,981✔
1056
    {
1057
        CPLFree(pszSpheroidName);
11✔
1058
        pszSpheroidName = CPLStrdup("unretrievable - using WGS84");
11✔
1059
        dfSemiMajor = SRS_WGS84_SEMIMAJOR;
11✔
1060
        dfInvFlattening = SRS_WGS84_INVFLATTENING;
11✔
1061
    }
1062
    else if (dfInvFlattening == 0.0 &&
7,970✔
1063
             ((psDefn->SemiMinor / psDefn->SemiMajor) < 0.99999999999999999 ||
7,917✔
1064
              (psDefn->SemiMinor / psDefn->SemiMajor) > 1.00000000000000001))
13✔
1065
    {
1066
        dfInvFlattening =
7,904✔
1067
            OSRCalcInvFlattening(psDefn->SemiMajor, psDefn->SemiMinor);
7,904✔
1068

1069
        /* Take official inverse flattening definition in the WGS84 case */
1070
        if (fabs(dfSemiMajor - SRS_WGS84_SEMIMAJOR) < 1e-10 &&
7,904✔
1071
            fabs(dfInvFlattening - SRS_WGS84_INVFLATTENING) < 1e-10)
4,605✔
1072
            dfInvFlattening = SRS_WGS84_INVFLATTENING;
4,529✔
1073
    }
1074
    if (!pszGeogName || strlen(pszGeogName) == 0)
7,981✔
1075
    {
1076
        CPLFree(pszGeogName);
13✔
1077
        pszGeogName = CPLStrdup(pszDatumName ? pszDatumName : "unknown");
13✔
1078
    }
1079

1080
    oSRS.SetGeogCS(pszGeogName, pszDatumName, pszSpheroidName, dfSemiMajor,
7,981✔
1081
                   dfInvFlattening, pszPMName, psDefn->PMLongToGreenwich,
1082
                   pszAngularUnits,
1083
                   psDefn->UOMAngleInDegrees * CPLAtof(SRS_UA_DEGREE_CONV));
7,981✔
1084

1085
    bool bGeog3DCRS = false;
7,980✔
1086
    bool bSetDatumEllipsoidCode = true;
7,980✔
1087
    bool bHasWarnedInconsistentGeogCRSEPSG = false;
7,980✔
1088
    {
1089
        const int nGCS = psDefn->GCS;
7,980✔
1090
        if (nGCS != KvUserDefined && nGCS > 0 &&
7,980✔
1091
            psDefn->Model != ModelTypeGeocentric)
7,881✔
1092
        {
1093
            OGRSpatialReference oSRSGeog;
15,763✔
1094
            const bool bGCSCodeValid =
1095
                oSRSGeog.importFromEPSG(nGCS) == OGRERR_NONE;
7,881✔
1096

1097
            const std::string osGTiffSRSSource =
1098
                CPLGetConfigOption("GTIFF_SRS_SOURCE", "");
15,764✔
1099

1100
            // GeoTIFF 1.0 might put a Geographic 3D code in GeodeticCRSGeoKey
1101
            bool bTryCompareToEPSG = oSRSGeog.GetAxesCount() == 2;
7,882✔
1102

1103
            if (psDefn->Datum != KvUserDefined)
7,882✔
1104
            {
1105
                char szCode[12];
1106
                snprintf(szCode, sizeof(szCode), "%d", psDefn->Datum);
7,882✔
1107
                auto ctx = static_cast<PJ_CONTEXT *>(
1108
                    GTIFGetPROJContext(hGTIF, true, nullptr));
7,882✔
1109
                auto datum = proj_create_from_database(
7,882✔
1110
                    ctx, "EPSG", szCode, PJ_CATEGORY_DATUM, 0, nullptr);
1111
                if (datum)
7,882✔
1112
                {
1113
                    if (proj_get_type(datum) ==
7,881✔
1114
                        PJ_TYPE_DYNAMIC_GEODETIC_REFERENCE_FRAME)
1115
                    {
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
                        bTryCompareToEPSG = false;
×
1120
                    }
1121
                    proj_destroy(datum);
7,881✔
1122
                }
1123
            }
1124

1125
            if (bTryCompareToEPSG && !oSRSGeog.IsSameGeogCS(&oSRS) &&
7,886✔
1126
                osGTiffSRSSource.empty())
4✔
1127
            {
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
                CPLError(
2✔
1133
                    CE_Warning, CPLE_AppDefined,
1134
                    "The definition of geographic CRS EPSG:%d got from GeoTIFF "
1135
                    "keys "
1136
                    "is not the same as the one from the EPSG registry, "
1137
                    "which may cause issues during reprojection operations. "
1138
                    "Set GTIFF_SRS_SOURCE configuration option to EPSG to "
1139
                    "use official parameters (overriding the ones from GeoTIFF "
1140
                    "keys), "
1141
                    "or to GEOKEYS to use custom values from GeoTIFF keys "
1142
                    "and drop the EPSG code.",
1143
                    nGCS);
1144
                bHasWarnedInconsistentGeogCRSEPSG = true;
2✔
1145
            }
1146
            if (EQUAL(osGTiffSRSSource.c_str(), "EPSG"))
7,882✔
1147
            {
1148
                oSRS.CopyGeogCSFrom(&oSRSGeog);
1✔
1149
            }
1150
            else if (osGTiffSRSSource.empty() && oSRSGeog.IsDynamic() &&
12,394✔
1151
                     psDefn->Model == ModelTypeGeographic)
4,513✔
1152
            {
1153
                // We should perhaps be more careful and detect overrides
1154
                // of geokeys...
1155
                oSRS = oSRSGeog;
4,220✔
1156
                bSetDatumEllipsoidCode = false;
4,220✔
1157
            }
1158
            else if (bGCSCodeValid && osGTiffSRSSource.empty())
3,661✔
1159
            {
1160
                oSRS.SetAuthority("GEOGCS", "EPSG", nGCS);
3,660✔
1161
            }
1162
            else
1163
            {
1164
                bSetDatumEllipsoidCode = false;
1✔
1165
            }
1166

1167
            int nVertSRSCode = verticalCSType;
7,882✔
1168
            if (verticalDatum == 6030 && nGCS == 4326)  // DatumE_WGS84
7,882✔
1169
            {
1170
                nVertSRSCode = 4979;
1✔
1171
            }
1172

1173
            // Try to reconstruct a Geographic3D CRS from the
1174
            // GeodeticCRSGeoKey and the VerticalGeoKey, when they are
1175
            // consistent
1176
            if (nVertSRSCode > 0 && nVertSRSCode != KvUserDefined)
7,882✔
1177
            {
1178
                OGRSpatialReference oTmpVertSRS;
52✔
1179
                if (oSRSGeog.IsGeographic() && oSRSGeog.GetAxesCount() == 2 &&
52✔
1180
                    oTmpVertSRS.importFromEPSG(nVertSRSCode) == OGRERR_NONE &&
26✔
1181
                    oTmpVertSRS.IsGeographic() &&
58✔
1182
                    oTmpVertSRS.GetAxesCount() == 3)
6✔
1183
                {
1184
                    const char *pszTmpCode =
1185
                        oSRSGeog.GetAuthorityCode("GEOGCS|DATUM");
6✔
1186
                    const char *pszTmpVertCode =
1187
                        oTmpVertSRS.GetAuthorityCode("GEOGCS|DATUM");
6✔
1188
                    if (pszTmpCode && pszTmpVertCode &&
6✔
1189
                        atoi(pszTmpCode) == atoi(pszTmpVertCode))
6✔
1190
                    {
1191
                        verticalCSType = 0;
6✔
1192
                        verticalDatum = 0;
6✔
1193
                        verticalUnits = 0;
6✔
1194
                        oSRS.CopyGeogCSFrom(&oTmpVertSRS);
6✔
1195
                        bSetDatumEllipsoidCode = false;
6✔
1196
                        bGeog3DCRS = true;
6✔
1197
                    }
1198
                }
1199
            }
1200
        }
1201
    }
1202
    if (bSetDatumEllipsoidCode)
7,981✔
1203
    {
1204
        if (psDefn->Datum != KvUserDefined)
3,758✔
1205
            oSRS.SetAuthority("DATUM", "EPSG", psDefn->Datum);
3,681✔
1206

1207
        if (psDefn->Ellipsoid != KvUserDefined)
3,758✔
1208
            oSRS.SetAuthority("SPHEROID", "EPSG", psDefn->Ellipsoid);
3,684✔
1209
    }
1210

1211
    CPLFree(pszGeogName);
7,981✔
1212
    CPLFree(pszDatumName);
7,981✔
1213
    CPLFree(pszSpheroidName);
7,981✔
1214
    CPLFree(pszPMName);
7,981✔
1215
    CPLFree(pszAngularUnits);
7,981✔
1216

1217
    /* -------------------------------------------------------------------- */
1218
    /*      Set projection units if not yet done                            */
1219
    /* -------------------------------------------------------------------- */
1220
    if (psDefn->Model == ModelTypeProjected && !linearUnitIsSet)
7,981✔
1221
    {
1222
        char *pszUnitsName = nullptr;
3,709✔
1223

1224
        if (psDefn->UOMLength != KvUserDefined)
3,709✔
1225
        {
1226
#if LIBGEOTIFF_VERSION >= 1600
1227
            GTIFGetUOMLengthInfoEx(projContext,
3,702✔
1228
#else
1229
            GTIFGetUOMLengthInfo(
1230
#endif
1231
                                   psDefn->UOMLength, &pszUnitsName, nullptr);
3,702✔
1232
        }
1233

1234
        if (pszUnitsName != nullptr)
3,709✔
1235
        {
1236
            char szUOMLength[12];
1237
            snprintf(szUOMLength, sizeof(szUOMLength), "%d", psDefn->UOMLength);
3,702✔
1238
            oSRS.SetTargetLinearUnits(nullptr, pszUnitsName,
3,702✔
1239
                                      psDefn->UOMLengthInMeters, "EPSG",
1240
                                      szUOMLength);
1241
        }
1242
        else
1243
            oSRS.SetLinearUnits("unknown", psDefn->UOMLengthInMeters);
7✔
1244

1245
        GTIFFreeMemory(pszUnitsName);
3,709✔
1246
    }
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
    unsigned short tmp = 0;
7,981✔
1255
    bool bGotFromEPSG = false;
7,981✔
1256
    if (psDefn->Model == ModelTypeProjected && psDefn->PCS != KvUserDefined &&
3,716✔
1257
        GDALGTIFKeyGetSHORT(hGTIF, ProjectionGeoKey, &tmp, 0, 1) == 0 &&
3,619✔
1258
        GDALGTIFKeyGetSHORT(hGTIF, ProjCoordTransGeoKey, &tmp, 0, 1) == 0 &&
3,618✔
1259
        GDALGTIFKeyGetSHORT(hGTIF, GeographicTypeGeoKey, &tmp, 0, 1) == 0 &&
3,618✔
1260
        GDALGTIFKeyGetSHORT(hGTIF, GeogGeodeticDatumGeoKey, &tmp, 0, 1) == 0 &&
3,523✔
1261
        GDALGTIFKeyGetSHORT(hGTIF, GeogEllipsoidGeoKey, &tmp, 0, 1) == 0 &&
15,218✔
1262
        CPLTestBool(CPLGetConfigOption("GTIFF_IMPORT_FROM_EPSG", "YES")))
3,521✔
1263
    {
1264
        // Save error state as importFromEPSGA() will call CPLReset()
1265
        CPLErrorNum errNo = CPLGetLastErrorNo();
3,521✔
1266
        CPLErr eErr = CPLGetLastErrorType();
3,521✔
1267
        const char *pszTmp = CPLGetLastErrorMsg();
3,521✔
1268
        char *pszLastErrorMsg = CPLStrdup(pszTmp ? pszTmp : "");
3,521✔
1269
        CPLPushErrorHandler(CPLQuietErrorHandler);
3,521✔
1270
        OGRSpatialReference oSRSTmp;
7,042✔
1271
        OGRErr eImportErr = oSRSTmp.importFromEPSG(psDefn->PCS);
3,521✔
1272
        CPLPopErrorHandler();
3,521✔
1273
        // Restore error state
1274
        CPLErrorSetState(eErr, errNo, pszLastErrorMsg);
3,521✔
1275
        CPLFree(pszLastErrorMsg);
3,521✔
1276
        bGotFromEPSG = eImportErr == OGRERR_NONE;
3,521✔
1277

1278
        if (bGotFromEPSG)
3,521✔
1279
        {
1280
            // See #6210. In case there's an overridden linear units, take it
1281
            // into account
1282
            const char *pszUnitsName = nullptr;
3,521✔
1283
            double dfUOMLengthInMeters = oSRS.GetLinearUnits(&pszUnitsName);
3,521✔
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
            if (fabs(dfUOMLengthInMeters - oSRSTmp.GetLinearUnits(nullptr)) >
3,521✔
1289
                1e-15 * dfUOMLengthInMeters)
3,521✔
1290
            {
1291
                CPLDebug("GTiff", "Modify EPSG:%d to have %s linear units...",
2✔
1292
                         psDefn->PCS, pszUnitsName ? pszUnitsName : "unknown");
4✔
1293

1294
                const char *pszUnitAuthorityCode =
1295
                    oSRS.GetAuthorityCode("PROJCS|UNIT");
2✔
1296
                const char *pszUnitAuthorityName =
1297
                    oSRS.GetAuthorityName("PROJCS|UNIT");
2✔
1298

1299
                if (pszUnitsName)
2✔
1300
                    oSRSTmp.SetLinearUnitsAndUpdateParameters(
2✔
1301
                        pszUnitsName, dfUOMLengthInMeters, pszUnitAuthorityCode,
1302
                        pszUnitAuthorityName);
1303
            }
1304

1305
            if (bGeog3DCRS)
3,521✔
1306
            {
1307
                oSRSTmp.CopyGeogCSFrom(&oSRS);
1✔
1308
                oSRSTmp.UpdateCoordinateSystemFromGeogCRS();
1✔
1309
            }
1310
            oSRS = std::move(oSRSTmp);
3,521✔
1311
        }
1312
    }
1313

1314
#if !defined(GEO_NORMALIZE_DISABLE_TOWGS84)
1315
    if (psDefn->TOWGS84Count > 0 && bGotFromEPSG &&
8,004✔
1316
        CPLTestBool(CPLGetConfigOption("OSR_STRIP_TOWGS84", "YES")))
23✔
1317
    {
1318
        CPLDebug("OSR", "TOWGS84 information has been removed. "
20✔
1319
                        "It can be kept by setting the OSR_STRIP_TOWGS84 "
1320
                        "configuration option to NO");
1321
    }
1322
    else if (psDefn->TOWGS84Count > 0 &&
7,981✔
1323
             (!bGotFromEPSG ||
20✔
1324
              !CPLTestBool(CPLGetConfigOption("OSR_STRIP_TOWGS84", "YES"))))
3✔
1325
    {
1326
        if (bGotFromEPSG)
20✔
1327
        {
1328
            double adfTOWGS84[7] = {0.0};
3✔
1329
            CPL_IGNORE_RET_VAL(oSRS.GetTOWGS84(adfTOWGS84));
3✔
1330
            bool bSame = true;
3✔
1331
            for (int i = 0; i < 7; i++)
3✔
1332
            {
1333
                if (fabs(adfTOWGS84[i] - psDefn->TOWGS84[i]) > 1e-5)
3✔
1334
                {
1335
                    bSame = false;
3✔
1336
                    break;
3✔
1337
                }
1338
            }
1339
            if (!bSame)
3✔
1340
            {
1341
                CPLDebug("GTiff",
3✔
1342
                         "Modify EPSG:%d to have "
1343
                         "TOWGS84=%f,%f,%f,%f,%f,%f,%f "
1344
                         "coming from GeogTOWGS84GeoKey, instead of "
1345
                         "%f,%f,%f,%f,%f,%f,%f coming from EPSG",
1346
                         psDefn->PCS, psDefn->TOWGS84[0], psDefn->TOWGS84[1],
3✔
1347
                         psDefn->TOWGS84[2], psDefn->TOWGS84[3],
1348
                         psDefn->TOWGS84[4], psDefn->TOWGS84[5],
1349
                         psDefn->TOWGS84[6], adfTOWGS84[0], adfTOWGS84[1],
1350
                         adfTOWGS84[2], adfTOWGS84[3], adfTOWGS84[4],
1351
                         adfTOWGS84[5], adfTOWGS84[6]);
1352
            }
1353
        }
1354

1355
        oSRS.SetTOWGS84(psDefn->TOWGS84[0], psDefn->TOWGS84[1],
20✔
1356
                        psDefn->TOWGS84[2], psDefn->TOWGS84[3],
1357
                        psDefn->TOWGS84[4], psDefn->TOWGS84[5],
1358
                        psDefn->TOWGS84[6]);
1359
    }
1360
#endif
1361

1362
    /* ==================================================================== */
1363
    /*      Handle projection parameters.                                   */
1364
    /* ==================================================================== */
1365
    if (psDefn->Model == ModelTypeProjected && !bGotFromEPSG)
7,981✔
1366
    {
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
        double adfParam[10] = {0.0};
195✔
1375
        int i = 0;  // Used after for.
195✔
1376

1377
        for (; i < std::min(10, psDefn->nParms); i++)
1,546✔
1378
            adfParam[i] = psDefn->ProjParm[i];
1,351✔
1379

1380
        for (; i < 10; i++)
794✔
1381
            adfParam[i] = 0.0;
599✔
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
        if (GDALGTIFKeyGetSHORT(hGTIF, ProjCoordTransGeoKey, &tmp, 0, 1) &&
195✔
1409
            psDefn->UOMAngleInDegrees != 0 && psDefn->UOMAngleInDegrees != 1 &&
198✔
1410
            CPLTestBool(CPLGetConfigOption(
3✔
1411
                "GTIFF_READ_ANGULAR_PARAMS_IN_DEGREE", "NO")))
1412
        {
1413
            adfParam[0] /= psDefn->UOMAngleInDegrees;
1✔
1414
            adfParam[1] /= psDefn->UOMAngleInDegrees;
1✔
1415
            adfParam[2] /= psDefn->UOMAngleInDegrees;
1✔
1416
            adfParam[3] /= psDefn->UOMAngleInDegrees;
1✔
1417
        }
1418
#endif
1419

1420
        /* --------------------------------------------------------------------
1421
         */
1422
        /*      Translation the fundamental projection. */
1423
        /* --------------------------------------------------------------------
1424
         */
1425
        switch (psDefn->CTProjection)
195✔
1426
        {
1427
            case CT_TransverseMercator:
103✔
1428
                oSRS.SetTM(adfParam[0], adfParam[1], adfParam[4], adfParam[5],
103✔
1429
                           adfParam[6]);
1430
                break;
103✔
1431

1432
            case CT_TransvMercator_SouthOriented:
×
1433
                oSRS.SetTMSO(adfParam[0], adfParam[1], adfParam[4], adfParam[5],
×
1434
                             adfParam[6]);
1435
                break;
×
1436

1437
            case CT_Mercator:
15✔
1438
                // If a lat_ts was specified use 2SP, otherwise use 1SP.
1439
                if (psDefn->ProjParmId[2] == ProjStdParallel1GeoKey)
15✔
1440
                {
1441
                    if (psDefn->ProjParmId[4] == ProjScaleAtNatOriginGeoKey)
3✔
1442
                        CPLError(CE_Warning, CPLE_AppDefined,
2✔
1443
                                 "Mercator projection should not define "
1444
                                 "both StdParallel1 and ScaleAtNatOrigin.  "
1445
                                 "Using StdParallel1 and ignoring "
1446
                                 "ScaleAtNatOrigin.");
1447
                    oSRS.SetMercator2SP(adfParam[2], adfParam[0], adfParam[1],
3✔
1448
                                        adfParam[5], adfParam[6]);
1449
                }
1450
                else
1451
                    oSRS.SetMercator(adfParam[0], adfParam[1], adfParam[4],
12✔
1452
                                     adfParam[5], adfParam[6]);
1453

1454
                // Override hack for google mercator.
1455
                if (psDefn->Projection == 1024 || psDefn->Projection == 9841)
15✔
1456
                {
1457
                    oSRS.SetExtension(
7✔
1458
                        "PROJCS", "PROJ4",
1459
                        "+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 "
1460
                        "+lon_0=0.0 "
1461
                        "+x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null "
1462
                        "+wktext  +no_defs");  // TODO(schwehr): Why 2 spaces?
1463
                }
1464
                break;
15✔
1465

1466
            case CT_ObliqueStereographic:
4✔
1467
                oSRS.SetOS(adfParam[0], adfParam[1], adfParam[4], adfParam[5],
4✔
1468
                           adfParam[6]);
1469
                break;
4✔
1470

1471
            case CT_Stereographic:
×
1472
                oSRS.SetStereographic(adfParam[0], adfParam[1], adfParam[4],
×
1473
                                      adfParam[5], adfParam[6]);
1474
                break;
×
1475

1476
            case CT_ObliqueMercator:  // Hotine.
1✔
1477
                oSRS.SetHOM(adfParam[0], adfParam[1], adfParam[2], adfParam[3],
1✔
1478
                            adfParam[4], adfParam[5], adfParam[6]);
1479
                break;
1✔
1480

1481
            case CT_HotineObliqueMercatorAzimuthCenter:
2✔
1482
                oSRS.SetHOMAC(adfParam[0], adfParam[1], adfParam[2],
2✔
1483
                              adfParam[3], adfParam[4], adfParam[5],
1484
                              adfParam[6]);
1485
                break;
2✔
1486

1487
            case CT_ObliqueMercator_Laborde:
×
1488
                oSRS.SetLOM(adfParam[0], adfParam[1], adfParam[2], adfParam[4],
×
1489
                            adfParam[5], adfParam[6]);
1490
                break;
×
1491

1492
            case CT_EquidistantConic:
1✔
1493
                oSRS.SetEC(adfParam[0], adfParam[1], adfParam[2], adfParam[3],
1✔
1494
                           adfParam[5], adfParam[6]);
1495
                break;
1✔
1496

1497
            case CT_CassiniSoldner:
1✔
1498
                oSRS.SetCS(adfParam[0], adfParam[1], adfParam[5], adfParam[6]);
1✔
1499
                break;
1✔
1500

1501
            case CT_Polyconic:
1✔
1502
                oSRS.SetPolyconic(adfParam[0], adfParam[1], adfParam[5],
1✔
1503
                                  adfParam[6]);
1504
                break;
1✔
1505

1506
            case CT_AzimuthalEquidistant:
16✔
1507
                oSRS.SetAE(adfParam[0], adfParam[1], adfParam[5], adfParam[6]);
16✔
1508
                break;
16✔
1509

1510
            case CT_MillerCylindrical:
1✔
1511
                oSRS.SetMC(adfParam[0], adfParam[1], adfParam[5], adfParam[6]);
1✔
1512
                break;
1✔
1513

1514
            case CT_Equirectangular:
5✔
1515
                oSRS.SetEquirectangular2(adfParam[0], adfParam[1], adfParam[2],
5✔
1516
                                         adfParam[5], adfParam[6]);
1517
                break;
5✔
1518

1519
            case CT_Gnomonic:
1✔
1520
                oSRS.SetGnomonic(adfParam[0], adfParam[1], adfParam[5],
1✔
1521
                                 adfParam[6]);
1522
                break;
1✔
1523

1524
            case CT_LambertAzimEqualArea:
2✔
1525
                oSRS.SetLAEA(adfParam[0], adfParam[1], adfParam[5],
2✔
1526
                             adfParam[6]);
1527
                break;
2✔
1528

1529
            case CT_Orthographic:
×
1530
                oSRS.SetOrthographic(adfParam[0], adfParam[1], adfParam[5],
×
1531
                                     adfParam[6]);
1532
                break;
×
1533

1534
            case CT_Robinson:
1✔
1535
                oSRS.SetRobinson(adfParam[1], adfParam[5], adfParam[6]);
1✔
1536
                break;
1✔
1537

1538
            case CT_Sinusoidal:
5✔
1539
                oSRS.SetSinusoidal(adfParam[1], adfParam[5], adfParam[6]);
5✔
1540
                break;
5✔
1541

1542
            case CT_VanDerGrinten:
1✔
1543
                oSRS.SetVDG(adfParam[1], adfParam[5], adfParam[6]);
1✔
1544
                break;
1✔
1545

1546
            case CT_PolarStereographic:
6✔
1547
                oSRS.SetPS(adfParam[0], adfParam[1], adfParam[4], adfParam[5],
6✔
1548
                           adfParam[6]);
1549
                break;
6✔
1550

1551
            case CT_LambertConfConic_2SP:
15✔
1552
                oSRS.SetLCC(adfParam[2], adfParam[3], adfParam[0], adfParam[1],
15✔
1553
                            adfParam[5], adfParam[6]);
1554
                break;
15✔
1555

1556
            case CT_LambertConfConic_1SP:
5✔
1557
                oSRS.SetLCC1SP(adfParam[0], adfParam[1], adfParam[4],
5✔
1558
                               adfParam[5], adfParam[6]);
1559
                break;
5✔
1560

1561
            case CT_AlbersEqualArea:
5✔
1562
                oSRS.SetACEA(adfParam[0], adfParam[1], adfParam[2], adfParam[3],
5✔
1563
                             adfParam[5], adfParam[6]);
1564
                break;
5✔
1565

1566
            case CT_NewZealandMapGrid:
1✔
1567
                oSRS.SetNZMG(adfParam[0], adfParam[1], adfParam[5],
1✔
1568
                             adfParam[6]);
1569
                break;
1✔
1570

1571
            case CT_CylindricalEqualArea:
1✔
1572
                oSRS.SetCEA(adfParam[0], adfParam[1], adfParam[5], adfParam[6]);
1✔
1573
                break;
1✔
1574
            default:
2✔
1575
                if (oSRS.IsProjected())
2✔
1576
                {
1577
                    const char *pszName = oSRS.GetName();
2✔
1578
                    std::string osName(pszName ? pszName : "unnamed");
4✔
1579
                    oSRS.Clear();
2✔
1580
                    oSRS.SetLocalCS(osName.c_str());
2✔
1581
                }
1582
                break;
2✔
1583
        }
1584
    }
1585

1586
    if (psDefn->Model == ModelTypeProjected && psDefn->PCS != KvUserDefined &&
7,981✔
1587
        !bGotFromEPSG)
3,619✔
1588
    {
1589
        OGRSpatialReference oSRSTest(oSRS);
196✔
1590
        OGRSpatialReference oSRSTmp;
196✔
1591

1592
        const bool bPCSCodeValid =
1593
            oSRSTmp.importFromEPSG(psDefn->PCS) == OGRERR_NONE;
98✔
1594
        oSRSTmp.SetAxisMappingStrategy(OAMS_TRADITIONAL_GIS_ORDER);
98✔
1595

1596
        // Force axis to avoid issues with SRS with northing, easting order
1597
        oSRSTest.SetAxes(nullptr, "X", OAO_East, "Y", OAO_North);
98✔
1598
        oSRSTmp.SetAxes(nullptr, "X", OAO_East, "Y", OAO_North);
98✔
1599

1600
        const std::string osGTiffSRSSource =
1601
            CPLGetConfigOption("GTIFF_SRS_SOURCE", "");
196✔
1602
        const char *const apszOptions[] = {
98✔
1603
            "IGNORE_DATA_AXIS_TO_SRS_AXIS_MAPPING=YES", nullptr};
1604
        if (!bHasWarnedInconsistentGeogCRSEPSG &&
292✔
1605
            !oSRSTmp.IsSame(&oSRS, apszOptions) &&
99✔
1606
            EQUAL(osGTiffSRSSource.c_str(), ""))
1✔
1607
        {
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
            CPLError(
×
1613
                CE_Warning, CPLE_AppDefined,
1614
                "The definition of projected CRS EPSG:%d got from GeoTIFF keys "
1615
                "is not the same as the one from the EPSG registry, "
1616
                "which may cause issues during reprojection operations. "
1617
                "Set GTIFF_SRS_SOURCE configuration option to EPSG to "
1618
                "use official parameters (overriding the ones from GeoTIFF "
1619
                "keys), "
1620
                "or to GEOKEYS to use custom values from GeoTIFF keys "
1621
                "and drop the EPSG code.",
1622
                psDefn->PCS);
×
1623
        }
1624
        if (EQUAL(osGTiffSRSSource.c_str(), "EPSG"))
98✔
1625
        {
1626
            oSRS = std::move(oSRSTmp);
1✔
1627
        }
1628
        else if (bPCSCodeValid && EQUAL(osGTiffSRSSource.c_str(), ""))
97✔
1629
        {
1630
            oSRS.SetAuthority(nullptr, "EPSG", psDefn->PCS);
96✔
1631
        }
1632
    }
1633

1634
    if (oSRS.IsProjected() && oSRS.GetAxesCount() == 2)
7,981✔
1635
    {
1636
        const char *pszProjCRSName = oSRS.GetAttrValue("PROJCS");
3,713✔
1637
        if (pszProjCRSName)
3,713✔
1638
        {
1639
            // Hack to be able to read properly what we have written for
1640
            // ESRI:102113 (ESRI ancient WebMercator).
1641
            if (EQUAL(pszProjCRSName, "WGS_1984_Web_Mercator"))
3,713✔
1642
                oSRS.SetFromUserInput("ESRI:102113");
2✔
1643
            // And for EPSG:900913
1644
            else if (EQUAL(pszProjCRSName, "Google Maps Global Mercator"))
3,711✔
1645
                oSRS.importFromEPSG(900913);
×
1646
            else if (strchr(pszProjCRSName, '_'))
3,711✔
1647
            {
1648
                // Morph from potential ESRI name
1649
                char *pszOfficialName = GTIFGetEPSGOfficialName(
11✔
1650
                    hGTIF, PJ_TYPE_PROJECTED_CRS, pszProjCRSName);
1651
                if (pszOfficialName)
11✔
1652
                {
1653
                    oSRS.SetProjCS(pszOfficialName);
×
1654
                    CPLFree(pszOfficialName);
×
1655
                }
1656
            }
1657
        }
1658
    }
1659

1660
    /* ==================================================================== */
1661
    /*      Handle vertical coordinate system information if we have it.    */
1662
    /* ==================================================================== */
1663
    bool bNeedManualVertCS = false;
7,980✔
1664
    char citation[2048] = {'\0'};
7,980✔
1665

1666
    // See https://github.com/OSGeo/gdal/pull/4197
1667
    if (verticalCSType > KvUserDefined || verticalDatum > KvUserDefined ||
7,980✔
1668
        verticalUnits > KvUserDefined)
7,981✔
1669
    {
1670
        CPLError(CE_Warning, CPLE_AppDefined,
×
1671
                 "At least one of VerticalCSTypeGeoKey, VerticalDatumGeoKey or "
1672
                 "VerticalUnitsGeoKey has a value in the private user range. "
1673
                 "Ignoring vertical information.");
1674
        verticalCSType = 0;
1✔
1675
        verticalDatum = 0;
1✔
1676
        verticalUnits = 0;
1✔
1677
    }
1678

1679
    if ((verticalCSType != 0 || verticalDatum != 0 || verticalUnits != 0) &&
8,013✔
1680
        (oSRS.IsGeographic() || oSRS.IsProjected() || oSRS.IsLocal()))
32✔
1681
    {
1682
        std::string osVertCRSName;
64✔
1683
        if (GDALGTIFKeyGetASCII(hGTIF, VerticalCitationGeoKey, citation,
32✔
1684
                                sizeof(citation)))
32✔
1685
        {
1686
            if (STARTS_WITH_CI(citation, "VCS Name = "))
7✔
1687
            {
1688
                memmove(citation, citation + strlen("VCS Name = "),
1✔
1689
                        strlen(citation + strlen("VCS Name = ")) + 1);
1✔
1690
                char *pszPipeChar = strchr(citation, '|');
1✔
1691
                if (pszPipeChar)
1✔
1692
                    *pszPipeChar = '\0';
1✔
1693
                osVertCRSName = citation;
1✔
1694
            }
1695
        }
1696

1697
        OGRSpatialReference oVertSRS;
64✔
1698
        bool bCanBuildCompoundCRS = oSRS.GetRoot() != nullptr;
32✔
1699
        if (verticalCSType != KvUserDefined && verticalCSType > 0)
32✔
1700
        {
1701
            if (!(oVertSRS.importFromEPSG(verticalCSType) == OGRERR_NONE &&
40✔
1702
                  oVertSRS.IsVertical()))
20✔
1703
            {
1704
                bCanBuildCompoundCRS = false;
×
1705
            }
1706
            else
1707
            {
1708
                osVertCRSName = oVertSRS.GetName();
20✔
1709
            }
1710
        }
1711
        if (osVertCRSName.empty())
32✔
1712
            osVertCRSName = "unknown";
11✔
1713

1714
        if (bCanBuildCompoundCRS)
32✔
1715
        {
1716
            const bool bHorizontalHasCode =
1717
                oSRS.GetAuthorityCode(nullptr) != nullptr;
32✔
1718
            const char *pszHorizontalName = oSRS.GetName();
32✔
1719
            const std::string osHorizontalName(
1720
                pszHorizontalName ? pszHorizontalName : "unnamed");
64✔
1721
            /* --------------------------------------------------------------------
1722
             */
1723
            /*      Promote to being a compound coordinate system. */
1724
            /* --------------------------------------------------------------------
1725
             */
1726
            OGR_SRSNode *poOldRoot = oSRS.GetRoot()->Clone();
32✔
1727

1728
            oSRS.Clear();
32✔
1729

1730
            /* --------------------------------------------------------------------
1731
             */
1732
            /*      Set COMPD_CS name. */
1733
            /* --------------------------------------------------------------------
1734
             */
1735
            char szCTString[512];
1736
            szCTString[0] = '\0';
32✔
1737
            if (GDALGTIFKeyGetASCII(hGTIF, GTCitationGeoKey, szCTString,
32✔
1738
                                    sizeof(szCTString)) &&
57✔
1739
                strstr(szCTString, " = ") == nullptr)
25✔
1740
            {
1741
                oSRS.SetNode("COMPD_CS", szCTString);
25✔
1742
            }
1743
            else
1744
            {
1745
                oSRS.SetNode(
7✔
1746
                    "COMPD_CS",
1747
                    (osHorizontalName + " + " + osVertCRSName).c_str());
14✔
1748
            }
1749

1750
            oSRS.GetRoot()->AddChild(poOldRoot);
32✔
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
            bNeedManualVertCS = true;
32✔
1760

1761
            if (!oVertSRS.IsEmpty())
32✔
1762
            {
1763
                oSRS.GetRoot()->AddChild(oVertSRS.GetRoot()->Clone());
20✔
1764
                bNeedManualVertCS = false;
20✔
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
                if (bHorizontalHasCode && verticalCSType != KvUserDefined &&
20✔
1771
                    verticalCSType > 0)
20✔
1772
                {
1773
                    const auto *poSRSMatch = oSRS.FindBestMatch(100);
20✔
1774
                    if (poSRSMatch)
20✔
1775
                        oSRS = *poSRSMatch;
3✔
1776
                    delete poSRSMatch;
20✔
1777
                }
1778
            }
1779
        }
1780
    }
1781

1782
    /* -------------------------------------------------------------------- */
1783
    /*      Collect some information from the VerticalCS if not provided    */
1784
    /*      via geokeys.                                                    */
1785
    /* -------------------------------------------------------------------- */
1786
    if (bNeedManualVertCS)
7,981✔
1787
    {
1788
        FillCompoundCRSWithManualVertCS(hGTIF, oSRS, citation, verticalDatum,
12✔
1789
                                        verticalUnits);
1790
    }
1791

1792
    // Hack for tiff_read.py:test_tiff_grads so as to normalize angular
1793
    // parameters to grad
1794
    if (psDefn->UOMAngleInDegrees != 1.0)
7,981✔
1795
    {
1796
        char *pszWKT = nullptr;
21✔
1797
        const char *const apszOptions[] = {
21✔
1798
            "FORMAT=WKT1", "ADD_TOWGS84_ON_EXPORT_TO_WKT1=NO", nullptr};
1799
        if (oSRS.exportToWkt(&pszWKT, apszOptions) == OGRERR_NONE)
21✔
1800
        {
1801
            const char *const apszOptionsImport[] = {
21✔
1802
#if PROJ_AT_LEAST_VERSION(9, 1, 0)
1803
                "UNSET_IDENTIFIERS_IF_INCOMPATIBLE_DEF=NO",
1804
#endif
1805
                nullptr
1806
            };
1807
            oSRS.importFromWkt(pszWKT, apszOptionsImport);
21✔
1808
        }
1809
        CPLFree(pszWKT);
21✔
1810
    }
1811

1812
    oSRS.StripTOWGS84IfKnownDatumAndAllowed();
7,981✔
1813

1814
    double dfCoordinateEpoch = 0.0;
7,981✔
1815
    if (GDALGTIFKeyGetDOUBLE(hGTIF, CoordinateEpochGeoKey, &dfCoordinateEpoch,
7,981✔
1816
                             0, 1))
7,981✔
1817
    {
1818
        oSRS.SetCoordinateEpoch(dfCoordinateEpoch);
2✔
1819
    }
1820

1821
    return OGRSpatialReference::ToHandle(oSRS.Clone());
7,981✔
1822
}
1823

1824
/************************************************************************/
1825
/*                          GTIFGetOGISDefn()                           */
1826
/************************************************************************/
1827

1828
char *GTIFGetOGISDefn(GTIF *hGTIF, GTIFDefn *psDefn)
×
1829
{
1830
    OGRSpatialReferenceH hSRS = GTIFGetOGISDefnAsOSR(hGTIF, psDefn);
×
1831

1832
    char *pszWKT = nullptr;
×
1833
    if (hSRS && OGRSpatialReference::FromHandle(hSRS)->exportToWkt(&pszWKT) ==
×
1834
                    OGRERR_NONE)
1835
    {
1836
        OSRDestroySpatialReference(hSRS);
×
1837
        return pszWKT;
×
1838
    }
1839
    CPLFree(pszWKT);
×
1840
    OSRDestroySpatialReference(hSRS);
×
1841

1842
    return nullptr;
×
1843
}
1844

1845
/************************************************************************/
1846
/*                     OGCDatumName2EPSGDatumCode()                     */
1847
/************************************************************************/
1848

1849
static int OGCDatumName2EPSGDatumCode(GTIF *psGTIF, const char *pszOGCName)
508✔
1850

1851
{
1852
    int nReturn = KvUserDefined;
508✔
1853

1854
    /* -------------------------------------------------------------------- */
1855
    /*      Do we know it as a built in?                                    */
1856
    /* -------------------------------------------------------------------- */
1857
    if (EQUAL(pszOGCName, "NAD27") ||
508✔
1858
        EQUAL(pszOGCName, "North_American_Datum_1927"))
508✔
1859
        return Datum_North_American_Datum_1927;
×
1860
    else if (EQUAL(pszOGCName, "NAD83") ||
508✔
1861
             EQUAL(pszOGCName, "North_American_Datum_1983"))
508✔
1862
        return Datum_North_American_Datum_1983;
×
1863
    else if (EQUAL(pszOGCName, "WGS84") || EQUAL(pszOGCName, "WGS_1984") ||
508✔
1864
             EQUAL(pszOGCName, "WGS 84"))
46✔
1865
        return Datum_WGS84;
462✔
1866
    else if (EQUAL(pszOGCName, "WGS72") || EQUAL(pszOGCName, "WGS_1972"))
46✔
1867
        return Datum_WGS72;
×
1868

1869
    /* Search in database */
1870
    auto ctx =
1871
        static_cast<PJ_CONTEXT *>(GTIFGetPROJContext(psGTIF, true, nullptr));
46✔
1872
    const PJ_TYPE searchType = PJ_TYPE_GEODETIC_REFERENCE_FRAME;
46✔
1873
    auto list = proj_create_from_name(ctx, "EPSG", pszOGCName, &searchType, 1,
46✔
1874
                                      true, /* approximate match */
1875
                                      10, nullptr);
1876
    if (list)
46✔
1877
    {
1878
        const auto listSize = proj_list_get_count(list);
46✔
1879
        for (int i = 0; nReturn == KvUserDefined && i < listSize; i++)
59✔
1880
        {
1881
            auto datum = proj_list_get(ctx, list, i);
13✔
1882
            if (datum)
13✔
1883
            {
1884
                const char *pszDatumName = proj_get_name(datum);
13✔
1885
                if (pszDatumName)
13✔
1886
                {
1887
                    char *pszTmp = CPLStrdup(pszDatumName);
13✔
1888
                    WKTMassageDatum(&pszTmp);
13✔
1889
                    if (EQUAL(pszTmp, pszOGCName))
13✔
1890
                    {
1891
                        const char *pszCode = proj_get_id_code(datum, 0);
3✔
1892
                        if (pszCode)
3✔
1893
                        {
1894
                            nReturn = atoi(pszCode);
3✔
1895
                        }
1896
                    }
1897
                    CPLFree(pszTmp);
13✔
1898
                }
1899
            }
1900
            proj_destroy(datum);
13✔
1901
        }
1902
        proj_list_destroy(list);
46✔
1903
    }
1904

1905
    return nReturn;
46✔
1906
}
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
{
1917
    /* -------------------------------------------------------------------- */
1918
    /*      Create an OGRSpatialReference object corresponding to the       */
1919
    /*      string.                                                         */
1920
    /* -------------------------------------------------------------------- */
1921

1922
    OGRSpatialReference oSRS;
×
1923
    if (oSRS.importFromWkt(pszOGCWKT) != OGRERR_NONE)
×
1924
    {
1925
        return FALSE;
×
1926
    }
1927
    return GTIFSetFromOGISDefnEx(psGTIF, OGRSpatialReference::ToHandle(&oSRS),
×
1928
                                 GEOTIFF_KEYS_STANDARD, GEOTIFF_VERSION_1_0);
×
1929
}
1930

1931
int GTIFSetFromOGISDefnEx(GTIF *psGTIF, OGRSpatialReferenceH hSRS,
3,161✔
1932
                          GTIFFKeysFlavorEnum eFlavor,
1933
                          GeoTIFFVersionEnum eVersion)
1934
{
1935
    std::map<geokey_t, std::string> oMapAsciiKeys;
6,322✔
1936

1937
    GTIFKeySet(psGTIF, GTRasterTypeGeoKey, TYPE_SHORT, 1, RasterPixelIsArea);
3,161✔
1938

1939
    const OGRSpatialReference *poSRS = OGRSpatialReference::FromHandle(hSRS);
3,161✔
1940

1941
    /* -------------------------------------------------------------------- */
1942
    /*      Set version number.                                             */
1943
    /* -------------------------------------------------------------------- */
1944
    if (eVersion == GEOTIFF_VERSION_AUTO)
3,161✔
1945
    {
1946
        if (poSRS->IsCompound() ||
5,847✔
1947
            (poSRS->IsGeographic() && poSRS->GetAxesCount() == 3))
2,917✔
1948
        {
1949
            eVersion = GEOTIFF_VERSION_1_1;
15✔
1950
        }
1951
        else
1952
        {
1953
            eVersion = GEOTIFF_VERSION_1_0;
2,915✔
1954
        }
1955
    }
1956
    CPLAssert(eVersion == GEOTIFF_VERSION_1_0 ||
3,161✔
1957
              eVersion == GEOTIFF_VERSION_1_1);
1958
    if (eVersion >= GEOTIFF_VERSION_1_1)
3,161✔
1959
    {
1960
#if LIBGEOTIFF_VERSION >= 1600
1961
        GTIFSetVersionNumbers(psGTIF, GEOTIFF_SPEC_1_1_VERSION,
18✔
1962
                              GEOTIFF_SPEC_1_1_KEY_REVISION,
1963
                              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
    }
1971

1972
    /* -------------------------------------------------------------------- */
1973
    /*      Get the ellipsoid definition.                                   */
1974
    /* -------------------------------------------------------------------- */
1975
    short nSpheroid = KvUserDefined;
3,161✔
1976

1977
    if (poSRS->GetAuthorityName("PROJCS|GEOGCS|DATUM|SPHEROID") != nullptr &&
4,116✔
1978
        EQUAL(poSRS->GetAuthorityName("PROJCS|GEOGCS|DATUM|SPHEROID"), "EPSG"))
955✔
1979
    {
1980
        nSpheroid = static_cast<short>(
955✔
1981
            atoi(poSRS->GetAuthorityCode("PROJCS|GEOGCS|DATUM|SPHEROID")));
955✔
1982
    }
1983
    else if (poSRS->GetAuthorityName("GEOGCS|DATUM|SPHEROID") != nullptr &&
3,790✔
1984
             EQUAL(poSRS->GetAuthorityName("GEOGCS|DATUM|SPHEROID"), "EPSG"))
1,584✔
1985
    {
1986
        nSpheroid = static_cast<short>(
1,584✔
1987
            atoi(poSRS->GetAuthorityCode("GEOGCS|DATUM|SPHEROID")));
1,584✔
1988
    }
1989

1990
    OGRErr eErr = OGRERR_NONE;
3,161✔
1991
    double dfSemiMajor = 0;
3,161✔
1992
    double dfInvFlattening = 0;
3,161✔
1993
    bool bHasEllipsoid = false;
3,161✔
1994
    if (!poSRS->IsLocal())
3,161✔
1995
    {
1996
        bHasEllipsoid = true;
3,159✔
1997
        if (poSRS->IsCompound())
3,159✔
1998
        {
1999
            OGRSpatialReference oSRSTmp(*poSRS);
26✔
2000
            oSRSTmp.StripVertical();
13✔
2001
            bHasEllipsoid = CPL_TO_BOOL(!oSRSTmp.IsLocal());
13✔
2002
        }
2003
        if (bHasEllipsoid)
3,159✔
2004
        {
2005
            dfSemiMajor = poSRS->GetSemiMajor(&eErr);
3,158✔
2006
            dfInvFlattening = poSRS->GetInvFlattening(&eErr);
3,158✔
2007
            if (eErr != OGRERR_NONE)
3,158✔
2008
            {
2009
                dfSemiMajor = 0.0;
111✔
2010
                dfInvFlattening = 0.0;
111✔
2011
            }
2012
        }
2013
    }
2014

2015
    /* -------------------------------------------------------------------- */
2016
    /*      Get the Datum so we can special case a few PCS codes.           */
2017
    /* -------------------------------------------------------------------- */
2018
    int nDatum = KvUserDefined;
3,161✔
2019

2020
    if (poSRS->GetAuthorityName("PROJCS|GEOGCS|DATUM") != nullptr &&
4,117✔
2021
        EQUAL(poSRS->GetAuthorityName("PROJCS|GEOGCS|DATUM"), "EPSG"))
956✔
2022
        nDatum = atoi(poSRS->GetAuthorityCode("PROJCS|GEOGCS|DATUM"));
956✔
2023
    else if (poSRS->GetAuthorityName("GEOGCS|DATUM") != nullptr &&
3,788✔
2024
             EQUAL(poSRS->GetAuthorityName("GEOGCS|DATUM"), "EPSG"))
1,583✔
2025
        nDatum = atoi(poSRS->GetAuthorityCode("GEOGCS|DATUM"));
1,583✔
2026
    else if (poSRS->GetAttrValue("DATUM") != nullptr)
622✔
2027
        nDatum =
2028
            OGCDatumName2EPSGDatumCode(psGTIF, poSRS->GetAttrValue("DATUM"));
508✔
2029

2030
    /* -------------------------------------------------------------------- */
2031
    /*      Get the GCS if possible.                                        */
2032
    /* -------------------------------------------------------------------- */
2033
    int nGCS = KvUserDefined;
3,161✔
2034

2035
    if (poSRS->GetAuthorityName("PROJCS|GEOGCS") != nullptr &&
4,093✔
2036
        EQUAL(poSRS->GetAuthorityName("PROJCS|GEOGCS"), "EPSG"))
932✔
2037
        nGCS = atoi(poSRS->GetAuthorityCode("PROJCS|GEOGCS"));
932✔
2038
    else if (poSRS->GetAuthorityName("GEOGCS") != nullptr &&
3,812✔
2039
             EQUAL(poSRS->GetAuthorityName("GEOGCS"), "EPSG"))
1,583✔
2040
        nGCS = atoi(poSRS->GetAuthorityCode("GEOGCS"));
1,582✔
2041

2042
    int nVerticalCSKeyValue = 0;
3,161✔
2043
    bool hasEllipsoidHeight = !poSRS->IsCompound() && poSRS->IsGeographic() &&
5,209✔
2044
                              poSRS->GetAxesCount() == 3;
2,048✔
2045
    if (nGCS == 4937 && eVersion >= GEOTIFF_VERSION_1_1)
3,161✔
2046
    {
2047
        // Workaround a bug of PROJ 6.3.0
2048
        // See https://github.com/OSGeo/PROJ/pull/1880
2049
        // EPSG:4937 = ETRS89 3D
2050
        hasEllipsoidHeight = true;
1✔
2051
        nVerticalCSKeyValue = nGCS;
1✔
2052
        nGCS = 4258;  // ETRS89 2D
1✔
2053
    }
2054
    else if (nGCS != KvUserDefined)
3,160✔
2055
    {
2056
        OGRSpatialReference oGeogCRS;
5,026✔
2057
        if (oGeogCRS.importFromEPSG(nGCS) == OGRERR_NONE &&
2,513✔
2058
            oGeogCRS.IsGeographic() && oGeogCRS.GetAxesCount() == 3)
2,513✔
2059
        {
2060
            hasEllipsoidHeight = true;
3✔
2061
            if (eVersion >= GEOTIFF_VERSION_1_1)
3✔
2062
            {
2063
                const auto candidate_nVerticalCSKeyValue = nGCS;
2✔
2064
                nGCS = KvUserDefined;
2✔
2065

2066
                // In case of a geographic 3D CRS, find the corresponding
2067
                // geographic 2D CRS
2068
                auto ctx = static_cast<PJ_CONTEXT *>(
2069
                    GTIFGetPROJContext(psGTIF, true, nullptr));
2✔
2070
                const auto type = PJ_TYPE_GEOGRAPHIC_2D_CRS;
2✔
2071
                auto list = proj_create_from_name(ctx, "EPSG",
2✔
2072
                                                  oGeogCRS.GetName(), &type, 1,
2073
                                                  false,  // exact match
2074
                                                  1,  // result set limit size,
2075
                                                  nullptr);
2076
                if (list && proj_list_get_count(list) == 1)
2✔
2077
                {
2078
                    auto crs2D = proj_list_get(ctx, list, 0);
2✔
2079
                    if (crs2D)
2✔
2080
                    {
2081
                        const char *pszCode = proj_get_id_code(crs2D, 0);
2✔
2082
                        if (pszCode)
2✔
2083
                        {
2084
                            nVerticalCSKeyValue = candidate_nVerticalCSKeyValue;
2✔
2085
                            nGCS = atoi(pszCode);
2✔
2086
                        }
2087
                        proj_destroy(crs2D);
2✔
2088
                    }
2089
                }
2090
                proj_list_destroy(list);
2✔
2091
            }
2092
        }
2093
    }
2094

2095
    // Deprecated way of encoding ellipsoidal height
2096
    if (hasEllipsoidHeight && nVerticalCSKeyValue == 0)
3,161✔
2097
    {
2098
        if (nGCS == 4979 || nDatum == 6326 || nSpheroid == 7030)
1✔
2099
        {
2100
            nVerticalCSKeyValue = 5030;  // WGS_84_ellipsoid
×
2101
            if (nGCS == 4979 || nDatum == 6326)
×
2102
            {
2103
                nGCS = 4326;
×
2104
            }
2105
        }
2106
        else if (nDatum >= 6001 && nDatum <= 6033)
1✔
2107
        {
2108
            nVerticalCSKeyValue = nDatum - 1000;
×
2109
        }
2110
        else if (nSpheroid >= 7001 && nSpheroid <= 7033)
1✔
2111
        {
2112
            nVerticalCSKeyValue = nSpheroid - 2000;
1✔
2113
        }
2114
    }
2115

2116
    if (nGCS > 32767)
3,161✔
2117
        nGCS = KvUserDefined;
×
2118

2119
    /* -------------------------------------------------------------------- */
2120
    /*      Get the linear units.                                           */
2121
    /* -------------------------------------------------------------------- */
2122
    const char *pszLinearUOMNameTmp = nullptr;
3,161✔
2123
    const double dfLinearUOM = poSRS->GetLinearUnits(&pszLinearUOMNameTmp);
3,161✔
2124
    const std::string osLinearUOMName(pszLinearUOMNameTmp ? pszLinearUOMNameTmp
2125
                                                          : "");
6,322✔
2126
    int nUOMLengthCode = 9001;  // Meters.
3,161✔
2127

2128
    if (poSRS->GetAuthorityName("PROJCS|UNIT") != nullptr &&
3,161✔
2129
        EQUAL(poSRS->GetAuthorityName("PROJCS|UNIT"), "EPSG") &&
4,144✔
2130
        poSRS->GetAttrNode("PROJCS|UNIT") != poSRS->GetAttrNode("GEOGCS|UNIT"))
983✔
2131
        nUOMLengthCode = atoi(poSRS->GetAuthorityCode("PROJCS|UNIT"));
983✔
2132
    else if (EQUAL(osLinearUOMName.c_str(), SRS_UL_FOOT) ||
4,356✔
2133
             fabs(dfLinearUOM - GTIFAtof(SRS_UL_FOOT_CONV)) < 0.0000001)
2,178✔
2134
        nUOMLengthCode = 9002;  // International foot.
×
2135
    else if (EQUAL(osLinearUOMName.c_str(), SRS_UL_US_FOOT) ||
4,355✔
2136
             std::abs(dfLinearUOM - GTIFAtof(SRS_UL_US_FOOT_CONV)) < 0.0000001)
2,177✔
2137
        nUOMLengthCode = 9003;  // US survey foot.
2✔
2138
    else if (fabs(dfLinearUOM - 1.0) > 0.00000001)
2,176✔
2139
        nUOMLengthCode = KvUserDefined;
2✔
2140

2141
    /* -------------------------------------------------------------------- */
2142
    /*      Get some authority values.                                      */
2143
    /* -------------------------------------------------------------------- */
2144
    int nPCS = KvUserDefined;
3,161✔
2145

2146
    if (poSRS->GetAuthorityName("PROJCS") != nullptr &&
4,086✔
2147
        EQUAL(poSRS->GetAuthorityName("PROJCS"), "EPSG"))
925✔
2148
    {
2149
        nPCS = atoi(poSRS->GetAuthorityCode("PROJCS"));
923✔
2150
        if (nPCS > 32767)
923✔
2151
            nPCS = KvUserDefined;
×
2152
    }
2153

2154
    /* -------------------------------------------------------------------- */
2155
    /*      Handle the projection transformation.                           */
2156
    /* -------------------------------------------------------------------- */
2157
    bool bWritePEString = false;
3,161✔
2158
    bool bUnknownProjection = false;
3,161✔
2159

2160
    const OGRSpatialReference *poSRSCompatibleOfWKT1 = poSRS;
3,161✔
2161
#if PROJ_AT_LEAST_VERSION(6, 3, 0)
2162
    OGRSpatialReference oSRSCompatibleOfWKT1;
6,322✔
2163
    if (poSRS->IsProjected() && !poSRS->IsCompound() &&
4,146✔
2164
        poSRS->GetAxesCount() == 3)
985✔
2165
    {
2166
        oSRSCompatibleOfWKT1 = *poSRS;
×
2167
        oSRSCompatibleOfWKT1.DemoteTo2D(nullptr);
×
2168
        poSRSCompatibleOfWKT1 = &oSRSCompatibleOfWKT1;
×
2169
        bWritePEString = true;
×
2170
    }
2171
#endif
2172

2173
    double dfAngUnitValue = 0;
3,161✔
2174
    std::string osAngUnitName;
3,161✔
2175
    if (bHasEllipsoid)
3,161✔
2176
    {
2177
        const char *angUnitNameTmp = "";
3,158✔
2178
        dfAngUnitValue = poSRS->GetAngularUnits(&angUnitNameTmp);
3,158✔
2179
        osAngUnitName = angUnitNameTmp;
3,158✔
2180
    }
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
    const auto ConvertAngularParam = [dfAngUnitValue](double dfValInDeg)
262✔
2188
    {
2189
        constexpr double DEG_TO_RAD = M_PI / 180.0;
130✔
2190
        return dfAngUnitValue != 0 &&
130✔
2191
                       std::fabs(dfAngUnitValue - DEG_TO_RAD) > 1e-10 &&
130✔
2192
                       !CPLTestBool(CPLGetConfigOption(
4✔
2193
                           "GTIFF_WRITE_ANGULAR_PARAMS_IN_DEGREE", "NO"))
2194
                   ? dfValInDeg * DEG_TO_RAD / dfAngUnitValue
260✔
2195
                   : dfValInDeg;
130✔
2196
    };
3,161✔
2197

2198
    const char *pszProjection =
2199
        poSRSCompatibleOfWKT1->GetAttrValue("PROJECTION");
3,161✔
2200
    if (nPCS != KvUserDefined)
3,161✔
2201
    {
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
        if (eFlavor == GEOTIFF_KEYS_ESRI_PE && nPCS == 3857)
923✔
2206
        {
2207
            bWritePEString = true;
1✔
2208
        }
2209
        else
2210
        {
2211
            GTIFKeySet(psGTIF, GTModelTypeGeoKey, TYPE_SHORT, 1,
922✔
2212
                       ModelTypeProjected);
2213
            GTIFKeySet(psGTIF, ProjectedCSTypeGeoKey, TYPE_SHORT, 1, nPCS);
922✔
2214
        }
2215
    }
2216
    else if (poSRSCompatibleOfWKT1->IsGeocentric())
2,238✔
2217
    {
2218
        GTIFKeySet(psGTIF, GTModelTypeGeoKey, TYPE_SHORT, 1,
2✔
2219
                   ModelTypeGeocentric);
2220
    }
2221
    else if (pszProjection == nullptr)
2,236✔
2222
    {
2223
        if (poSRSCompatibleOfWKT1->IsGeographic())
2,167✔
2224
            GTIFKeySet(psGTIF, GTModelTypeGeoKey, TYPE_SHORT, 1,
2,053✔
2225
                       ModelTypeGeographic);
2226
        // Otherwise, presumably something like LOCAL_CS.
2227
    }
2228
    else if (EQUAL(pszProjection, SRS_PT_ALBERS_CONIC_EQUAL_AREA))
69✔
2229
    {
2230
        GTIFKeySet(psGTIF, GTModelTypeGeoKey, TYPE_SHORT, 1,
3✔
2231
                   ModelTypeProjected);
2232
        GTIFKeySet(psGTIF, ProjectedCSTypeGeoKey, TYPE_SHORT, 1, KvUserDefined);
3✔
2233
        GTIFKeySet(psGTIF, ProjectionGeoKey, TYPE_SHORT, 1, KvUserDefined);
3✔
2234

2235
        GTIFKeySet(psGTIF, ProjCoordTransGeoKey, TYPE_SHORT, 1,
3✔
2236
                   CT_AlbersEqualArea);
2237

2238
        GTIFKeySet(psGTIF, ProjStdParallelGeoKey, TYPE_DOUBLE, 1,
3✔
2239
                   ConvertAngularParam(poSRSCompatibleOfWKT1->GetNormProjParm(
2240
                       SRS_PP_STANDARD_PARALLEL_1, 0.0)));
2241

2242
        GTIFKeySet(psGTIF, ProjStdParallel2GeoKey, TYPE_DOUBLE, 1,
3✔
2243
                   ConvertAngularParam(poSRSCompatibleOfWKT1->GetNormProjParm(
2244
                       SRS_PP_STANDARD_PARALLEL_2, 0.0)));
2245

2246
        GTIFKeySet(psGTIF, ProjNatOriginLatGeoKey, TYPE_DOUBLE, 1,
3✔
2247
                   ConvertAngularParam(poSRSCompatibleOfWKT1->GetNormProjParm(
2248
                       SRS_PP_LATITUDE_OF_CENTER, 0.0)));
2249

2250
        GTIFKeySet(psGTIF, ProjNatOriginLongGeoKey, TYPE_DOUBLE, 1,
3✔
2251
                   ConvertAngularParam(poSRSCompatibleOfWKT1->GetNormProjParm(
2252
                       SRS_PP_LONGITUDE_OF_CENTER, 0.0)));
2253

2254
        GTIFKeySet(
3✔
2255
            psGTIF, ProjFalseEastingGeoKey, TYPE_DOUBLE, 1,
2256
            poSRSCompatibleOfWKT1->GetProjParm(SRS_PP_FALSE_EASTING, 0.0));
2257

2258
        GTIFKeySet(
3✔
2259
            psGTIF, ProjFalseNorthingGeoKey, TYPE_DOUBLE, 1,
2260
            poSRSCompatibleOfWKT1->GetProjParm(SRS_PP_FALSE_NORTHING, 0.0));
2261
    }
2262

2263
    else if (poSRSCompatibleOfWKT1->GetUTMZone() != 0)
66✔
2264
    {
2265
        GTIFKeySet(psGTIF, GTModelTypeGeoKey, TYPE_SHORT, 1,
10✔
2266
                   ModelTypeProjected);
2267

2268
        int bNorth = 0;
10✔
2269
        const int nZone = poSRSCompatibleOfWKT1->GetUTMZone(&bNorth);
10✔
2270

2271
        if (nDatum == Datum_North_American_Datum_1983 && nZone >= 3 &&
10✔
2272
            nZone <= 22 && bNorth && nUOMLengthCode == 9001)
1✔
2273
        {
2274
            nPCS = 26900 + nZone;
1✔
2275

2276
            GTIFKeySet(psGTIF, ProjectedCSTypeGeoKey, TYPE_SHORT, 1, nPCS);
1✔
2277
        }
2278
        else if (nDatum == Datum_North_American_Datum_1927 && nZone >= 3 &&
9✔
2279
                 nZone <= 22 && bNorth && nUOMLengthCode == 9001)
1✔
2280
        {
2281
            nPCS = 26700 + nZone;
1✔
2282

2283
            GTIFKeySet(psGTIF, ProjectedCSTypeGeoKey, TYPE_SHORT, 1, nPCS);
1✔
2284
        }
2285
        else if (nDatum == Datum_WGS84 && nUOMLengthCode == 9001)
8✔
2286
        {
2287
            if (bNorth)
2✔
2288
                nPCS = 32600 + nZone;
1✔
2289
            else
2290
                nPCS = 32700 + nZone;
1✔
2291

2292
            GTIFKeySet(psGTIF, ProjectedCSTypeGeoKey, TYPE_SHORT, 1, nPCS);
2✔
2293
        }
2294
        else
2295
        {
2296
            const int nProjection = nZone + (bNorth ? 16000 : 16100);
6✔
2297
            GTIFKeySet(psGTIF, ProjectedCSTypeGeoKey, TYPE_SHORT, 1,
6✔
2298
                       KvUserDefined);
2299

2300
            GTIFKeySet(psGTIF, ProjectionGeoKey, TYPE_SHORT, 1, nProjection);
6✔
2301
        }
2302
    }
2303

2304
    else if (EQUAL(pszProjection, SRS_PT_TRANSVERSE_MERCATOR))
56✔
2305
    {
2306
        GTIFKeySet(psGTIF, GTModelTypeGeoKey, TYPE_SHORT, 1,
6✔
2307
                   ModelTypeProjected);
2308
        GTIFKeySet(psGTIF, ProjectedCSTypeGeoKey, TYPE_SHORT, 1, KvUserDefined);
6✔
2309
        GTIFKeySet(psGTIF, ProjectionGeoKey, TYPE_SHORT, 1, KvUserDefined);
6✔
2310

2311
        GTIFKeySet(psGTIF, ProjCoordTransGeoKey, TYPE_SHORT, 1,
6✔
2312
                   CT_TransverseMercator);
2313

2314
        GTIFKeySet(psGTIF, ProjNatOriginLatGeoKey, TYPE_DOUBLE, 1,
6✔
2315
                   ConvertAngularParam(poSRSCompatibleOfWKT1->GetNormProjParm(
2316
                       SRS_PP_LATITUDE_OF_ORIGIN, 0.0)));
2317

2318
        GTIFKeySet(psGTIF, ProjNatOriginLongGeoKey, TYPE_DOUBLE, 1,
6✔
2319
                   ConvertAngularParam(poSRSCompatibleOfWKT1->GetNormProjParm(
2320
                       SRS_PP_CENTRAL_MERIDIAN, 0.0)));
2321

2322
        GTIFKeySet(
6✔
2323
            psGTIF, ProjScaleAtNatOriginGeoKey, TYPE_DOUBLE, 1,
2324
            poSRSCompatibleOfWKT1->GetNormProjParm(SRS_PP_SCALE_FACTOR, 1.0));
2325

2326
        GTIFKeySet(
6✔
2327
            psGTIF, ProjFalseEastingGeoKey, TYPE_DOUBLE, 1,
2328
            poSRSCompatibleOfWKT1->GetProjParm(SRS_PP_FALSE_EASTING, 0.0));
2329

2330
        GTIFKeySet(
6✔
2331
            psGTIF, ProjFalseNorthingGeoKey, TYPE_DOUBLE, 1,
2332
            poSRSCompatibleOfWKT1->GetProjParm(SRS_PP_FALSE_NORTHING, 0.0));
2333
    }
2334

2335
    else if (EQUAL(pszProjection, SRS_PT_TRANSVERSE_MERCATOR_SOUTH_ORIENTED))
50✔
2336
    {
2337
        GTIFKeySet(psGTIF, GTModelTypeGeoKey, TYPE_SHORT, 1,
×
2338
                   ModelTypeProjected);
2339
        GTIFKeySet(psGTIF, ProjectedCSTypeGeoKey, TYPE_SHORT, 1, KvUserDefined);
×
2340
        GTIFKeySet(psGTIF, ProjectionGeoKey, TYPE_SHORT, 1, KvUserDefined);
×
2341

2342
        GTIFKeySet(psGTIF, ProjCoordTransGeoKey, TYPE_SHORT, 1,
×
2343
                   CT_TransvMercator_SouthOriented);
2344

2345
        GTIFKeySet(psGTIF, ProjNatOriginLatGeoKey, TYPE_DOUBLE, 1,
×
2346
                   ConvertAngularParam(poSRSCompatibleOfWKT1->GetNormProjParm(
2347
                       SRS_PP_LATITUDE_OF_ORIGIN, 0.0)));
2348

2349
        GTIFKeySet(psGTIF, ProjNatOriginLongGeoKey, TYPE_DOUBLE, 1,
×
2350
                   ConvertAngularParam(poSRSCompatibleOfWKT1->GetNormProjParm(
2351
                       SRS_PP_CENTRAL_MERIDIAN, 0.0)));
2352

2353
        GTIFKeySet(
×
2354
            psGTIF, ProjScaleAtNatOriginGeoKey, TYPE_DOUBLE, 1,
2355
            poSRSCompatibleOfWKT1->GetNormProjParm(SRS_PP_SCALE_FACTOR, 1.0));
2356

2357
        GTIFKeySet(
×
2358
            psGTIF, ProjFalseEastingGeoKey, TYPE_DOUBLE, 1,
2359
            poSRSCompatibleOfWKT1->GetProjParm(SRS_PP_FALSE_EASTING, 0.0));
2360

2361
        GTIFKeySet(
×
2362
            psGTIF, ProjFalseNorthingGeoKey, TYPE_DOUBLE, 1,
2363
            poSRSCompatibleOfWKT1->GetProjParm(SRS_PP_FALSE_NORTHING, 0.0));
2364
    }
2365

2366
    else if (EQUAL(pszProjection, SRS_PT_MERCATOR_2SP) ||
50✔
2367
             EQUAL(pszProjection, SRS_PT_MERCATOR_1SP))
49✔
2368

2369
    {
2370
        GTIFKeySet(psGTIF, GTModelTypeGeoKey, TYPE_SHORT, 1,
5✔
2371
                   ModelTypeProjected);
2372
        GTIFKeySet(psGTIF, ProjectedCSTypeGeoKey, TYPE_SHORT, 1, KvUserDefined);
5✔
2373
        GTIFKeySet(psGTIF, ProjectionGeoKey, TYPE_SHORT, 1, KvUserDefined);
5✔
2374

2375
        GTIFKeySet(psGTIF, ProjCoordTransGeoKey, TYPE_SHORT, 1, CT_Mercator);
5✔
2376

2377
        GTIFKeySet(psGTIF, ProjNatOriginLatGeoKey, TYPE_DOUBLE, 1,
5✔
2378
                   ConvertAngularParam(poSRSCompatibleOfWKT1->GetNormProjParm(
2379
                       SRS_PP_LATITUDE_OF_ORIGIN, 0.0)));
2380

2381
        GTIFKeySet(psGTIF, ProjNatOriginLongGeoKey, TYPE_DOUBLE, 1,
5✔
2382
                   ConvertAngularParam(poSRSCompatibleOfWKT1->GetNormProjParm(
2383
                       SRS_PP_CENTRAL_MERIDIAN, 0.0)));
2384

2385
        if (EQUAL(pszProjection, SRS_PT_MERCATOR_2SP))
5✔
2386
            GTIFKeySet(
1✔
2387
                psGTIF, ProjStdParallel1GeoKey, TYPE_DOUBLE, 1,
2388
                ConvertAngularParam(poSRSCompatibleOfWKT1->GetNormProjParm(
2389
                    SRS_PP_STANDARD_PARALLEL_1, 0.0)));
2390
        else
2391
            GTIFKeySet(psGTIF, ProjScaleAtNatOriginGeoKey, TYPE_DOUBLE, 1,
4✔
2392
                       poSRSCompatibleOfWKT1->GetNormProjParm(
2393
                           SRS_PP_SCALE_FACTOR, 1.0));
2394

2395
        GTIFKeySet(
5✔
2396
            psGTIF, ProjFalseEastingGeoKey, TYPE_DOUBLE, 1,
2397
            poSRSCompatibleOfWKT1->GetProjParm(SRS_PP_FALSE_EASTING, 0.0));
2398

2399
        GTIFKeySet(
5✔
2400
            psGTIF, ProjFalseNorthingGeoKey, TYPE_DOUBLE, 1,
2401
            poSRSCompatibleOfWKT1->GetProjParm(SRS_PP_FALSE_NORTHING, 0.0));
2402
    }
2403

2404
    else if (EQUAL(pszProjection, SRS_PT_OBLIQUE_STEREOGRAPHIC))
45✔
2405
    {
2406
        GTIFKeySet(psGTIF, GTModelTypeGeoKey, TYPE_SHORT, 1,
1✔
2407
                   ModelTypeProjected);
2408
        GTIFKeySet(psGTIF, ProjectedCSTypeGeoKey, TYPE_SHORT, 1, KvUserDefined);
1✔
2409
        GTIFKeySet(psGTIF, ProjectionGeoKey, TYPE_SHORT, 1, KvUserDefined);
1✔
2410

2411
        GTIFKeySet(psGTIF, ProjCoordTransGeoKey, TYPE_SHORT, 1,
1✔
2412
                   CT_ObliqueStereographic);
2413

2414
        GTIFKeySet(psGTIF, ProjNatOriginLatGeoKey, TYPE_DOUBLE, 1,
1✔
2415
                   ConvertAngularParam(poSRSCompatibleOfWKT1->GetNormProjParm(
2416
                       SRS_PP_LATITUDE_OF_ORIGIN, 0.0)));
2417

2418
        GTIFKeySet(psGTIF, ProjNatOriginLongGeoKey, TYPE_DOUBLE, 1,
1✔
2419
                   ConvertAngularParam(poSRSCompatibleOfWKT1->GetNormProjParm(
2420
                       SRS_PP_CENTRAL_MERIDIAN, 0.0)));
2421

2422
        GTIFKeySet(
1✔
2423
            psGTIF, ProjScaleAtNatOriginGeoKey, TYPE_DOUBLE, 1,
2424
            poSRSCompatibleOfWKT1->GetNormProjParm(SRS_PP_SCALE_FACTOR, 1.0));
2425

2426
        GTIFKeySet(
1✔
2427
            psGTIF, ProjFalseEastingGeoKey, TYPE_DOUBLE, 1,
2428
            poSRSCompatibleOfWKT1->GetProjParm(SRS_PP_FALSE_EASTING, 0.0));
2429

2430
        GTIFKeySet(
1✔
2431
            psGTIF, ProjFalseNorthingGeoKey, TYPE_DOUBLE, 1,
2432
            poSRSCompatibleOfWKT1->GetProjParm(SRS_PP_FALSE_NORTHING, 0.0));
2433
    }
2434

2435
    else if (EQUAL(pszProjection, SRS_PT_STEREOGRAPHIC))
44✔
2436
    {
2437
        GTIFKeySet(psGTIF, GTModelTypeGeoKey, TYPE_SHORT, 1,
×
2438
                   ModelTypeProjected);
2439
        GTIFKeySet(psGTIF, ProjectedCSTypeGeoKey, TYPE_SHORT, 1, KvUserDefined);
×
2440
        GTIFKeySet(psGTIF, ProjectionGeoKey, TYPE_SHORT, 1, KvUserDefined);
×
2441

2442
        GTIFKeySet(psGTIF, ProjCoordTransGeoKey, TYPE_SHORT, 1,
×
2443
                   CT_Stereographic);
2444

2445
        GTIFKeySet(psGTIF, ProjCenterLatGeoKey, TYPE_DOUBLE, 1,
×
2446
                   ConvertAngularParam(poSRSCompatibleOfWKT1->GetNormProjParm(
2447
                       SRS_PP_LATITUDE_OF_ORIGIN, 0.0)));
2448

2449
        GTIFKeySet(psGTIF, ProjCenterLongGeoKey, TYPE_DOUBLE, 1,
×
2450
                   ConvertAngularParam(poSRSCompatibleOfWKT1->GetNormProjParm(
2451
                       SRS_PP_CENTRAL_MERIDIAN, 0.0)));
2452

2453
        GTIFKeySet(
×
2454
            psGTIF, ProjScaleAtNatOriginGeoKey, TYPE_DOUBLE, 1,
2455
            poSRSCompatibleOfWKT1->GetNormProjParm(SRS_PP_SCALE_FACTOR, 1.0));
2456

2457
        GTIFKeySet(
×
2458
            psGTIF, ProjFalseEastingGeoKey, TYPE_DOUBLE, 1,
2459
            poSRSCompatibleOfWKT1->GetProjParm(SRS_PP_FALSE_EASTING, 0.0));
2460

2461
        GTIFKeySet(
×
2462
            psGTIF, ProjFalseNorthingGeoKey, TYPE_DOUBLE, 1,
2463
            poSRSCompatibleOfWKT1->GetProjParm(SRS_PP_FALSE_NORTHING, 0.0));
2464
    }
2465

2466
    else if (EQUAL(pszProjection, SRS_PT_POLAR_STEREOGRAPHIC))
44✔
2467
    {
2468
        GTIFKeySet(psGTIF, GTModelTypeGeoKey, TYPE_SHORT, 1,
4✔
2469
                   ModelTypeProjected);
2470
        GTIFKeySet(psGTIF, ProjectedCSTypeGeoKey, TYPE_SHORT, 1, KvUserDefined);
4✔
2471
        GTIFKeySet(psGTIF, ProjectionGeoKey, TYPE_SHORT, 1, KvUserDefined);
4✔
2472

2473
        GTIFKeySet(psGTIF, ProjCoordTransGeoKey, TYPE_SHORT, 1,
4✔
2474
                   CT_PolarStereographic);
2475

2476
        GTIFKeySet(psGTIF, ProjNatOriginLatGeoKey, TYPE_DOUBLE, 1,
4✔
2477
                   ConvertAngularParam(poSRSCompatibleOfWKT1->GetNormProjParm(
2478
                       SRS_PP_LATITUDE_OF_ORIGIN, 0.0)));
2479

2480
        GTIFKeySet(psGTIF, ProjStraightVertPoleLongGeoKey, TYPE_DOUBLE, 1,
4✔
2481
                   ConvertAngularParam(poSRSCompatibleOfWKT1->GetNormProjParm(
2482
                       SRS_PP_CENTRAL_MERIDIAN, 0.0)));
2483

2484
        GTIFKeySet(
4✔
2485
            psGTIF, ProjScaleAtNatOriginGeoKey, TYPE_DOUBLE, 1,
2486
            poSRSCompatibleOfWKT1->GetNormProjParm(SRS_PP_SCALE_FACTOR, 1.0));
2487

2488
        GTIFKeySet(
4✔
2489
            psGTIF, ProjFalseEastingGeoKey, TYPE_DOUBLE, 1,
2490
            poSRSCompatibleOfWKT1->GetProjParm(SRS_PP_FALSE_EASTING, 0.0));
2491

2492
        GTIFKeySet(
4✔
2493
            psGTIF, ProjFalseNorthingGeoKey, TYPE_DOUBLE, 1,
2494
            poSRSCompatibleOfWKT1->GetProjParm(SRS_PP_FALSE_NORTHING, 0.0));
2495
    }
2496

2497
    else if (EQUAL(pszProjection, SRS_PT_HOTINE_OBLIQUE_MERCATOR))
40✔
2498
    {
2499
        GTIFKeySet(psGTIF, GTModelTypeGeoKey, TYPE_SHORT, 1,
1✔
2500
                   ModelTypeProjected);
2501
        GTIFKeySet(psGTIF, ProjectedCSTypeGeoKey, TYPE_SHORT, 1, KvUserDefined);
1✔
2502
        GTIFKeySet(psGTIF, ProjectionGeoKey, TYPE_SHORT, 1, KvUserDefined);
1✔
2503

2504
        GTIFKeySet(psGTIF, ProjCoordTransGeoKey, TYPE_SHORT, 1,
1✔
2505
                   CT_ObliqueMercator);
2506

2507
        GTIFKeySet(psGTIF, ProjCenterLatGeoKey, TYPE_DOUBLE, 1,
1✔
2508
                   ConvertAngularParam(poSRSCompatibleOfWKT1->GetNormProjParm(
2509
                       SRS_PP_LATITUDE_OF_CENTER, 0.0)));
2510

2511
        GTIFKeySet(psGTIF, ProjCenterLongGeoKey, TYPE_DOUBLE, 1,
1✔
2512
                   ConvertAngularParam(poSRSCompatibleOfWKT1->GetNormProjParm(
2513
                       SRS_PP_LONGITUDE_OF_CENTER, 0.0)));
2514

2515
        GTIFKeySet(psGTIF, ProjAzimuthAngleGeoKey, TYPE_DOUBLE, 1,
1✔
2516
                   poSRSCompatibleOfWKT1->GetNormProjParm(SRS_PP_AZIMUTH, 0.0));
2517

2518
        GTIFKeySet(psGTIF, ProjRectifiedGridAngleGeoKey, TYPE_DOUBLE, 1,
1✔
2519
                   ConvertAngularParam(poSRSCompatibleOfWKT1->GetNormProjParm(
2520
                       SRS_PP_RECTIFIED_GRID_ANGLE, 0.0)));
2521

2522
        GTIFKeySet(
1✔
2523
            psGTIF, ProjScaleAtCenterGeoKey, TYPE_DOUBLE, 1,
2524
            poSRSCompatibleOfWKT1->GetNormProjParm(SRS_PP_SCALE_FACTOR, 1.0));
2525

2526
        GTIFKeySet(
1✔
2527
            psGTIF, ProjFalseEastingGeoKey, TYPE_DOUBLE, 1,
2528
            poSRSCompatibleOfWKT1->GetProjParm(SRS_PP_FALSE_EASTING, 0.0));
2529

2530
        GTIFKeySet(
1✔
2531
            psGTIF, ProjFalseNorthingGeoKey, TYPE_DOUBLE, 1,
2532
            poSRSCompatibleOfWKT1->GetProjParm(SRS_PP_FALSE_NORTHING, 0.0));
2533
    }
2534

2535
    else if (EQUAL(pszProjection,
39✔
2536
                   SRS_PT_HOTINE_OBLIQUE_MERCATOR_AZIMUTH_CENTER))
2537
    {
2538
        GTIFKeySet(psGTIF, GTModelTypeGeoKey, TYPE_SHORT, 1,
2✔
2539
                   ModelTypeProjected);
2540
        GTIFKeySet(psGTIF, ProjectedCSTypeGeoKey, TYPE_SHORT, 1, KvUserDefined);
2✔
2541
        GTIFKeySet(psGTIF, ProjectionGeoKey, TYPE_SHORT, 1, KvUserDefined);
2✔
2542

2543
        GTIFKeySet(psGTIF, ProjCoordTransGeoKey, TYPE_SHORT, 1,
2✔
2544
                   CT_HotineObliqueMercatorAzimuthCenter);
2545

2546
        GTIFKeySet(psGTIF, ProjCenterLatGeoKey, TYPE_DOUBLE, 1,
2✔
2547
                   ConvertAngularParam(poSRSCompatibleOfWKT1->GetNormProjParm(
2548
                       SRS_PP_LATITUDE_OF_CENTER, 0.0)));
2549

2550
        GTIFKeySet(psGTIF, ProjCenterLongGeoKey, TYPE_DOUBLE, 1,
2✔
2551
                   ConvertAngularParam(poSRSCompatibleOfWKT1->GetNormProjParm(
2552
                       SRS_PP_LONGITUDE_OF_CENTER, 0.0)));
2553

2554
        GTIFKeySet(psGTIF, ProjAzimuthAngleGeoKey, TYPE_DOUBLE, 1,
2✔
2555
                   poSRSCompatibleOfWKT1->GetNormProjParm(SRS_PP_AZIMUTH, 0.0));
2556

2557
        GTIFKeySet(psGTIF, ProjRectifiedGridAngleGeoKey, TYPE_DOUBLE, 1,
2✔
2558
                   ConvertAngularParam(poSRSCompatibleOfWKT1->GetNormProjParm(
2559
                       SRS_PP_RECTIFIED_GRID_ANGLE, 0.0)));
2560

2561
        GTIFKeySet(
2✔
2562
            psGTIF, ProjScaleAtCenterGeoKey, TYPE_DOUBLE, 1,
2563
            poSRSCompatibleOfWKT1->GetNormProjParm(SRS_PP_SCALE_FACTOR, 1.0));
2564

2565
        GTIFKeySet(
2✔
2566
            psGTIF, ProjFalseEastingGeoKey, TYPE_DOUBLE, 1,
2567
            poSRSCompatibleOfWKT1->GetProjParm(SRS_PP_FALSE_EASTING, 0.0));
2568

2569
        GTIFKeySet(
2✔
2570
            psGTIF, ProjFalseNorthingGeoKey, TYPE_DOUBLE, 1,
2571
            poSRSCompatibleOfWKT1->GetProjParm(SRS_PP_FALSE_NORTHING, 0.0));
2572
    }
2573

2574
    else if (EQUAL(pszProjection, "Laborde_Oblique_Mercator"))
37✔
2575
    {
2576
        GTIFKeySet(psGTIF, GTModelTypeGeoKey, TYPE_SHORT, 1,
×
2577
                   ModelTypeProjected);
2578
        GTIFKeySet(psGTIF, ProjectedCSTypeGeoKey, TYPE_SHORT, 1, KvUserDefined);
×
2579
        GTIFKeySet(psGTIF, ProjectionGeoKey, TYPE_SHORT, 1, KvUserDefined);
×
2580

2581
        GTIFKeySet(psGTIF, ProjCoordTransGeoKey, TYPE_SHORT, 1,
×
2582
                   CT_ObliqueMercator_Laborde);
2583

2584
        GTIFKeySet(psGTIF, ProjCenterLatGeoKey, TYPE_DOUBLE, 1,
×
2585
                   ConvertAngularParam(poSRSCompatibleOfWKT1->GetNormProjParm(
2586
                       SRS_PP_LATITUDE_OF_CENTER, 0.0)));
2587

2588
        GTIFKeySet(psGTIF, ProjCenterLongGeoKey, TYPE_DOUBLE, 1,
×
2589
                   ConvertAngularParam(poSRSCompatibleOfWKT1->GetNormProjParm(
2590
                       SRS_PP_LONGITUDE_OF_CENTER, 0.0)));
2591

2592
        GTIFKeySet(psGTIF, ProjAzimuthAngleGeoKey, TYPE_DOUBLE, 1,
×
2593
                   poSRSCompatibleOfWKT1->GetNormProjParm(SRS_PP_AZIMUTH, 0.0));
2594

2595
        GTIFKeySet(
×
2596
            psGTIF, ProjScaleAtCenterGeoKey, TYPE_DOUBLE, 1,
2597
            poSRSCompatibleOfWKT1->GetNormProjParm(SRS_PP_SCALE_FACTOR, 1.0));
2598

2599
        GTIFKeySet(
×
2600
            psGTIF, ProjFalseEastingGeoKey, TYPE_DOUBLE, 1,
2601
            poSRSCompatibleOfWKT1->GetProjParm(SRS_PP_FALSE_EASTING, 0.0));
2602

2603
        GTIFKeySet(
×
2604
            psGTIF, ProjFalseNorthingGeoKey, TYPE_DOUBLE, 1,
2605
            poSRSCompatibleOfWKT1->GetProjParm(SRS_PP_FALSE_NORTHING, 0.0));
2606
    }
2607

2608
    else if (EQUAL(pszProjection, SRS_PT_CASSINI_SOLDNER))
37✔
2609
    {
2610
        GTIFKeySet(psGTIF, GTModelTypeGeoKey, TYPE_SHORT, 1,
1✔
2611
                   ModelTypeProjected);
2612
        GTIFKeySet(psGTIF, ProjectedCSTypeGeoKey, TYPE_SHORT, 1, KvUserDefined);
1✔
2613
        GTIFKeySet(psGTIF, ProjectionGeoKey, TYPE_SHORT, 1, KvUserDefined);
1✔
2614

2615
        GTIFKeySet(psGTIF, ProjCoordTransGeoKey, TYPE_SHORT, 1,
1✔
2616
                   CT_CassiniSoldner);
2617

2618
        GTIFKeySet(psGTIF, ProjNatOriginLatGeoKey, TYPE_DOUBLE, 1,
1✔
2619
                   ConvertAngularParam(poSRSCompatibleOfWKT1->GetNormProjParm(
2620
                       SRS_PP_LATITUDE_OF_ORIGIN, 0.0)));
2621

2622
        GTIFKeySet(psGTIF, ProjNatOriginLongGeoKey, TYPE_DOUBLE, 1,
1✔
2623
                   ConvertAngularParam(poSRSCompatibleOfWKT1->GetNormProjParm(
2624
                       SRS_PP_CENTRAL_MERIDIAN, 0.0)));
2625

2626
        GTIFKeySet(
1✔
2627
            psGTIF, ProjFalseEastingGeoKey, TYPE_DOUBLE, 1,
2628
            poSRSCompatibleOfWKT1->GetProjParm(SRS_PP_FALSE_EASTING, 0.0));
2629

2630
        GTIFKeySet(
1✔
2631
            psGTIF, ProjFalseNorthingGeoKey, TYPE_DOUBLE, 1,
2632
            poSRSCompatibleOfWKT1->GetProjParm(SRS_PP_FALSE_NORTHING, 0.0));
2633
    }
2634

2635
    else if (EQUAL(pszProjection, SRS_PT_EQUIDISTANT_CONIC))
36✔
2636
    {
2637
        GTIFKeySet(psGTIF, GTModelTypeGeoKey, TYPE_SHORT, 1,
1✔
2638
                   ModelTypeProjected);
2639
        GTIFKeySet(psGTIF, ProjectedCSTypeGeoKey, TYPE_SHORT, 1, KvUserDefined);
1✔
2640
        GTIFKeySet(psGTIF, ProjectionGeoKey, TYPE_SHORT, 1, KvUserDefined);
1✔
2641

2642
        GTIFKeySet(psGTIF, ProjCoordTransGeoKey, TYPE_SHORT, 1,
1✔
2643
                   CT_EquidistantConic);
2644

2645
        GTIFKeySet(psGTIF, ProjStdParallel1GeoKey, TYPE_DOUBLE, 1,
1✔
2646
                   ConvertAngularParam(poSRSCompatibleOfWKT1->GetNormProjParm(
2647
                       SRS_PP_STANDARD_PARALLEL_1, 0.0)));
2648

2649
        GTIFKeySet(psGTIF, ProjStdParallel2GeoKey, TYPE_DOUBLE, 1,
1✔
2650
                   ConvertAngularParam(poSRSCompatibleOfWKT1->GetNormProjParm(
2651
                       SRS_PP_STANDARD_PARALLEL_2, 0.0)));
2652

2653
        GTIFKeySet(psGTIF, ProjNatOriginLatGeoKey, TYPE_DOUBLE, 1,
1✔
2654
                   ConvertAngularParam(poSRSCompatibleOfWKT1->GetNormProjParm(
2655
                       SRS_PP_LATITUDE_OF_CENTER, 0.0)));
2656

2657
        GTIFKeySet(psGTIF, ProjNatOriginLongGeoKey, TYPE_DOUBLE, 1,
1✔
2658
                   ConvertAngularParam(poSRSCompatibleOfWKT1->GetNormProjParm(
2659
                       SRS_PP_LONGITUDE_OF_CENTER, 0.0)));
2660

2661
        GTIFKeySet(
1✔
2662
            psGTIF, ProjFalseEastingGeoKey, TYPE_DOUBLE, 1,
2663
            poSRSCompatibleOfWKT1->GetProjParm(SRS_PP_FALSE_EASTING, 0.0));
2664

2665
        GTIFKeySet(
1✔
2666
            psGTIF, ProjFalseNorthingGeoKey, TYPE_DOUBLE, 1,
2667
            poSRSCompatibleOfWKT1->GetProjParm(SRS_PP_FALSE_NORTHING, 0.0));
2668
    }
2669

2670
    else if (EQUAL(pszProjection, SRS_PT_POLYCONIC))
35✔
2671
    {
2672
        GTIFKeySet(psGTIF, GTModelTypeGeoKey, TYPE_SHORT, 1,
1✔
2673
                   ModelTypeProjected);
2674
        GTIFKeySet(psGTIF, ProjectedCSTypeGeoKey, TYPE_SHORT, 1, KvUserDefined);
1✔
2675
        GTIFKeySet(psGTIF, ProjectionGeoKey, TYPE_SHORT, 1, KvUserDefined);
1✔
2676

2677
        GTIFKeySet(psGTIF, ProjCoordTransGeoKey, TYPE_SHORT, 1, CT_Polyconic);
1✔
2678

2679
        GTIFKeySet(psGTIF, ProjNatOriginLatGeoKey, TYPE_DOUBLE, 1,
1✔
2680
                   poSRSCompatibleOfWKT1->GetNormProjParm(
2681
                       SRS_PP_LATITUDE_OF_ORIGIN, 0.0));
2682

2683
        GTIFKeySet(psGTIF, ProjNatOriginLongGeoKey, TYPE_DOUBLE, 1,
1✔
2684
                   poSRSCompatibleOfWKT1->GetNormProjParm(
2685
                       SRS_PP_CENTRAL_MERIDIAN, 0.0));
2686

2687
        GTIFKeySet(
1✔
2688
            psGTIF, ProjScaleAtNatOriginGeoKey, TYPE_DOUBLE, 1,
2689
            poSRSCompatibleOfWKT1->GetNormProjParm(SRS_PP_SCALE_FACTOR, 1.0));
2690

2691
        GTIFKeySet(
1✔
2692
            psGTIF, ProjFalseEastingGeoKey, TYPE_DOUBLE, 1,
2693
            poSRSCompatibleOfWKT1->GetProjParm(SRS_PP_FALSE_EASTING, 0.0));
2694

2695
        GTIFKeySet(
1✔
2696
            psGTIF, ProjFalseNorthingGeoKey, TYPE_DOUBLE, 1,
2697
            poSRSCompatibleOfWKT1->GetProjParm(SRS_PP_FALSE_NORTHING, 0.0));
2698
    }
2699

2700
    else if (EQUAL(pszProjection, SRS_PT_AZIMUTHAL_EQUIDISTANT))
34✔
2701
    {
2702
        GTIFKeySet(psGTIF, GTModelTypeGeoKey, TYPE_SHORT, 1,
3✔
2703
                   ModelTypeProjected);
2704
        GTIFKeySet(psGTIF, ProjectedCSTypeGeoKey, TYPE_SHORT, 1, KvUserDefined);
3✔
2705
        GTIFKeySet(psGTIF, ProjectionGeoKey, TYPE_SHORT, 1, KvUserDefined);
3✔
2706

2707
        GTIFKeySet(psGTIF, ProjCoordTransGeoKey, TYPE_SHORT, 1,
3✔
2708
                   CT_AzimuthalEquidistant);
2709

2710
        GTIFKeySet(psGTIF, ProjCenterLatGeoKey, TYPE_DOUBLE, 1,
3✔
2711
                   ConvertAngularParam(poSRSCompatibleOfWKT1->GetNormProjParm(
2712
                       SRS_PP_LATITUDE_OF_CENTER, 0.0)));
2713

2714
        GTIFKeySet(psGTIF, ProjCenterLongGeoKey, TYPE_DOUBLE, 1,
3✔
2715
                   ConvertAngularParam(poSRSCompatibleOfWKT1->GetNormProjParm(
2716
                       SRS_PP_LONGITUDE_OF_CENTER, 0.0)));
2717

2718
        GTIFKeySet(
3✔
2719
            psGTIF, ProjFalseEastingGeoKey, TYPE_DOUBLE, 1,
2720
            poSRSCompatibleOfWKT1->GetProjParm(SRS_PP_FALSE_EASTING, 0.0));
2721

2722
        GTIFKeySet(
3✔
2723
            psGTIF, ProjFalseNorthingGeoKey, TYPE_DOUBLE, 1,
2724
            poSRSCompatibleOfWKT1->GetProjParm(SRS_PP_FALSE_NORTHING, 0.0));
2725
    }
2726

2727
    else if (EQUAL(pszProjection, SRS_PT_MILLER_CYLINDRICAL))
31✔
2728
    {
2729
        GTIFKeySet(psGTIF, GTModelTypeGeoKey, TYPE_SHORT, 1,
1✔
2730
                   ModelTypeProjected);
2731
        GTIFKeySet(psGTIF, ProjectedCSTypeGeoKey, TYPE_SHORT, 1, KvUserDefined);
1✔
2732
        GTIFKeySet(psGTIF, ProjectionGeoKey, TYPE_SHORT, 1, KvUserDefined);
1✔
2733

2734
        GTIFKeySet(psGTIF, ProjCoordTransGeoKey, TYPE_SHORT, 1,
1✔
2735
                   CT_MillerCylindrical);
2736

2737
        GTIFKeySet(psGTIF, ProjCenterLatGeoKey, TYPE_DOUBLE, 1,
1✔
2738
                   ConvertAngularParam(poSRSCompatibleOfWKT1->GetNormProjParm(
2739
                       SRS_PP_LATITUDE_OF_CENTER, 0.0)));
2740

2741
        GTIFKeySet(psGTIF, ProjCenterLongGeoKey, TYPE_DOUBLE, 1,
1✔
2742
                   ConvertAngularParam(poSRSCompatibleOfWKT1->GetNormProjParm(
2743
                       SRS_PP_LONGITUDE_OF_CENTER, 0.0)));
2744

2745
        GTIFKeySet(
1✔
2746
            psGTIF, ProjFalseEastingGeoKey, TYPE_DOUBLE, 1,
2747
            poSRSCompatibleOfWKT1->GetProjParm(SRS_PP_FALSE_EASTING, 0.0));
2748

2749
        GTIFKeySet(
1✔
2750
            psGTIF, ProjFalseNorthingGeoKey, TYPE_DOUBLE, 1,
2751
            poSRSCompatibleOfWKT1->GetProjParm(SRS_PP_FALSE_NORTHING, 0.0));
2752
    }
2753

2754
    else if (EQUAL(pszProjection, SRS_PT_EQUIRECTANGULAR))
30✔
2755
    {
2756
        GTIFKeySet(psGTIF, GTModelTypeGeoKey, TYPE_SHORT, 1,
9✔
2757
                   ModelTypeProjected);
2758
        GTIFKeySet(psGTIF, ProjectedCSTypeGeoKey, TYPE_SHORT, 1, KvUserDefined);
9✔
2759
        GTIFKeySet(psGTIF, ProjectionGeoKey, TYPE_SHORT, 1, KvUserDefined);
9✔
2760

2761
        GTIFKeySet(psGTIF, ProjCoordTransGeoKey, TYPE_SHORT, 1,
9✔
2762
                   CT_Equirectangular);
2763

2764
        GTIFKeySet(psGTIF, ProjCenterLatGeoKey, TYPE_DOUBLE, 1,
9✔
2765
                   ConvertAngularParam(poSRSCompatibleOfWKT1->GetNormProjParm(
2766
                       SRS_PP_LATITUDE_OF_ORIGIN, 0.0)));
2767

2768
        GTIFKeySet(psGTIF, ProjCenterLongGeoKey, TYPE_DOUBLE, 1,
9✔
2769
                   ConvertAngularParam(poSRSCompatibleOfWKT1->GetNormProjParm(
2770
                       SRS_PP_CENTRAL_MERIDIAN, 0.0)));
2771

2772
        GTIFKeySet(psGTIF, ProjStdParallel1GeoKey, TYPE_DOUBLE, 1,
9✔
2773
                   ConvertAngularParam(poSRSCompatibleOfWKT1->GetNormProjParm(
2774
                       SRS_PP_STANDARD_PARALLEL_1, 0.0)));
2775

2776
        GTIFKeySet(
9✔
2777
            psGTIF, ProjFalseEastingGeoKey, TYPE_DOUBLE, 1,
2778
            poSRSCompatibleOfWKT1->GetProjParm(SRS_PP_FALSE_EASTING, 0.0));
2779

2780
        GTIFKeySet(
9✔
2781
            psGTIF, ProjFalseNorthingGeoKey, TYPE_DOUBLE, 1,
2782
            poSRSCompatibleOfWKT1->GetProjParm(SRS_PP_FALSE_NORTHING, 0.0));
2783
    }
2784

2785
    else if (EQUAL(pszProjection, SRS_PT_GNOMONIC))
21✔
2786
    {
2787
        GTIFKeySet(psGTIF, GTModelTypeGeoKey, TYPE_SHORT, 1,
1✔
2788
                   ModelTypeProjected);
2789
        GTIFKeySet(psGTIF, ProjectedCSTypeGeoKey, TYPE_SHORT, 1, KvUserDefined);
1✔
2790
        GTIFKeySet(psGTIF, ProjectionGeoKey, TYPE_SHORT, 1, KvUserDefined);
1✔
2791

2792
        GTIFKeySet(psGTIF, ProjCoordTransGeoKey, TYPE_SHORT, 1, CT_Gnomonic);
1✔
2793

2794
        GTIFKeySet(psGTIF, ProjCenterLatGeoKey, TYPE_DOUBLE, 1,
1✔
2795
                   ConvertAngularParam(poSRSCompatibleOfWKT1->GetNormProjParm(
2796
                       SRS_PP_LATITUDE_OF_ORIGIN, 0.0)));
2797

2798
        GTIFKeySet(psGTIF, ProjCenterLongGeoKey, TYPE_DOUBLE, 1,
1✔
2799
                   ConvertAngularParam(poSRSCompatibleOfWKT1->GetNormProjParm(
2800
                       SRS_PP_CENTRAL_MERIDIAN, 0.0)));
2801

2802
        GTIFKeySet(
1✔
2803
            psGTIF, ProjFalseEastingGeoKey, TYPE_DOUBLE, 1,
2804
            poSRSCompatibleOfWKT1->GetProjParm(SRS_PP_FALSE_EASTING, 0.0));
2805

2806
        GTIFKeySet(
1✔
2807
            psGTIF, ProjFalseNorthingGeoKey, TYPE_DOUBLE, 1,
2808
            poSRSCompatibleOfWKT1->GetProjParm(SRS_PP_FALSE_NORTHING, 0.0));
2809
    }
2810

2811
    else if (EQUAL(pszProjection, SRS_PT_LAMBERT_AZIMUTHAL_EQUAL_AREA))
20✔
2812
    {
2813
        GTIFKeySet(psGTIF, GTModelTypeGeoKey, TYPE_SHORT, 1,
2✔
2814
                   ModelTypeProjected);
2815
        GTIFKeySet(psGTIF, ProjectedCSTypeGeoKey, TYPE_SHORT, 1, KvUserDefined);
2✔
2816
        GTIFKeySet(psGTIF, ProjectionGeoKey, TYPE_SHORT, 1, KvUserDefined);
2✔
2817

2818
        GTIFKeySet(psGTIF, ProjCoordTransGeoKey, TYPE_SHORT, 1,
2✔
2819
                   CT_LambertAzimEqualArea);
2820

2821
        GTIFKeySet(psGTIF, ProjCenterLatGeoKey, TYPE_DOUBLE, 1,
2✔
2822
                   ConvertAngularParam(poSRSCompatibleOfWKT1->GetNormProjParm(
2823
                       SRS_PP_LATITUDE_OF_CENTER, 0.0)));
2824

2825
        GTIFKeySet(psGTIF, ProjCenterLongGeoKey, TYPE_DOUBLE, 1,
2✔
2826
                   ConvertAngularParam(poSRSCompatibleOfWKT1->GetNormProjParm(
2827
                       SRS_PP_LONGITUDE_OF_CENTER, 0.0)));
2828

2829
        GTIFKeySet(
2✔
2830
            psGTIF, ProjFalseEastingGeoKey, TYPE_DOUBLE, 1,
2831
            poSRSCompatibleOfWKT1->GetProjParm(SRS_PP_FALSE_EASTING, 0.0));
2832

2833
        GTIFKeySet(
2✔
2834
            psGTIF, ProjFalseNorthingGeoKey, TYPE_DOUBLE, 1,
2835
            poSRSCompatibleOfWKT1->GetProjParm(SRS_PP_FALSE_NORTHING, 0.0));
2836
    }
2837

2838
    else if (EQUAL(pszProjection, SRS_PT_ORTHOGRAPHIC))
18✔
2839
    {
2840
        GTIFKeySet(psGTIF, GTModelTypeGeoKey, TYPE_SHORT, 1,
×
2841
                   ModelTypeProjected);
2842
        GTIFKeySet(psGTIF, ProjectedCSTypeGeoKey, TYPE_SHORT, 1, KvUserDefined);
×
2843
        GTIFKeySet(psGTIF, ProjectionGeoKey, TYPE_SHORT, 1, KvUserDefined);
×
2844

2845
        GTIFKeySet(psGTIF, ProjCoordTransGeoKey, TYPE_SHORT, 1,
×
2846
                   CT_Orthographic);
2847

2848
        GTIFKeySet(psGTIF, ProjCenterLatGeoKey, TYPE_DOUBLE, 1,
×
2849
                   ConvertAngularParam(poSRSCompatibleOfWKT1->GetNormProjParm(
2850
                       SRS_PP_LATITUDE_OF_ORIGIN, 0.0)));
2851

2852
        GTIFKeySet(psGTIF, ProjCenterLongGeoKey, TYPE_DOUBLE, 1,
×
2853
                   ConvertAngularParam(poSRSCompatibleOfWKT1->GetNormProjParm(
2854
                       SRS_PP_CENTRAL_MERIDIAN, 0.0)));
2855

2856
        GTIFKeySet(
×
2857
            psGTIF, ProjFalseEastingGeoKey, TYPE_DOUBLE, 1,
2858
            poSRSCompatibleOfWKT1->GetProjParm(SRS_PP_FALSE_EASTING, 0.0));
2859

2860
        GTIFKeySet(
×
2861
            psGTIF, ProjFalseNorthingGeoKey, TYPE_DOUBLE, 1,
2862
            poSRSCompatibleOfWKT1->GetProjParm(SRS_PP_FALSE_NORTHING, 0.0));
2863
    }
2864

2865
    else if (EQUAL(pszProjection, SRS_PT_NEW_ZEALAND_MAP_GRID))
18✔
2866
    {
2867
        GTIFKeySet(psGTIF, GTModelTypeGeoKey, TYPE_SHORT, 1,
1✔
2868
                   ModelTypeProjected);
2869
        GTIFKeySet(psGTIF, ProjectedCSTypeGeoKey, TYPE_SHORT, 1, KvUserDefined);
1✔
2870
        GTIFKeySet(psGTIF, ProjectionGeoKey, TYPE_SHORT, 1, KvUserDefined);
1✔
2871

2872
        GTIFKeySet(psGTIF, ProjCoordTransGeoKey, TYPE_SHORT, 1,
1✔
2873
                   CT_NewZealandMapGrid);
2874

2875
        GTIFKeySet(psGTIF, ProjNatOriginLatGeoKey, TYPE_DOUBLE, 1,
1✔
2876
                   ConvertAngularParam(poSRSCompatibleOfWKT1->GetNormProjParm(
2877
                       SRS_PP_LATITUDE_OF_ORIGIN, 0.0)));
2878

2879
        GTIFKeySet(psGTIF, ProjNatOriginLongGeoKey, TYPE_DOUBLE, 1,
1✔
2880
                   ConvertAngularParam(poSRSCompatibleOfWKT1->GetNormProjParm(
2881
                       SRS_PP_CENTRAL_MERIDIAN, 0.0)));
2882

2883
        GTIFKeySet(
1✔
2884
            psGTIF, ProjFalseEastingGeoKey, TYPE_DOUBLE, 1,
2885
            poSRSCompatibleOfWKT1->GetProjParm(SRS_PP_FALSE_EASTING, 0.0));
2886

2887
        GTIFKeySet(
1✔
2888
            psGTIF, ProjFalseNorthingGeoKey, TYPE_DOUBLE, 1,
2889
            poSRSCompatibleOfWKT1->GetProjParm(SRS_PP_FALSE_NORTHING, 0.0));
2890
    }
2891

2892
    else if (EQUAL(pszProjection, SRS_PT_ROBINSON))
17✔
2893
    {
2894
        GTIFKeySet(psGTIF, GTModelTypeGeoKey, TYPE_SHORT, 1,
1✔
2895
                   ModelTypeProjected);
2896
        GTIFKeySet(psGTIF, ProjectedCSTypeGeoKey, TYPE_SHORT, 1, KvUserDefined);
1✔
2897
        GTIFKeySet(psGTIF, ProjectionGeoKey, TYPE_SHORT, 1, KvUserDefined);
1✔
2898

2899
        GTIFKeySet(psGTIF, ProjCoordTransGeoKey, TYPE_SHORT, 1, CT_Robinson);
1✔
2900

2901
        GTIFKeySet(psGTIF, ProjCenterLongGeoKey, TYPE_DOUBLE, 1,
1✔
2902
                   ConvertAngularParam(poSRSCompatibleOfWKT1->GetNormProjParm(
2903
                       SRS_PP_LONGITUDE_OF_CENTER, 0.0)));
2904

2905
        GTIFKeySet(
1✔
2906
            psGTIF, ProjFalseEastingGeoKey, TYPE_DOUBLE, 1,
2907
            poSRSCompatibleOfWKT1->GetProjParm(SRS_PP_FALSE_EASTING, 0.0));
2908

2909
        GTIFKeySet(
1✔
2910
            psGTIF, ProjFalseNorthingGeoKey, TYPE_DOUBLE, 1,
2911
            poSRSCompatibleOfWKT1->GetProjParm(SRS_PP_FALSE_NORTHING, 0.0));
2912
    }
2913

2914
    else if (EQUAL(pszProjection, SRS_PT_SINUSOIDAL))
16✔
2915
    {
2916
        GTIFKeySet(psGTIF, GTModelTypeGeoKey, TYPE_SHORT, 1,
4✔
2917
                   ModelTypeProjected);
2918
        GTIFKeySet(psGTIF, ProjectedCSTypeGeoKey, TYPE_SHORT, 1, KvUserDefined);
4✔
2919
        GTIFKeySet(psGTIF, ProjectionGeoKey, TYPE_SHORT, 1, KvUserDefined);
4✔
2920

2921
        GTIFKeySet(psGTIF, ProjCoordTransGeoKey, TYPE_SHORT, 1, CT_Sinusoidal);
4✔
2922

2923
        GTIFKeySet(psGTIF, ProjCenterLongGeoKey, TYPE_DOUBLE, 1,
4✔
2924
                   ConvertAngularParam(poSRSCompatibleOfWKT1->GetNormProjParm(
2925
                       SRS_PP_LONGITUDE_OF_CENTER, 0.0)));
2926

2927
        GTIFKeySet(
4✔
2928
            psGTIF, ProjFalseEastingGeoKey, TYPE_DOUBLE, 1,
2929
            poSRSCompatibleOfWKT1->GetProjParm(SRS_PP_FALSE_EASTING, 0.0));
2930

2931
        GTIFKeySet(
4✔
2932
            psGTIF, ProjFalseNorthingGeoKey, TYPE_DOUBLE, 1,
2933
            poSRSCompatibleOfWKT1->GetProjParm(SRS_PP_FALSE_NORTHING, 0.0));
2934
    }
2935

2936
    else if (EQUAL(pszProjection, SRS_PT_VANDERGRINTEN))
12✔
2937
    {
2938
        GTIFKeySet(psGTIF, GTModelTypeGeoKey, TYPE_SHORT, 1,
2✔
2939
                   ModelTypeProjected);
2940
        GTIFKeySet(psGTIF, ProjectedCSTypeGeoKey, TYPE_SHORT, 1, KvUserDefined);
2✔
2941
        GTIFKeySet(psGTIF, ProjectionGeoKey, TYPE_SHORT, 1, KvUserDefined);
2✔
2942

2943
        GTIFKeySet(psGTIF, ProjCoordTransGeoKey, TYPE_SHORT, 1,
2✔
2944
                   CT_VanDerGrinten);
2945

2946
        GTIFKeySet(psGTIF, ProjCenterLongGeoKey, TYPE_DOUBLE, 1,
2✔
2947
                   ConvertAngularParam(poSRSCompatibleOfWKT1->GetNormProjParm(
2948
                       SRS_PP_CENTRAL_MERIDIAN, 0.0)));
2949

2950
        GTIFKeySet(
2✔
2951
            psGTIF, ProjFalseEastingGeoKey, TYPE_DOUBLE, 1,
2952
            poSRSCompatibleOfWKT1->GetProjParm(SRS_PP_FALSE_EASTING, 0.0));
2953

2954
        GTIFKeySet(
2✔
2955
            psGTIF, ProjFalseNorthingGeoKey, TYPE_DOUBLE, 1,
2956
            poSRSCompatibleOfWKT1->GetProjParm(SRS_PP_FALSE_NORTHING, 0.0));
2957
    }
2958

2959
    else if (EQUAL(pszProjection, SRS_PT_LAMBERT_CONFORMAL_CONIC_2SP))
10✔
2960
    {
2961
        GTIFKeySet(psGTIF, GTModelTypeGeoKey, TYPE_SHORT, 1,
3✔
2962
                   ModelTypeProjected);
2963
        GTIFKeySet(psGTIF, ProjectedCSTypeGeoKey, TYPE_SHORT, 1, KvUserDefined);
3✔
2964
        GTIFKeySet(psGTIF, ProjectionGeoKey, TYPE_SHORT, 1, KvUserDefined);
3✔
2965

2966
        GTIFKeySet(psGTIF, ProjCoordTransGeoKey, TYPE_SHORT, 1,
3✔
2967
                   CT_LambertConfConic_2SP);
2968

2969
        GTIFKeySet(psGTIF, ProjFalseOriginLatGeoKey, TYPE_DOUBLE, 1,
3✔
2970
                   ConvertAngularParam(poSRSCompatibleOfWKT1->GetNormProjParm(
2971
                       SRS_PP_LATITUDE_OF_ORIGIN, 0.0)));
2972

2973
        GTIFKeySet(psGTIF, ProjFalseOriginLongGeoKey, TYPE_DOUBLE, 1,
3✔
2974
                   ConvertAngularParam(poSRSCompatibleOfWKT1->GetNormProjParm(
2975
                       SRS_PP_CENTRAL_MERIDIAN, 0.0)));
2976

2977
        GTIFKeySet(psGTIF, ProjStdParallel1GeoKey, TYPE_DOUBLE, 1,
3✔
2978
                   ConvertAngularParam(poSRSCompatibleOfWKT1->GetNormProjParm(
2979
                       SRS_PP_STANDARD_PARALLEL_1, 0.0)));
2980

2981
        GTIFKeySet(psGTIF, ProjStdParallel2GeoKey, TYPE_DOUBLE, 1,
3✔
2982
                   ConvertAngularParam(poSRSCompatibleOfWKT1->GetNormProjParm(
2983
                       SRS_PP_STANDARD_PARALLEL_2, 0.0)));
2984

2985
        GTIFKeySet(
3✔
2986
            psGTIF, ProjFalseOriginEastingGeoKey, TYPE_DOUBLE, 1,
2987
            poSRSCompatibleOfWKT1->GetProjParm(SRS_PP_FALSE_EASTING, 0.0));
2988

2989
        GTIFKeySet(
3✔
2990
            psGTIF, ProjFalseOriginNorthingGeoKey, TYPE_DOUBLE, 1,
2991
            poSRSCompatibleOfWKT1->GetProjParm(SRS_PP_FALSE_NORTHING, 0.0));
2992
    }
2993

2994
    else if (EQUAL(pszProjection, SRS_PT_LAMBERT_CONFORMAL_CONIC_1SP))
7✔
2995
    {
2996
        GTIFKeySet(psGTIF, GTModelTypeGeoKey, TYPE_SHORT, 1,
3✔
2997
                   ModelTypeProjected);
2998
        GTIFKeySet(psGTIF, ProjectedCSTypeGeoKey, TYPE_SHORT, 1, KvUserDefined);
3✔
2999
        GTIFKeySet(psGTIF, ProjectionGeoKey, TYPE_SHORT, 1, KvUserDefined);
3✔
3000

3001
        GTIFKeySet(psGTIF, ProjCoordTransGeoKey, TYPE_SHORT, 1,
3✔
3002
                   CT_LambertConfConic_1SP);
3003

3004
        GTIFKeySet(psGTIF, ProjNatOriginLatGeoKey, TYPE_DOUBLE, 1,
3✔
3005
                   ConvertAngularParam(poSRSCompatibleOfWKT1->GetNormProjParm(
3006
                       SRS_PP_LATITUDE_OF_ORIGIN, 0.0)));
3007

3008
        GTIFKeySet(psGTIF, ProjNatOriginLongGeoKey, TYPE_DOUBLE, 1,
3✔
3009
                   ConvertAngularParam(poSRSCompatibleOfWKT1->GetNormProjParm(
3010
                       SRS_PP_CENTRAL_MERIDIAN, 0.0)));
3011

3012
        GTIFKeySet(
3✔
3013
            psGTIF, ProjScaleAtNatOriginGeoKey, TYPE_DOUBLE, 1,
3014
            poSRSCompatibleOfWKT1->GetNormProjParm(SRS_PP_SCALE_FACTOR, 1.0));
3015

3016
        GTIFKeySet(
3✔
3017
            psGTIF, ProjFalseEastingGeoKey, TYPE_DOUBLE, 1,
3018
            poSRSCompatibleOfWKT1->GetProjParm(SRS_PP_FALSE_EASTING, 0.0));
3019

3020
        GTIFKeySet(
3✔
3021
            psGTIF, ProjFalseNorthingGeoKey, TYPE_DOUBLE, 1,
3022
            poSRSCompatibleOfWKT1->GetProjParm(SRS_PP_FALSE_NORTHING, 0.0));
3023
    }
3024

3025
    else if (EQUAL(pszProjection, SRS_PT_CYLINDRICAL_EQUAL_AREA))
4✔
3026
    {
3027
        GTIFKeySet(psGTIF, GTModelTypeGeoKey, TYPE_SHORT, 1,
1✔
3028
                   ModelTypeProjected);
3029
        GTIFKeySet(psGTIF, ProjectedCSTypeGeoKey, TYPE_SHORT, 1, KvUserDefined);
1✔
3030
        GTIFKeySet(psGTIF, ProjectionGeoKey, TYPE_SHORT, 1, KvUserDefined);
1✔
3031

3032
        GTIFKeySet(psGTIF, ProjCoordTransGeoKey, TYPE_SHORT, 1,
1✔
3033
                   CT_CylindricalEqualArea);
3034

3035
        GTIFKeySet(psGTIF, ProjNatOriginLongGeoKey, TYPE_DOUBLE, 1,
1✔
3036
                   ConvertAngularParam(poSRSCompatibleOfWKT1->GetNormProjParm(
3037
                       SRS_PP_CENTRAL_MERIDIAN, 0.0)));
3038

3039
        GTIFKeySet(psGTIF, ProjStdParallel1GeoKey, TYPE_DOUBLE, 1,
1✔
3040
                   ConvertAngularParam(poSRSCompatibleOfWKT1->GetNormProjParm(
3041
                       SRS_PP_STANDARD_PARALLEL_1, 0.0)));
3042

3043
        GTIFKeySet(
1✔
3044
            psGTIF, ProjFalseEastingGeoKey, TYPE_DOUBLE, 1,
3045
            poSRSCompatibleOfWKT1->GetProjParm(SRS_PP_FALSE_EASTING, 0.0));
3046

3047
        GTIFKeySet(
1✔
3048
            psGTIF, ProjFalseNorthingGeoKey, TYPE_DOUBLE, 1,
3049
            poSRSCompatibleOfWKT1->GetProjParm(SRS_PP_FALSE_NORTHING, 0.0));
3050
    }
3051
    else
3052
    {
3053
        bWritePEString = true;
3✔
3054
        bUnknownProjection = true;
3✔
3055
    }
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
    if (eFlavor == GEOTIFF_KEYS_ESRI_PE ||
6,321✔
3060
        poSRS->GetAttrValue("VERTCS") != nullptr)
3,160✔
3061
    {
3062
        bWritePEString = true;
1✔
3063
    }
3064

3065
    if (nPCS == KvUserDefined)
3,161✔
3066
    {
3067
        const char *pszPROJ4Ext =
3068
            poSRS->GetExtension("PROJCS", "PROJ4", nullptr);
2,234✔
3069
        if (pszPROJ4Ext &&
2,234✔
3070
            strstr(pszPROJ4Ext, "+proj=merc +a=6378137 +b=6378137"))
4✔
3071
        {
3072
            bWritePEString = true;
3✔
3073
        }
3074
    }
3075

3076
    bWritePEString &=
3,161✔
3077
        CPLTestBool(CPLGetConfigOption("GTIFF_ESRI_CITATION", "YES"));
3,161✔
3078

3079
    bool peStrStored = false;
3,161✔
3080

3081
    if (bWritePEString)
3,161✔
3082
    {
3083
        // Anything we can't map, store as an ESRI PE string with a citation
3084
        // key.
3085
        char *pszPEString = nullptr;
7✔
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
        const char *const apszOptionsDefault[] = {nullptr};
7✔
3089
        const char *const apszOptionsEsri[] = {"FORMAT=WKT1_ESRI", nullptr};
7✔
3090
        const char *const *papszOptions = apszOptionsDefault;
7✔
3091
        if (!(bUnknownProjection &&
17✔
3092
              poSRS->GetExtension("PROJCS", "PROJ4", nullptr) != nullptr) &&
14✔
3093
            !(poSRS->IsProjected() && !poSRS->IsCompound() &&
7✔
3094
              poSRS->GetAxesCount() == 3))
7✔
3095
        {
3096
            papszOptions = apszOptionsEsri;
7✔
3097
        }
3098
        poSRS->exportToWkt(&pszPEString, papszOptions);
7✔
3099
        const int peStrLen = static_cast<int>(strlen(pszPEString));
7✔
3100
        if (peStrLen > 0)
7✔
3101
        {
3102
            char *outPeStr = static_cast<char *>(
3103
                CPLMalloc(peStrLen + strlen("ESRI PE String = ") + 1));
7✔
3104
            strcpy(outPeStr, "ESRI PE String = ");
7✔
3105
            strcat(outPeStr, pszPEString);
7✔
3106
            oMapAsciiKeys[PCSCitationGeoKey] = outPeStr;
7✔
3107
            peStrStored = true;
7✔
3108
            CPLFree(outPeStr);
7✔
3109
        }
3110
        CPLFree(pszPEString);
7✔
3111
        GTIFKeySet(psGTIF, GTModelTypeGeoKey, TYPE_SHORT, 1, KvUserDefined);
7✔
3112

3113
        // Not completely sure we really need to imitate ArcGIS to that point
3114
        // but that cannot hurt.
3115
        if (nPCS == 3857)
7✔
3116
        {
3117
            oMapAsciiKeys[GTCitationGeoKey] =
1✔
3118
                "PCS Name = WGS_1984_Web_Mercator_Auxiliary_Sphere";
1✔
3119
            GTIFKeySet(psGTIF, GeographicTypeGeoKey, TYPE_SHORT, 1, GCS_WGS_84);
1✔
3120
            GTIFKeySet(psGTIF, GeogSemiMajorAxisGeoKey, TYPE_DOUBLE, 1,
1✔
3121
                       6378137.0);
3122
            GTIFKeySet(psGTIF, GeogInvFlatteningGeoKey, TYPE_DOUBLE, 1,
1✔
3123
                       298.257223563);
3124
        }
3125
    }
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
    double dfFE = 0.0;
3,161✔
3133
    double dfFN = 0.0;
3,161✔
3134

3135
    if (eVersion == GEOTIFF_VERSION_1_0 &&
6,304✔
3136
        (GDALGTIFKeyGetDOUBLE(psGTIF, ProjFalseEastingGeoKey, &dfFE, 0, 1) ||
3,143✔
3137
         GDALGTIFKeyGetDOUBLE(psGTIF, ProjFalseNorthingGeoKey, &dfFN, 0, 1) ||
3,090✔
3138
         GDALGTIFKeyGetDOUBLE(psGTIF, ProjFalseOriginEastingGeoKey, &dfFE, 0,
3,090✔
3139
                              1) ||
3,087✔
3140
         GDALGTIFKeyGetDOUBLE(psGTIF, ProjFalseOriginNorthingGeoKey, &dfFN, 0,
3,087✔
3141
                              1)) &&
56✔
3142
        (dfFE != 0.0 || dfFN != 0.0) && nUOMLengthCode != 9001)
6,304✔
3143
    {
3144
        GTIFKeySet(psGTIF, ProjLinearUnitsInterpCorrectGeoKey, TYPE_SHORT, 1,
3✔
3145
                   static_cast<short>(1));
3146
    }
3147

3148
    /* -------------------------------------------------------------------- */
3149
    /*      Write linear units information.                                 */
3150
    /* -------------------------------------------------------------------- */
3151
    if (poSRS->IsGeocentric())
3,161✔
3152
    {
3153
        GTIFKeySet(psGTIF, GeogLinearUnitsGeoKey, TYPE_SHORT, 1,
2✔
3154
                   nUOMLengthCode);
3155
        if (nUOMLengthCode == KvUserDefined)
2✔
3156
            GTIFKeySet(psGTIF, GeogLinearUnitSizeGeoKey, TYPE_DOUBLE, 1,
×
3157
                       dfLinearUOM);
3158
    }
3159
    else if (!poSRS->IsGeographic() &&
4,086✔
3160
             (nPCS == KvUserDefined || eVersion == GEOTIFF_VERSION_1_0))
927✔
3161
    {
3162
        GTIFKeySet(psGTIF, ProjLinearUnitsGeoKey, TYPE_SHORT, 1,
1,098✔
3163
                   nUOMLengthCode);
3164
        if (nUOMLengthCode == KvUserDefined)
1,098✔
3165
            GTIFKeySet(psGTIF, ProjLinearUnitSizeGeoKey, TYPE_DOUBLE, 1,
2✔
3166
                       dfLinearUOM);
3167

3168
        // If linear units name is available and user defined, store it as
3169
        // citation.
3170
        if (!peStrStored && nUOMLengthCode == KvUserDefined &&
1,091✔
3171
            !osLinearUOMName.empty() &&
2,191✔
3172
            CPLTestBool(CPLGetConfigOption("GTIFF_ESRI_CITATION", "YES")))
2✔
3173
        {
3174
            SetLinearUnitCitation(oMapAsciiKeys, osLinearUOMName.c_str());
2✔
3175
        }
3176
    }
3177

3178
    /* -------------------------------------------------------------------- */
3179
    /*      Write angular units.                                            */
3180
    /* -------------------------------------------------------------------- */
3181

3182
    if (bHasEllipsoid &&
3,161✔
3183
        (nGCS == KvUserDefined || eVersion == GEOTIFF_VERSION_1_0))
2,514✔
3184
    {
3185
        if (EQUAL(osAngUnitName.c_str(), "Degree"))
3,141✔
3186
            GTIFKeySet(psGTIF, GeogAngularUnitsGeoKey, TYPE_SHORT, 1,
3,130✔
3187
                       Angular_Degree);
3188
        else if (EQUAL(osAngUnitName.c_str(), "arc-second"))
11✔
3189
            GTIFKeySet(psGTIF, GeogAngularUnitsGeoKey, TYPE_SHORT, 1,
1✔
3190
                       Angular_Arc_Second);
3191
        else if (EQUAL(osAngUnitName.c_str(), "arc-minute"))
10✔
3192
            GTIFKeySet(psGTIF, GeogAngularUnitsGeoKey, TYPE_SHORT, 1,
1✔
3193
                       Angular_Arc_Minute);
3194
        else if (EQUAL(osAngUnitName.c_str(), "grad"))
9✔
3195
            GTIFKeySet(psGTIF, GeogAngularUnitsGeoKey, TYPE_SHORT, 1,
4✔
3196
                       Angular_Grad);
3197
        else if (EQUAL(osAngUnitName.c_str(), "gon"))
5✔
3198
            GTIFKeySet(psGTIF, GeogAngularUnitsGeoKey, TYPE_SHORT, 1,
1✔
3199
                       Angular_Gon);
3200
        else if (EQUAL(osAngUnitName.c_str(), "radian"))
4✔
3201
            GTIFKeySet(psGTIF, GeogAngularUnitsGeoKey, TYPE_SHORT, 1,
1✔
3202
                       Angular_Radian);
3203
        // else if (EQUAL(osAngUnitName.c_str(), "microradian"))
3204
        //    GTIFKeySet(psGTIF, GeogAngularUnitsGeoKey, TYPE_SHORT, 1,
3205
        //               9109);
3206
        else
3207
        {
3208
            // GeogCitationGeoKey may be rewritten if the gcs is user defined.
3209
            oMapAsciiKeys[GeogCitationGeoKey] = osAngUnitName;
3✔
3210
            GTIFKeySet(psGTIF, GeogAngularUnitSizeGeoKey, TYPE_DOUBLE, 1,
3✔
3211
                       dfAngUnitValue);
3212
        }
3213
    }
3214

3215
    /* -------------------------------------------------------------------- */
3216
    /*      Try to write a citation from the main coordinate system         */
3217
    /*      name.                                                           */
3218
    /* -------------------------------------------------------------------- */
3219
    if (poSRS->GetName() != nullptr &&
6,211✔
3220
        ((poSRS->IsProjected() &&
3,050✔
3221
          (nPCS == KvUserDefined || eVersion == GEOTIFF_VERSION_1_0)) ||
927✔
3222
         poSRS->IsCompound() || poSRS->IsLocal() ||
2,066✔
3223
         (poSRS->IsGeocentric() &&
2,051✔
3224
          (nGCS == KvUserDefined || eVersion == GEOTIFF_VERSION_1_0))))
×
3225
    {
3226
        if (!(bWritePEString && nPCS == 3857))
1,001✔
3227
        {
3228
            oMapAsciiKeys[GTCitationGeoKey] = poSRS->GetName();
1,000✔
3229
        }
3230
    }
3231

3232
    /* -------------------------------------------------------------------- */
3233
    /*      Try to write a GCS citation.                                    */
3234
    /* -------------------------------------------------------------------- */
3235
    if (nGCS == KvUserDefined || eVersion == GEOTIFF_VERSION_1_0)
3,161✔
3236
    {
3237
        const OGR_SRSNode *poGCS = poSRS->GetAttrNode("GEOGCS");
3,144✔
3238

3239
        if (poGCS != nullptr && poGCS->GetChild(0) != nullptr)
3,144✔
3240
        {
3241
            oMapAsciiKeys[GeogCitationGeoKey] = poGCS->GetChild(0)->GetValue();
3,027✔
3242
        }
3243
    }
3244

3245
    /* -------------------------------------------------------------------- */
3246
    /*      Try to identify the GCS/datum, scanning the EPSG datum file for */
3247
    /*      a match.                                                        */
3248
    /* -------------------------------------------------------------------- */
3249
    if (nPCS == KvUserDefined)
3,161✔
3250
    {
3251
        if (nGCS == KvUserDefined && poSRS->IsGeographic() &&
2,705✔
3252
            std::fabs(poSRS->GetAngularUnits() - CPLAtof(SRS_UA_DEGREE_CONV)) <
471✔
3253
                1e-9)
3254
        {
3255
            if (nDatum == Datum_North_American_Datum_1927)
464✔
3256
                nGCS = GCS_NAD27;
×
3257
            else if (nDatum == Datum_North_American_Datum_1983)
464✔
3258
                nGCS = GCS_NAD83;
×
3259
            else if (nDatum == Datum_WGS84 || nDatum == DatumE_WGS84)
464✔
3260
                nGCS = GCS_WGS_84;
459✔
3261
        }
3262

3263
        if (nGCS != KvUserDefined)
2,234✔
3264
        {
3265
            GTIFKeySet(psGTIF, GeographicTypeGeoKey, TYPE_SHORT, 1, nGCS);
2,051✔
3266
        }
3267
        else if (nDatum != KvUserDefined)
183✔
3268
        {
3269
            GTIFKeySet(psGTIF, GeographicTypeGeoKey, TYPE_SHORT, 1,
27✔
3270
                       KvUserDefined);
3271
            GTIFKeySet(psGTIF, GeogGeodeticDatumGeoKey, TYPE_SHORT, 1, nDatum);
27✔
3272
        }
3273
        else if (nSpheroid != KvUserDefined)
156✔
3274
        {
3275
            GTIFKeySet(psGTIF, GeographicTypeGeoKey, TYPE_SHORT, 1,
2✔
3276
                       KvUserDefined);
3277
            GTIFKeySet(psGTIF, GeogGeodeticDatumGeoKey, TYPE_SHORT, 1,
2✔
3278
                       KvUserDefined);
3279
            GTIFKeySet(psGTIF, GeogEllipsoidGeoKey, TYPE_SHORT, 1, nSpheroid);
2✔
3280
        }
3281
        else if (dfSemiMajor != 0.0)
154✔
3282
        {
3283
            GTIFKeySet(psGTIF, GeographicTypeGeoKey, TYPE_SHORT, 1,
40✔
3284
                       KvUserDefined);
3285
            GTIFKeySet(psGTIF, GeogGeodeticDatumGeoKey, TYPE_SHORT, 1,
40✔
3286
                       KvUserDefined);
3287
            GTIFKeySet(psGTIF, GeogEllipsoidGeoKey, TYPE_SHORT, 1,
40✔
3288
                       KvUserDefined);
3289
            GTIFKeySet(psGTIF, GeogSemiMajorAxisGeoKey, TYPE_DOUBLE, 1,
40✔
3290
                       dfSemiMajor);
3291
            if (dfInvFlattening == 0.0)
40✔
3292
                GTIFKeySet(psGTIF, GeogSemiMinorAxisGeoKey, TYPE_DOUBLE, 1,
13✔
3293
                           dfSemiMajor);
3294
            else
3295
                GTIFKeySet(psGTIF, GeogInvFlatteningGeoKey, TYPE_DOUBLE, 1,
27✔
3296
                           dfInvFlattening);
3297
        }
3298
        else if (poSRS->GetAttrValue("DATUM") != nullptr &&
114✔
3299
                 strstr(poSRS->GetAttrValue("DATUM"), "unknown") == nullptr &&
114✔
3300
                 strstr(poSRS->GetAttrValue("DATUM"), "unnamed") == nullptr)
×
3301

3302
        {
3303
            CPLError(CE_Warning, CPLE_AppDefined,
×
3304
                     "Couldn't translate `%s' to a GeoTIFF datum.",
3305
                     poSRS->GetAttrValue("DATUM"));
3306
        }
3307

3308
        if (nGCS == KvUserDefined || eVersion == GEOTIFF_VERSION_1_0)
2,234✔
3309
        {
3310
            // Always set InvFlattening if it is available.
3311
            // So that it doesn't need to calculate from SemiMinor.
3312
            if (dfInvFlattening != 0.0)
2,225✔
3313
                GTIFKeySet(psGTIF, GeogInvFlatteningGeoKey, TYPE_DOUBLE, 1,
2,096✔
3314
                           dfInvFlattening);
3315
            // Always set SemiMajor to keep the precision and in case of
3316
            // editing.
3317
            if (dfSemiMajor != 0.0)
2,225✔
3318
                GTIFKeySet(psGTIF, GeogSemiMajorAxisGeoKey, TYPE_DOUBLE, 1,
2,111✔
3319
                           dfSemiMajor);
3320

3321
            if (nGCS == KvUserDefined &&
2,408✔
3322
                CPLTestBool(CPLGetConfigOption("GTIFF_ESRI_CITATION", "YES")))
183✔
3323
            {
3324
                SetGeogCSCitation(psGTIF, oMapAsciiKeys, poSRS,
183✔
3325
                                  osAngUnitName.c_str(), nDatum, nSpheroid);
3326
            }
3327
        }
3328
    }
3329

3330
/* -------------------------------------------------------------------- */
3331
/*      Do we have TOWGS84 parameters?                                  */
3332
/* -------------------------------------------------------------------- */
3333
#if !defined(GEO_NORMALIZE_DISABLE_TOWGS84)
3334
    double adfTOWGS84[7] = {0.0};
3,161✔
3335

3336
    if ((nGCS == KvUserDefined || eVersion == GEOTIFF_VERSION_1_0) &&
6,305✔
3337
        poSRS->GetTOWGS84(adfTOWGS84) == OGRERR_NONE)
3,144✔
3338
    {
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
        bool bUseReferenceTOWGS84 = false;
23✔
3343
        const char *pszAuthName = poSRS->GetAuthorityName(nullptr);
23✔
3344
        const char *pszAuthCode = poSRS->GetAuthorityCode(nullptr);
23✔
3345
        if (pszAuthName && EQUAL(pszAuthName, "EPSG") && pszAuthCode)
23✔
3346
        {
3347
            CPLErrorHandlerPusher oErrorHandler(CPLQuietErrorHandler);
24✔
3348
            OGRSpatialReference oRefSRS;
24✔
3349
            double adfRefTOWGS84[7] = {0.0};
12✔
3350
            if (oRefSRS.importFromEPSG(atoi(pszAuthCode)) == OGRERR_NONE)
12✔
3351
            {
3352
                oRefSRS.AddGuessedTOWGS84();
12✔
3353
                if (oRefSRS.GetTOWGS84(adfRefTOWGS84) == OGRERR_NONE &&
15✔
3354
                    memcmp(adfRefTOWGS84, adfTOWGS84, sizeof(adfTOWGS84)) == 0)
3✔
3355
                {
3356
                    bUseReferenceTOWGS84 = true;
2✔
3357
                }
3358
            }
3359
        }
3360
        const char *pszWriteTOWGS84 =
3361
            CPLGetConfigOption("GTIFF_WRITE_TOWGS84", "AUTO");
23✔
3362
        if ((EQUAL(pszWriteTOWGS84, "YES") || EQUAL(pszWriteTOWGS84, "TRUE") ||
23✔
3363
             EQUAL(pszWriteTOWGS84, "ON")) ||
22✔
3364
            (!bUseReferenceTOWGS84 && EQUAL(pszWriteTOWGS84, "AUTO")))
22✔
3365
        {
3366
            if (adfTOWGS84[3] == 0.0 && adfTOWGS84[4] == 0.0 &&
22✔
3367
                adfTOWGS84[5] == 0.0 && adfTOWGS84[6] == 0.0)
18✔
3368
            {
3369
                if (nGCS == GCS_WGS_84 && adfTOWGS84[0] == 0.0 &&
15✔
3370
                    adfTOWGS84[1] == 0.0 && adfTOWGS84[2] == 0.0)
9✔
3371
                {
3372
                    ;  // Do nothing.
3373
                }
3374
                else
3375
                    GTIFKeySet(psGTIF, GeogTOWGS84GeoKey, TYPE_DOUBLE, 3,
6✔
3376
                               adfTOWGS84);
3377
            }
3378
            else
3379
                GTIFKeySet(psGTIF, GeogTOWGS84GeoKey, TYPE_DOUBLE, 7,
7✔
3380
                           adfTOWGS84);
3381
        }
3382
    }
3383
#endif
3384

3385
    /* -------------------------------------------------------------------- */
3386
    /*      Do we have vertical information to set?                         */
3387
    /* -------------------------------------------------------------------- */
3388
    if (poSRS->GetAttrValue("COMPD_CS|VERT_CS") != nullptr)
3,161✔
3389
    {
3390
        bool bGotVertCSCode = false;
13✔
3391
        const char *pszVertCSCode = poSRS->GetAuthorityCode("COMPD_CS|VERT_CS");
13✔
3392
        const char *pszVertCSAuthName =
3393
            poSRS->GetAuthorityName("COMPD_CS|VERT_CS");
13✔
3394
        if (pszVertCSCode && pszVertCSAuthName && atoi(pszVertCSCode) &&
13✔
3395
            EQUAL(pszVertCSAuthName, "EPSG"))
11✔
3396
        {
3397
            bGotVertCSCode = true;
10✔
3398
            GTIFKeySet(psGTIF, VerticalCSTypeGeoKey, TYPE_SHORT, 1,
10✔
3399
                       atoi(pszVertCSCode));
3400
        }
3401
        else if (eVersion >= GEOTIFF_VERSION_1_1)
3✔
3402
        {
3403
            GTIFKeySet(psGTIF, VerticalCSTypeGeoKey, TYPE_SHORT, 1,
3✔
3404
                       KvUserDefined);
3405
        }
3406

3407
        if (eVersion == GEOTIFF_VERSION_1_0 || !bGotVertCSCode)
13✔
3408
        {
3409
            oMapAsciiKeys[VerticalCitationGeoKey] =
×
3410
                poSRS->GetAttrValue("COMPD_CS|VERT_CS");
3✔
3411

3412
            const char *pszVertDatumCode =
3413
                poSRS->GetAuthorityCode("COMPD_CS|VERT_CS|VERT_DATUM");
3✔
3414
            const char *pszVertDatumAuthName =
3415
                poSRS->GetAuthorityName("COMPD_CS|VERT_CS|VERT_DATUM");
3✔
3416
            if (pszVertDatumCode && pszVertDatumAuthName &&
3✔
3417
                atoi(pszVertDatumCode) && EQUAL(pszVertDatumAuthName, "EPSG"))
2✔
3418
            {
3419
                GTIFKeySet(psGTIF, VerticalDatumGeoKey, TYPE_SHORT, 1,
1✔
3420
                           atoi(pszVertDatumCode));
3421
            }
3422
            else if (eVersion >= GEOTIFF_VERSION_1_1)
2✔
3423
            {
3424
                GTIFKeySet(psGTIF, VerticalDatumGeoKey, TYPE_SHORT, 1,
2✔
3425
                           KvUserDefined);
3426
            }
3427

3428
            const char *pszVertUnitCode =
3429
                poSRS->GetAuthorityCode("COMPD_CS|VERT_CS|UNIT");
3✔
3430
            const char *pszVertUnitAuthName =
3431
                poSRS->GetAuthorityName("COMPD_CS|VERT_CS|UNIT");
3✔
3432
            if (pszVertUnitCode && pszVertUnitAuthName &&
3✔
3433
                atoi(pszVertUnitCode) && EQUAL(pszVertUnitAuthName, "EPSG"))
3✔
3434
            {
3435
                GTIFKeySet(psGTIF, VerticalUnitsGeoKey, TYPE_SHORT, 1,
3✔
3436
                           atoi(pszVertUnitCode));
3437
            }
3438
            else if (eVersion >= GEOTIFF_VERSION_1_1)
×
3439
            {
3440
                GTIFKeySet(psGTIF, VerticalUnitsGeoKey, TYPE_SHORT, 1,
×
3441
                           KvUserDefined);
3442
            }
3443
        }
3444
    }
3445
    else if (eVersion >= GEOTIFF_VERSION_1_1 && nVerticalCSKeyValue != 0)
3,148✔
3446
    {
3447
        GTIFKeySet(psGTIF, VerticalCSTypeGeoKey, TYPE_SHORT, 1,
3✔
3448
                   nVerticalCSKeyValue);
3449
    }
3450

3451
    const double dfCoordinateEpoch = poSRS->GetCoordinateEpoch();
3,161✔
3452
    if (dfCoordinateEpoch > 0)
3,161✔
3453
    {
3454
        GTIFKeySet(psGTIF, CoordinateEpochGeoKey, TYPE_DOUBLE, 1,
2✔
3455
                   dfCoordinateEpoch);
3456
    }
3457

3458
    /* -------------------------------------------------------------------- */
3459
    /*      Write all ascii keys                                            */
3460
    /* -------------------------------------------------------------------- */
3461
    for (const auto &oIter : oMapAsciiKeys)
7,201✔
3462
    {
3463
        GTIFKeySet(psGTIF, oIter.first, TYPE_ASCII, 0, oIter.second.c_str());
4,040✔
3464
    }
3465

3466
    return TRUE;
6,322✔
3467
}
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
{
3477
    OGRSpatialReferenceH hSRS = nullptr;
×
3478
    if (ppszWKT)
×
3479
        *ppszWKT = nullptr;
×
3480
    CPLErr eErr =
3481
        GTIFWktFromMemBufEx(nSize, pabyBuffer, &hSRS, padfGeoTransform,
×
3482
                            pnGCPCount, ppasGCPList, nullptr, nullptr);
3483
    if (eErr == CE_None)
×
3484
    {
3485
        if (hSRS && ppszWKT)
×
3486
        {
3487
            OSRExportToWkt(hSRS, ppszWKT);
×
3488
        }
3489
    }
3490
    OSRDestroySpatialReference(hSRS);
×
3491
    return eErr;
×
3492
}
3493

3494
CPLErr GTIFWktFromMemBufEx(int nSize, unsigned char *pabyBuffer,
580✔
3495
                           OGRSpatialReferenceH *phSRS,
3496
                           double *padfGeoTransform, int *pnGCPCount,
3497
                           GDAL_GCP **ppasGCPList, int *pbPixelIsPoint,
3498
                           char ***ppapszRPCMD)
3499

3500
{
3501
    const std::string osFilename(
3502
        VSIMemGenerateHiddenFilename("wkt_from_mem_buf.tif"));
1,160✔
3503

3504
    /* -------------------------------------------------------------------- */
3505
    /*      Initialization of libtiff and libgeotiff.                       */
3506
    /* -------------------------------------------------------------------- */
3507
    GTiffOneTimeInit();  // For RPC tag.
580✔
3508
    LibgeotiffOneTimeInit();
580✔
3509

3510
    /* -------------------------------------------------------------------- */
3511
    /*      Create a memory file from the buffer.                           */
3512
    /* -------------------------------------------------------------------- */
3513
    VSILFILE *fp =
3514
        VSIFileFromMemBuffer(osFilename.c_str(), pabyBuffer, nSize, FALSE);
580✔
3515
    if (fp == nullptr)
580✔
3516
        return CE_Failure;
×
3517

3518
    /* -------------------------------------------------------------------- */
3519
    /*      Initialize access to the memory geotiff structure.              */
3520
    /* -------------------------------------------------------------------- */
3521
    TIFF *hTIFF = VSI_TIFFOpen(osFilename.c_str(), "rc", fp);
580✔
3522

3523
    if (hTIFF == nullptr)
580✔
3524
    {
3525
        CPLError(CE_Failure, CPLE_AppDefined,
×
3526
                 "TIFF/GeoTIFF structure is corrupt.");
3527
        VSIUnlink(osFilename.c_str());
×
3528
        CPL_IGNORE_RET_VAL(VSIFCloseL(fp));
×
3529
        return CE_Failure;
×
3530
    }
3531

3532
    /* -------------------------------------------------------------------- */
3533
    /*      Get the projection definition.                                  */
3534
    /* -------------------------------------------------------------------- */
3535
    bool bPixelIsPoint = false;
580✔
3536
    bool bPointGeoIgnore = false;
580✔
3537
    unsigned short nRasterType = 0;
580✔
3538

3539
    GTIF *hGTIF = GTIFNew(hTIFF);
580✔
3540
    if (hGTIF)
580✔
3541
        GTIFAttachPROJContext(hGTIF, OSRGetProjTLSContext());
580✔
3542

3543
    if (hGTIF != nullptr &&
1,160✔
3544
        GDALGTIFKeyGetSHORT(hGTIF, GTRasterTypeGeoKey, &nRasterType, 0, 1) ==
580✔
3545
            1 &&
1,160✔
3546
        nRasterType == static_cast<unsigned short>(RasterPixelIsPoint))
579✔
3547
    {
3548
        bPixelIsPoint = true;
5✔
3549
        bPointGeoIgnore =
3550
            CPLTestBool(CPLGetConfigOption("GTIFF_POINT_GEO_IGNORE", "FALSE"));
5✔
3551
    }
3552
    if (pbPixelIsPoint)
580✔
3553
        *pbPixelIsPoint = bPixelIsPoint;
580✔
3554
    if (ppapszRPCMD)
580✔
3555
        *ppapszRPCMD = nullptr;
580✔
3556

3557
    if (phSRS)
580✔
3558
    {
3559
        *phSRS = nullptr;
580✔
3560
        if (hGTIF != nullptr)
580✔
3561
        {
3562
            GTIFDefn *psGTIFDefn = GTIFAllocDefn();
580✔
3563
            if (GTIFGetDefn(hGTIF, psGTIFDefn))
580✔
3564
            {
3565
                *phSRS = GTIFGetOGISDefnAsOSR(hGTIF, psGTIFDefn);
579✔
3566
            }
3567
            GTIFFreeDefn(psGTIFDefn);
580✔
3568
        }
3569
    }
3570
    if (hGTIF)
580✔
3571
        GTIFFree(hGTIF);
580✔
3572

3573
    /* -------------------------------------------------------------------- */
3574
    /*      Get geotransform or tiepoints.                                  */
3575
    /* -------------------------------------------------------------------- */
3576
    double *padfTiePoints = nullptr;
580✔
3577
    double *padfScale = nullptr;
580✔
3578
    double *padfMatrix = nullptr;
580✔
3579
    int16_t nCount = 0;
580✔
3580

3581
    padfGeoTransform[0] = 0.0;
580✔
3582
    padfGeoTransform[1] = 1.0;
580✔
3583
    padfGeoTransform[2] = 0.0;
580✔
3584
    padfGeoTransform[3] = 0.0;
580✔
3585
    padfGeoTransform[4] = 0.0;
580✔
3586
    padfGeoTransform[5] = 1.0;
580✔
3587

3588
    *pnGCPCount = 0;
580✔
3589
    *ppasGCPList = nullptr;
580✔
3590

3591
    if (TIFFGetField(hTIFF, TIFFTAG_GEOPIXELSCALE, &nCount, &padfScale) &&
1,121✔
3592
        nCount >= 2)
541✔
3593
    {
3594
        padfGeoTransform[1] = padfScale[0];
541✔
3595
        padfGeoTransform[5] = -std::abs(padfScale[1]);
541✔
3596

3597
        if (TIFFGetField(hTIFF, TIFFTAG_GEOTIEPOINTS, &nCount,
541✔
3598
                         &padfTiePoints) &&
1,082✔
3599
            nCount >= 6)
541✔
3600
        {
3601
            padfGeoTransform[0] =
541✔
3602
                padfTiePoints[3] - padfTiePoints[0] * padfGeoTransform[1];
541✔
3603
            padfGeoTransform[3] =
541✔
3604
                padfTiePoints[4] - padfTiePoints[1] * padfGeoTransform[5];
541✔
3605

3606
            // Adjust for pixel is point in transform.
3607
            if (bPixelIsPoint && !bPointGeoIgnore)
541✔
3608
            {
3609
                padfGeoTransform[0] -=
4✔
3610
                    padfGeoTransform[1] * 0.5 + padfGeoTransform[2] * 0.5;
4✔
3611
                padfGeoTransform[3] -=
4✔
3612
                    padfGeoTransform[4] * 0.5 + padfGeoTransform[5] * 0.5;
4✔
3613
            }
3614
        }
3615
    }
3616
    else if (TIFFGetField(hTIFF, TIFFTAG_GEOTIEPOINTS, &nCount,
39✔
3617
                          &padfTiePoints) &&
52✔
3618
             nCount >= 6)
13✔
3619
    {
3620
        *pnGCPCount = nCount / 6;
13✔
3621
        *ppasGCPList =
13✔
3622
            static_cast<GDAL_GCP *>(CPLCalloc(sizeof(GDAL_GCP), *pnGCPCount));
13✔
3623

3624
        for (int iGCP = 0; iGCP < *pnGCPCount; iGCP++)
48✔
3625
        {
3626
            char szID[32] = {};
35✔
3627
            GDAL_GCP *psGCP = *ppasGCPList + iGCP;
35✔
3628

3629
            snprintf(szID, sizeof(szID), "%d", iGCP + 1);
35✔
3630
            psGCP->pszId = CPLStrdup(szID);
35✔
3631
            psGCP->pszInfo = CPLStrdup("");
35✔
3632
            psGCP->dfGCPPixel = padfTiePoints[iGCP * 6 + 0];
35✔
3633
            psGCP->dfGCPLine = padfTiePoints[iGCP * 6 + 1];
35✔
3634
            psGCP->dfGCPX = padfTiePoints[iGCP * 6 + 3];
35✔
3635
            psGCP->dfGCPY = padfTiePoints[iGCP * 6 + 4];
35✔
3636
            psGCP->dfGCPZ = padfTiePoints[iGCP * 6 + 5];
35✔
3637
        }
3638
    }
3639
    else if (TIFFGetField(hTIFF, TIFFTAG_GEOTRANSMATRIX, &nCount,
26✔
3640
                          &padfMatrix) &&
33✔
3641
             nCount == 16)
7✔
3642
    {
3643
        padfGeoTransform[0] = padfMatrix[3];
7✔
3644
        padfGeoTransform[1] = padfMatrix[0];
7✔
3645
        padfGeoTransform[2] = padfMatrix[1];
7✔
3646
        padfGeoTransform[3] = padfMatrix[7];
7✔
3647
        padfGeoTransform[4] = padfMatrix[4];
7✔
3648
        padfGeoTransform[5] = padfMatrix[5];
7✔
3649
    }
3650

3651
    /* -------------------------------------------------------------------- */
3652
    /*      Read RPC                                                        */
3653
    /* -------------------------------------------------------------------- */
3654
    if (ppapszRPCMD != nullptr)
580✔
3655
    {
3656
        *ppapszRPCMD = GTiffDatasetReadRPCTag(hTIFF);
580✔
3657
    }
3658

3659
    /* -------------------------------------------------------------------- */
3660
    /*      Cleanup.                                                        */
3661
    /* -------------------------------------------------------------------- */
3662
    XTIFFClose(hTIFF);
580✔
3663
    CPL_IGNORE_RET_VAL(VSIFCloseL(fp));
580✔
3664

3665
    VSIUnlink(osFilename.c_str());
580✔
3666

3667
    if (phSRS && *phSRS == nullptr)
580✔
3668
        return CE_Failure;
1✔
3669

3670
    return CE_None;
579✔
3671
}
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
{
3681
    OGRSpatialReference oSRS;
×
3682
    if (pszWKT != nullptr)
×
3683
        oSRS.importFromWkt(pszWKT);
×
3684
    return GTIFMemBufFromSRS(OGRSpatialReference::ToHandle(&oSRS),
×
3685
                             padfGeoTransform, nGCPCount, pasGCPList, pnSize,
3686
                             ppabyBuffer, FALSE, nullptr);
×
3687
}
3688

3689
CPLErr GTIFMemBufFromSRS(OGRSpatialReferenceH hSRS,
225✔
3690
                         const double *padfGeoTransform, int nGCPCount,
3691
                         const GDAL_GCP *pasGCPList, int *pnSize,
3692
                         unsigned char **ppabyBuffer, int bPixelIsPoint,
3693
                         char **papszRPCMD)
3694

3695
{
3696
    const std::string osFilename(
3697
        VSIMemGenerateHiddenFilename("wkt_from_mem_buf.tif"));
450✔
3698

3699
    /* -------------------------------------------------------------------- */
3700
    /*      Initialization of libtiff and libgeotiff.                       */
3701
    /* -------------------------------------------------------------------- */
3702
    GTiffOneTimeInit();  // For RPC tag.
225✔
3703
    LibgeotiffOneTimeInit();
225✔
3704

3705
    /* -------------------------------------------------------------------- */
3706
    /*      Initialize access to the memory geotiff structure.              */
3707
    /* -------------------------------------------------------------------- */
3708
    VSILFILE *fpL = VSIFOpenL(osFilename.c_str(), "w");
225✔
3709
    if (fpL == nullptr)
225✔
3710
        return CE_Failure;
×
3711

3712
    TIFF *hTIFF = VSI_TIFFOpen(osFilename.c_str(), "w", fpL);
225✔
3713

3714
    if (hTIFF == nullptr)
225✔
3715
    {
3716
        CPLError(CE_Failure, CPLE_AppDefined,
×
3717
                 "TIFF/GeoTIFF structure is corrupt.");
3718
        CPL_IGNORE_RET_VAL(VSIFCloseL(fpL));
×
3719
        VSIUnlink(osFilename.c_str());
×
3720
        return CE_Failure;
×
3721
    }
3722

3723
    /* -------------------------------------------------------------------- */
3724
    /*      Write some minimal set of image parameters.                     */
3725
    /* -------------------------------------------------------------------- */
3726
    TIFFSetField(hTIFF, TIFFTAG_IMAGEWIDTH, 1);
225✔
3727
    TIFFSetField(hTIFF, TIFFTAG_IMAGELENGTH, 1);
225✔
3728
    TIFFSetField(hTIFF, TIFFTAG_BITSPERSAMPLE, 8);
225✔
3729
    TIFFSetField(hTIFF, TIFFTAG_SAMPLESPERPIXEL, 1);
225✔
3730
    TIFFSetField(hTIFF, TIFFTAG_ROWSPERSTRIP, 1);
225✔
3731
    TIFFSetField(hTIFF, TIFFTAG_PLANARCONFIG, PLANARCONFIG_CONTIG);
225✔
3732
    TIFFSetField(hTIFF, TIFFTAG_PHOTOMETRIC, PHOTOMETRIC_MINISBLACK);
225✔
3733

3734
    /* -------------------------------------------------------------------- */
3735
    /*      Get the projection definition.                                  */
3736
    /* -------------------------------------------------------------------- */
3737

3738
    bool bPointGeoIgnore = false;
225✔
3739
    if (bPixelIsPoint)
225✔
3740
    {
3741
        bPointGeoIgnore =
3742
            CPLTestBool(CPLGetConfigOption("GTIFF_POINT_GEO_IGNORE", "FALSE"));
1✔
3743
    }
3744

3745
    GTIF *hGTIF = nullptr;
225✔
3746
    if (hSRS != nullptr || bPixelIsPoint)
225✔
3747
    {
3748
        hGTIF = GTIFNew(hTIFF);
225✔
3749
        if (hGTIF)
225✔
3750
        {
3751
            GTIFAttachPROJContext(hGTIF, OSRGetProjTLSContext());
225✔
3752

3753
            if (hSRS != nullptr)
225✔
3754
                GTIFSetFromOGISDefnEx(hGTIF, hSRS, GEOTIFF_KEYS_STANDARD,
225✔
3755
                                      GEOTIFF_VERSION_1_0);
3756

3757
            if (bPixelIsPoint)
225✔
3758
            {
3759
                GTIFKeySet(hGTIF, GTRasterTypeGeoKey, TYPE_SHORT, 1,
1✔
3760
                           RasterPixelIsPoint);
3761
            }
3762

3763
            GTIFWriteKeys(hGTIF);
225✔
3764
            GTIFFree(hGTIF);
225✔
3765
        }
3766
    }
3767

3768
    /* -------------------------------------------------------------------- */
3769
    /*      Set the geotransform, or GCPs.                                  */
3770
    /* -------------------------------------------------------------------- */
3771

3772
    if (padfGeoTransform[0] != 0.0 || padfGeoTransform[1] != 1.0 ||
28✔
3773
        padfGeoTransform[2] != 0.0 || padfGeoTransform[3] != 0.0 ||
28✔
3774
        padfGeoTransform[4] != 0.0 || std::abs(padfGeoTransform[5]) != 1.0)
253✔
3775
    {
3776

3777
        if (padfGeoTransform[2] == 0.0 && padfGeoTransform[4] == 0.0)
208✔
3778
        {
3779
            double adfPixelScale[3] = {padfGeoTransform[1],
205✔
3780
                                       fabs(padfGeoTransform[5]), 0.0};
205✔
3781

3782
            TIFFSetField(hTIFF, TIFFTAG_GEOPIXELSCALE, 3, adfPixelScale);
205✔
3783

3784
            double adfTiePoints[6] = {
205✔
3785
                0.0, 0.0, 0.0, padfGeoTransform[0], padfGeoTransform[3], 0.0};
205✔
3786

3787
            if (bPixelIsPoint && !bPointGeoIgnore)
205✔
3788
            {
3789
                adfTiePoints[3] +=
1✔
3790
                    padfGeoTransform[1] * 0.5 + padfGeoTransform[2] * 0.5;
1✔
3791
                adfTiePoints[4] +=
1✔
3792
                    padfGeoTransform[4] * 0.5 + padfGeoTransform[5] * 0.5;
1✔
3793
            }
3794

3795
            TIFFSetField(hTIFF, TIFFTAG_GEOTIEPOINTS, 6, adfTiePoints);
205✔
3796
        }
3797
        else
3798
        {
3799
            double adfMatrix[16] = {0.0};
3✔
3800

3801
            adfMatrix[0] = padfGeoTransform[1];
3✔
3802
            adfMatrix[1] = padfGeoTransform[2];
3✔
3803
            adfMatrix[3] = padfGeoTransform[0];
3✔
3804
            adfMatrix[4] = padfGeoTransform[4];
3✔
3805
            adfMatrix[5] = padfGeoTransform[5];
3✔
3806
            adfMatrix[7] = padfGeoTransform[3];
3✔
3807
            adfMatrix[15] = 1.0;
3✔
3808

3809
            if (bPixelIsPoint && !bPointGeoIgnore)
3✔
3810
            {
3811
                adfMatrix[3] +=
×
3812
                    padfGeoTransform[1] * 0.5 + padfGeoTransform[2] * 0.5;
×
3813
                adfMatrix[7] +=
×
3814
                    padfGeoTransform[4] * 0.5 + padfGeoTransform[5] * 0.5;
×
3815
            }
3816

3817
            TIFFSetField(hTIFF, TIFFTAG_GEOTRANSMATRIX, 16, adfMatrix);
3✔
3818
        }
3819
    }
3820

3821
    /* -------------------------------------------------------------------- */
3822
    /*      Otherwise write tiepoints if they are available.                */
3823
    /* -------------------------------------------------------------------- */
3824
    else if (nGCPCount > 0)
17✔
3825
    {
3826
        double *padfTiePoints =
3827
            static_cast<double *>(CPLMalloc(6 * sizeof(double) * nGCPCount));
5✔
3828

3829
        for (int iGCP = 0; iGCP < nGCPCount; iGCP++)
16✔
3830
        {
3831

3832
            padfTiePoints[iGCP * 6 + 0] = pasGCPList[iGCP].dfGCPPixel;
11✔
3833
            padfTiePoints[iGCP * 6 + 1] = pasGCPList[iGCP].dfGCPLine;
11✔
3834
            padfTiePoints[iGCP * 6 + 2] = 0;
11✔
3835
            padfTiePoints[iGCP * 6 + 3] = pasGCPList[iGCP].dfGCPX;
11✔
3836
            padfTiePoints[iGCP * 6 + 4] = pasGCPList[iGCP].dfGCPY;
11✔
3837
            padfTiePoints[iGCP * 6 + 5] = pasGCPList[iGCP].dfGCPZ;
11✔
3838
        }
3839

3840
        TIFFSetField(hTIFF, TIFFTAG_GEOTIEPOINTS, 6 * nGCPCount, padfTiePoints);
5✔
3841
        CPLFree(padfTiePoints);
5✔
3842
    }
3843

3844
    /* -------------------------------------------------------------------- */
3845
    /*      Write RPC                                                       */
3846
    /* -------------------------------------------------------------------- */
3847
    if (papszRPCMD != nullptr)
225✔
3848
    {
3849
        GTiffDatasetWriteRPCTag(hTIFF, papszRPCMD);
1✔
3850
    }
3851

3852
    /* -------------------------------------------------------------------- */
3853
    /*      Cleanup and return the created memory buffer.                   */
3854
    /* -------------------------------------------------------------------- */
3855
    GByte bySmallImage = 0;
225✔
3856

3857
    TIFFWriteEncodedStrip(hTIFF, 0, reinterpret_cast<char *>(&bySmallImage), 1);
225✔
3858
    TIFFWriteCheck(hTIFF, TIFFIsTiled(hTIFF), "GTIFMemBufFromWkt");
225✔
3859
    TIFFWriteDirectory(hTIFF);
225✔
3860

3861
    XTIFFClose(hTIFF);
225✔
3862
    CPL_IGNORE_RET_VAL(VSIFCloseL(fpL));
225✔
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
    GUIntBig nBigLength = 0;
225✔
3870

3871
    *ppabyBuffer = VSIGetMemFileBuffer(osFilename.c_str(), &nBigLength, TRUE);
225✔
3872
    *pnSize = static_cast<int>(nBigLength);
225✔
3873

3874
    return CE_None;
225✔
3875
}
STATUS · Troubleshooting · Open an Issue · Sales · Support · CAREERS · ENTERPRISE · START FREE · SCHEDULE DEMO
ANNOUNCEMENTS · TWITTER · TOS & SLA · Supported CI Services · What's a CI service? · Automated Testing

© 2025 Coveralls, Inc