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