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 10463 : void LibgeotiffOneTimeInit()
105 : {
106 10463 : std::lock_guard<std::mutex> oLock(oDeleteMutex);
107 :
108 : static bool bOneTimeInitDone = false;
109 :
110 10463 : if (bOneTimeInitDone)
111 9771 : return;
112 :
113 692 : bOneTimeInitDone = true;
114 :
115 : // This isn't thread-safe, so better do it now
116 692 : XTIFFInitialize();
117 : }
118 :
119 : /************************************************************************/
120 : /* GTIFToCPLRecyleString() */
121 : /* */
122 : /* This changes a string from the libgeotiff heap to the GDAL */
123 : /* heap. */
124 : /************************************************************************/
125 :
126 42578 : static void GTIFToCPLRecycleString(char **ppszTarget)
127 :
128 : {
129 42578 : if (*ppszTarget == nullptr)
130 54 : return;
131 :
132 42524 : char *pszTempString = CPLStrdup(*ppszTarget);
133 42524 : GTIFFreeMemory(*ppszTarget);
134 42524 : *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 258 : static PJ *GTIFGetOfficialCRS(GTIF *hGTIF, PJ_TYPE searchType,
460 : const char *pszAuthName, const char *pszName)
461 : {
462 258 : PJ *ret = nullptr;
463 : /* Search in database the corresponding EPSG 'official' name */
464 : auto ctx =
465 258 : static_cast<PJ_CONTEXT *>(GTIFGetPROJContext(hGTIF, true, nullptr));
466 : auto list =
467 258 : proj_create_from_name(ctx, pszAuthName, pszName, &searchType, 1,
468 : /* approximateMatch = */ false, 1, nullptr);
469 258 : if (list)
470 : {
471 258 : const auto listSize = proj_list_get_count(list);
472 258 : if (listSize == 1)
473 : {
474 29 : ret = proj_list_get(ctx, list, 0);
475 : }
476 258 : proj_list_destroy(list);
477 : }
478 258 : 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 8928 : OGRSpatialReferenceH GTIFGetOGISDefnAsOSR(GTIF *hGTIF, GTIFDefn *psDefn)
507 :
508 : {
509 17856 : OGRSpatialReference oSRS;
510 :
511 8928 : LibgeotiffOneTimeInit();
512 :
513 : #if LIBGEOTIFF_VERSION >= 1600
514 8928 : 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 8928 : 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 8928 : unsigned short verticalCSType = 0;
530 8928 : unsigned short verticalDatum = 0;
531 8928 : unsigned short verticalUnits = 0;
532 :
533 8928 : GDALGTIFKeyGetSHORT(hGTIF, VerticalCSTypeGeoKey, &verticalCSType, 0, 1);
534 8928 : GDALGTIFKeyGetSHORT(hGTIF, VerticalDatumGeoKey, &verticalDatum, 0, 1);
535 8928 : GDALGTIFKeyGetSHORT(hGTIF, VerticalUnitsGeoKey, &verticalUnits, 0, 1);
536 :
537 8928 : 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 8928 : if (psDefn->Model != ModelTypeProjected &&
598 4875 : 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 8584 : 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 8584 : 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 8584 : unsigned short bLinearUnitsMarkedCorrect = FALSE;
798 :
799 8584 : GDALGTIFKeyGetSHORT(hGTIF, ProjLinearUnitsInterpCorrectGeoKey,
800 : &bLinearUnitsMarkedCorrect, 0, 1);
801 :
802 8584 : 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 8584 : OGRBoolean linearUnitIsSet = FALSE;
837 8584 : if (psDefn->Model == ModelTypeProjected)
838 : {
839 4053 : char szCTString[512] = {'\0'};
840 4053 : if (psDefn->PCS != KvUserDefined)
841 : {
842 3954 : char *pszPCSName = nullptr;
843 :
844 : #if LIBGEOTIFF_VERSION >= 1600
845 3954 : GTIFGetPCSInfoEx(projContext,
846 : #else
847 : GTIFGetPCSInfo(
848 : #endif
849 3954 : psDefn->PCS, &pszPCSName, nullptr, nullptr,
850 : nullptr);
851 :
852 3954 : oSRS.SetProjCS(pszPCSName ? pszPCSName : "unnamed");
853 3954 : if (pszPCSName)
854 3954 : GTIFFreeMemory(pszPCSName);
855 :
856 3954 : 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 4053 : if (CheckCitationKeyForStatePlaneUTM(hGTIF, psDefn, &oSRS,
902 4053 : &linearUnitIsSet))
903 : {
904 0 : return OGRSpatialReference::ToHandle(oSRS.Clone());
905 : }
906 :
907 : /* Handle ESRI PE string in citation */
908 4053 : szCTString[0] = '\0';
909 4053 : if (GDALGTIFKeyGetASCII(hGTIF, GTCitationGeoKey, szCTString,
910 4053 : sizeof(szCTString)))
911 4029 : SetCitationToSRS(hGTIF, szCTString, sizeof(szCTString),
912 : GTCitationGeoKey, &oSRS, &linearUnitIsSet);
913 : }
914 :
915 : /* ==================================================================== */
916 : /* Setup the GeogCS */
917 : /* ==================================================================== */
918 8584 : char *pszGeogName = nullptr;
919 8584 : char *pszDatumName = nullptr;
920 8584 : char *pszPMName = nullptr;
921 8584 : char *pszSpheroidName = nullptr;
922 8584 : char *pszAngularUnits = nullptr;
923 8584 : char szGCSName[512] = {'\0'};
924 :
925 8584 : if (!
926 : #if LIBGEOTIFF_VERSION >= 1600
927 17168 : GTIFGetGCSInfoEx(projContext,
928 : #else
929 : GTIFGetGCSInfo(
930 : #endif
931 8584 : psDefn->GCS, &pszGeogName, nullptr, nullptr,
932 8689 : nullptr) &&
933 105 : GDALGTIFKeyGetASCII(hGTIF, GeogCitationGeoKey, szGCSName,
934 : sizeof(szGCSName)))
935 : {
936 92 : GetGeogCSFromCitation(szGCSName, sizeof(szGCSName), GeogCitationGeoKey,
937 : &pszGeogName, &pszDatumName, &pszPMName,
938 : &pszSpheroidName, &pszAngularUnits);
939 : }
940 : else
941 : {
942 8492 : GTIFToCPLRecycleString(&pszGeogName);
943 : }
944 :
945 8584 : 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 8584 : 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 8584 : 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 8584 : if (!pszDatumName)
982 : {
983 : #if LIBGEOTIFF_VERSION >= 1600
984 8518 : GTIFGetDatumInfoEx(projContext,
985 : #else
986 : GTIFGetDatumInfo(
987 : #endif
988 8518 : psDefn->Datum, &pszDatumName, nullptr);
989 8518 : GTIFToCPLRecycleString(&pszDatumName);
990 : }
991 :
992 8584 : double dfSemiMajor = 0.0;
993 8584 : double dfInvFlattening = 0.0;
994 8584 : if (!pszSpheroidName)
995 : {
996 : #if LIBGEOTIFF_VERSION >= 1600
997 8516 : GTIFGetEllipsoidInfoEx(projContext,
998 : #else
999 : GTIFGetEllipsoidInfo(
1000 : #endif
1001 8516 : psDefn->Ellipsoid, &pszSpheroidName, nullptr,
1002 : nullptr);
1003 8516 : 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 8584 : if (!pszPMName)
1019 : {
1020 : #if LIBGEOTIFF_VERSION >= 1600
1021 8494 : GTIFGetPMInfoEx(projContext,
1022 : #else
1023 : GTIFGetPMInfo(
1024 : #endif
1025 8494 : psDefn->PM, &pszPMName, nullptr);
1026 8494 : GTIFToCPLRecycleString(&pszPMName);
1027 : }
1028 : else
1029 : {
1030 90 : CPL_IGNORE_RET_VAL(
1031 90 : GDALGTIFKeyGetDOUBLE(hGTIF, GeogPrimeMeridianLongGeoKey,
1032 : &(psDefn->PMLongToGreenwich), 0, 1));
1033 : }
1034 :
1035 8584 : if (!pszAngularUnits)
1036 : {
1037 : #if LIBGEOTIFF_VERSION >= 1600
1038 8569 : GTIFGetUOMAngleInfoEx(projContext,
1039 : #else
1040 : GTIFGetUOMAngleInfo(
1041 : #endif
1042 8569 : psDefn->UOMAngle, &pszAngularUnits,
1043 : &psDefn->UOMAngleInDegrees);
1044 8569 : if (pszAngularUnits == nullptr)
1045 11 : pszAngularUnits = CPLStrdup("unknown");
1046 : else
1047 8558 : 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 8584 : 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 8584 : dfSemiMajor = psDefn->SemiMajor;
1068 8584 : 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 8573 : else if (dfInvFlattening == 0.0 &&
1076 8516 : ((psDefn->SemiMinor / psDefn->SemiMajor) < 0.99999999999999999 ||
1077 14 : (psDefn->SemiMinor / psDefn->SemiMajor) > 1.00000000000000001))
1078 : {
1079 8502 : dfInvFlattening =
1080 8502 : OSRCalcInvFlattening(psDefn->SemiMajor, psDefn->SemiMinor);
1081 :
1082 : /* Take official inverse flattening definition in the WGS84 case */
1083 8502 : if (fabs(dfSemiMajor - SRS_WGS84_SEMIMAJOR) < 1e-10 &&
1084 4909 : fabs(dfInvFlattening - SRS_WGS84_INVFLATTENING) < 1e-10)
1085 4832 : dfInvFlattening = SRS_WGS84_INVFLATTENING;
1086 : }
1087 8584 : if (!pszGeogName || strlen(pszGeogName) == 0)
1088 : {
1089 13 : CPLFree(pszGeogName);
1090 13 : pszGeogName = CPLStrdup(pszDatumName ? pszDatumName : "unknown");
1091 : }
1092 :
1093 8584 : oSRS.SetGeogCS(pszGeogName, pszDatumName, pszSpheroidName, dfSemiMajor,
1094 : dfInvFlattening, pszPMName, psDefn->PMLongToGreenwich,
1095 : pszAngularUnits,
1096 8584 : psDefn->UOMAngleInDegrees * CPLAtof(SRS_UA_DEGREE_CONV));
1097 :
1098 8584 : bool bGeog3DCRS = false;
1099 8584 : bool bSetDatumEllipsoidCode = true;
1100 8584 : bool bHasWarnedInconsistentGeogCRSEPSG = false;
1101 : {
1102 8584 : const int nGCS = psDefn->GCS;
1103 8584 : if (nGCS != KvUserDefined && nGCS > 0 &&
1104 8479 : psDefn->Model != ModelTypeGeocentric)
1105 : {
1106 16958 : OGRSpatialReference oSRSGeog;
1107 : const bool bGCSCodeValid =
1108 8479 : oSRSGeog.importFromEPSG(nGCS) == OGRERR_NONE;
1109 :
1110 : const std::string osGTiffSRSSource =
1111 16958 : CPLGetConfigOption("GTIFF_SRS_SOURCE", "");
1112 :
1113 : // GeoTIFF 1.0 might put a Geographic 3D code in GeodeticCRSGeoKey
1114 8479 : bool bTryCompareToEPSG = oSRSGeog.GetAxesCount() == 2;
1115 :
1116 8479 : if (psDefn->Datum != KvUserDefined)
1117 : {
1118 : char szCode[12];
1119 8479 : snprintf(szCode, sizeof(szCode), "%d", psDefn->Datum);
1120 : auto ctx = static_cast<PJ_CONTEXT *>(
1121 8479 : GTIFGetPROJContext(hGTIF, true, nullptr));
1122 8479 : auto datum = proj_create_from_database(
1123 : ctx, "EPSG", szCode, PJ_CATEGORY_DATUM, 0, nullptr);
1124 8479 : if (datum)
1125 : {
1126 8478 : 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 8478 : proj_destroy(datum);
1135 : }
1136 : }
1137 :
1138 8483 : 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 8479 : if (EQUAL(osGTiffSRSSource.c_str(), "EPSG"))
1160 : {
1161 1 : oSRS.CopyGeogCSFrom(&oSRSGeog);
1162 : }
1163 13294 : else if (osGTiffSRSSource.empty() && oSRSGeog.IsDynamic() &&
1164 4816 : psDefn->Model == ModelTypeGeographic)
1165 : {
1166 : // We should perhaps be more careful and detect overrides
1167 : // of geokeys...
1168 4485 : oSRS = oSRSGeog;
1169 4485 : bSetDatumEllipsoidCode = false;
1170 : }
1171 3993 : else if (bGCSCodeValid && osGTiffSRSSource.empty())
1172 : {
1173 3992 : oSRS.SetAuthority("GEOGCS", "EPSG", nGCS);
1174 : }
1175 : else
1176 : {
1177 1 : bSetDatumEllipsoidCode = false;
1178 : }
1179 :
1180 8479 : int nVertSRSCode = verticalCSType;
1181 8479 : 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 8479 : 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 8584 : if (bSetDatumEllipsoidCode)
1216 : {
1217 4096 : if (psDefn->Datum != KvUserDefined)
1218 4017 : oSRS.SetAuthority("DATUM", "EPSG", psDefn->Datum);
1219 :
1220 4096 : if (psDefn->Ellipsoid != KvUserDefined)
1221 4021 : oSRS.SetAuthority("SPHEROID", "EPSG", psDefn->Ellipsoid);
1222 : }
1223 :
1224 8584 : CPLFree(pszGeogName);
1225 8584 : CPLFree(pszDatumName);
1226 8584 : CPLFree(pszSpheroidName);
1227 8584 : CPLFree(pszPMName);
1228 8584 : CPLFree(pszAngularUnits);
1229 :
1230 : /* -------------------------------------------------------------------- */
1231 : /* Set projection units if not yet done */
1232 : /* -------------------------------------------------------------------- */
1233 8584 : if (psDefn->Model == ModelTypeProjected && !linearUnitIsSet)
1234 : {
1235 4046 : char *pszUnitsName = nullptr;
1236 :
1237 4046 : if (psDefn->UOMLength != KvUserDefined)
1238 : {
1239 : #if LIBGEOTIFF_VERSION >= 1600
1240 4039 : GTIFGetUOMLengthInfoEx(projContext,
1241 : #else
1242 : GTIFGetUOMLengthInfo(
1243 : #endif
1244 4039 : psDefn->UOMLength, &pszUnitsName, nullptr);
1245 : }
1246 :
1247 4046 : if (pszUnitsName != nullptr)
1248 : {
1249 : char szUOMLength[12];
1250 4039 : snprintf(szUOMLength, sizeof(szUOMLength), "%d", psDefn->UOMLength);
1251 4039 : oSRS.SetTargetLinearUnits(nullptr, pszUnitsName,
1252 : psDefn->UOMLengthInMeters, "EPSG",
1253 : szUOMLength);
1254 : }
1255 : else
1256 7 : oSRS.SetLinearUnits("unknown", psDefn->UOMLengthInMeters);
1257 :
1258 4046 : 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 8584 : unsigned short tmp = 0;
1268 8584 : bool bGotFromEPSG = false;
1269 4053 : if (psDefn->Model == ModelTypeProjected && psDefn->PCS != KvUserDefined &&
1270 3954 : GDALGTIFKeyGetSHORT(hGTIF, ProjectionGeoKey, &tmp, 0, 1) == 0 &&
1271 3953 : GDALGTIFKeyGetSHORT(hGTIF, ProjCoordTransGeoKey, &tmp, 0, 1) == 0 &&
1272 3953 : GDALGTIFKeyGetSHORT(hGTIF, GeographicTypeGeoKey, &tmp, 0, 1) == 0 &&
1273 3857 : GDALGTIFKeyGetSHORT(hGTIF, GeogGeodeticDatumGeoKey, &tmp, 0, 1) == 0 &&
1274 16492 : GDALGTIFKeyGetSHORT(hGTIF, GeogEllipsoidGeoKey, &tmp, 0, 1) == 0 &&
1275 3855 : CPLTestBool(CPLGetConfigOption("GTIFF_IMPORT_FROM_EPSG", "YES")))
1276 : {
1277 : // Save error state as importFromEPSGA() will call CPLReset()
1278 3855 : CPLErrorNum errNo = CPLGetLastErrorNo();
1279 3855 : CPLErr eErr = CPLGetLastErrorType();
1280 3855 : const char *pszTmp = CPLGetLastErrorMsg();
1281 3855 : char *pszLastErrorMsg = CPLStrdup(pszTmp ? pszTmp : "");
1282 3855 : CPLPushErrorHandler(CPLQuietErrorHandler);
1283 7710 : OGRSpatialReference oSRSTmp;
1284 3855 : OGRErr eImportErr = oSRSTmp.importFromEPSG(psDefn->PCS);
1285 3855 : CPLPopErrorHandler();
1286 : // Restore error state
1287 3855 : CPLErrorSetState(eErr, errNo, pszLastErrorMsg);
1288 3855 : CPLFree(pszLastErrorMsg);
1289 3855 : bGotFromEPSG = eImportErr == OGRERR_NONE;
1290 :
1291 3855 : if (bGotFromEPSG)
1292 : {
1293 : // See #6210. In case there's an overridden linear units, take it
1294 : // into account
1295 3855 : const char *pszUnitsName = nullptr;
1296 3855 : 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 3855 : if (fabs(dfUOMLengthInMeters - oSRSTmp.GetLinearUnits(nullptr)) >
1302 3855 : 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 3855 : if (bGeog3DCRS)
1319 : {
1320 1 : oSRSTmp.CopyGeogCSFrom(&oSRS);
1321 1 : oSRSTmp.UpdateCoordinateSystemFromGeogCRS();
1322 : }
1323 3855 : oSRS = std::move(oSRSTmp);
1324 : }
1325 : }
1326 :
1327 : #if !defined(GEO_NORMALIZE_DISABLE_TOWGS84)
1328 8607 : 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 8584 : 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 8584 : 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 8584 : if (psDefn->Model == ModelTypeProjected && psDefn->PCS != KvUserDefined &&
1600 3954 : !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 8584 : const bool bIs2DProjCRS = oSRS.IsProjected() && oSRS.GetAxesCount() == 2;
1648 8584 : if (bIs2DProjCRS)
1649 : {
1650 4050 : const char *pszProjCRSName = oSRS.GetAttrValue("PROJCS");
1651 4050 : if (pszProjCRSName)
1652 : {
1653 : // Hack to be able to read properly what we have written for
1654 : // ESRI:102113 (ESRI ancient WebMercator).
1655 4050 : if (EQUAL(pszProjCRSName, "WGS_1984_Web_Mercator"))
1656 2 : oSRS.SetFromUserInput("ESRI:102113");
1657 : // And for EPSG:900913
1658 4048 : else if (EQUAL(pszProjCRSName, "Google Maps Global Mercator"))
1659 0 : oSRS.importFromEPSG(900913);
1660 4048 : 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 4534 : if ((bIs2DProjCRS || oSRS.IsGeographic()) &&
1675 13246 : oSRS.GetAuthorityCode(nullptr) == nullptr &&
1676 128 : psDefn->UOMAngleInDegrees == 1.0)
1677 : {
1678 132 : const PJ_TYPE type = bIs2DProjCRS ? PJ_TYPE_PROJECTED_CRS
1679 18 : : oSRS.GetAxesCount() == 2
1680 18 : ? PJ_TYPE_GEOGRAPHIC_2D_CRS
1681 114 : : 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 114 : PJ *refCRS = GTIFGetOfficialCRS(hGTIF, type, "EPSG", oSRS.GetName());
1687 114 : if (!refCRS)
1688 93 : refCRS = GTIFGetOfficialCRS(hGTIF, type, nullptr, oSRS.GetName());
1689 114 : 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 8584 : bool bNeedManualVertCS = false;
1748 8584 : char citation[2048] = {'\0'};
1749 :
1750 : // See https://github.com/OSGeo/gdal/pull/4197
1751 8584 : if (verticalCSType > KvUserDefined || verticalDatum > KvUserDefined ||
1752 8584 : 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 8616 : 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 8584 : 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 8584 : 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 21 : oSRS.importFromWkt(pszWKT, apszOptionsImport);
1891 : }
1892 21 : CPLFree(pszWKT);
1893 : }
1894 :
1895 8584 : oSRS.StripTOWGS84IfKnownDatumAndAllowed();
1896 :
1897 8584 : double dfCoordinateEpoch = 0.0;
1898 8584 : if (GDALGTIFKeyGetDOUBLE(hGTIF, CoordinateEpochGeoKey, &dfCoordinateEpoch,
1899 8584 : 0, 1))
1900 : {
1901 2 : oSRS.SetCoordinateEpoch(dfCoordinateEpoch);
1902 : }
1903 :
1904 8584 : return OGRSpatialReference::ToHandle(oSRS.Clone());
1905 : }
1906 :
1907 : /************************************************************************/
1908 : /* GTIFGetOGISDefn() */
1909 : /************************************************************************/
1910 :
1911 0 : char *GTIFGetOGISDefn(GTIF *hGTIF, GTIFDefn *psDefn)
1912 : {
1913 0 : OGRSpatialReferenceH hSRS = GTIFGetOGISDefnAsOSR(hGTIF, psDefn);
1914 :
1915 0 : char *pszWKT = nullptr;
1916 0 : if (hSRS && OGRSpatialReference::FromHandle(hSRS)->exportToWkt(&pszWKT) ==
1917 : OGRERR_NONE)
1918 : {
1919 0 : OSRDestroySpatialReference(hSRS);
1920 0 : return pszWKT;
1921 : }
1922 0 : CPLFree(pszWKT);
1923 0 : OSRDestroySpatialReference(hSRS);
1924 :
1925 0 : return nullptr;
1926 : }
1927 :
1928 : /************************************************************************/
1929 : /* OGCDatumName2EPSGDatumCode() */
1930 : /************************************************************************/
1931 :
1932 578 : static int OGCDatumName2EPSGDatumCode(GTIF *psGTIF, const char *pszOGCName)
1933 :
1934 : {
1935 578 : int nReturn = KvUserDefined;
1936 :
1937 : /* -------------------------------------------------------------------- */
1938 : /* Do we know it as a built in? */
1939 : /* -------------------------------------------------------------------- */
1940 578 : if (EQUAL(pszOGCName, "NAD27") ||
1941 578 : EQUAL(pszOGCName, "North_American_Datum_1927"))
1942 0 : return Datum_North_American_Datum_1927;
1943 578 : else if (EQUAL(pszOGCName, "NAD83") ||
1944 578 : EQUAL(pszOGCName, "North_American_Datum_1983"))
1945 0 : return Datum_North_American_Datum_1983;
1946 578 : else if (EQUAL(pszOGCName, "WGS84") || EQUAL(pszOGCName, "WGS_1984") ||
1947 48 : EQUAL(pszOGCName, "WGS 84"))
1948 530 : return Datum_WGS84;
1949 48 : else if (EQUAL(pszOGCName, "WGS72") || EQUAL(pszOGCName, "WGS_1972"))
1950 0 : return Datum_WGS72;
1951 :
1952 : /* Search in database */
1953 : auto ctx =
1954 48 : static_cast<PJ_CONTEXT *>(GTIFGetPROJContext(psGTIF, true, nullptr));
1955 48 : const PJ_TYPE searchType = PJ_TYPE_GEODETIC_REFERENCE_FRAME;
1956 : auto list =
1957 48 : proj_create_from_name(ctx, "EPSG", pszOGCName, &searchType, 1,
1958 : /* approximateMatch = */ true, 10, nullptr);
1959 48 : if (list)
1960 : {
1961 48 : const auto listSize = proj_list_get_count(list);
1962 61 : for (int i = 0; nReturn == KvUserDefined && i < listSize; i++)
1963 : {
1964 13 : auto datum = proj_list_get(ctx, list, i);
1965 13 : if (datum)
1966 : {
1967 13 : const char *pszDatumName = proj_get_name(datum);
1968 13 : if (pszDatumName)
1969 : {
1970 13 : char *pszTmp = CPLStrdup(pszDatumName);
1971 13 : WKTMassageDatum(&pszTmp);
1972 13 : if (EQUAL(pszTmp, pszOGCName))
1973 : {
1974 3 : const char *pszCode = proj_get_id_code(datum, 0);
1975 3 : if (pszCode)
1976 : {
1977 3 : nReturn = atoi(pszCode);
1978 : }
1979 : }
1980 13 : CPLFree(pszTmp);
1981 : }
1982 : }
1983 13 : proj_destroy(datum);
1984 : }
1985 48 : proj_list_destroy(list);
1986 : }
1987 :
1988 48 : return nReturn;
1989 : }
1990 :
1991 : /************************************************************************/
1992 : /* GTIFSetFromOGISDefn() */
1993 : /* */
1994 : /* Write GeoTIFF projection tags from an OGC WKT definition. */
1995 : /************************************************************************/
1996 :
1997 0 : int GTIFSetFromOGISDefn(GTIF *psGTIF, const char *pszOGCWKT)
1998 :
1999 : {
2000 : /* -------------------------------------------------------------------- */
2001 : /* Create an OGRSpatialReference object corresponding to the */
2002 : /* string. */
2003 : /* -------------------------------------------------------------------- */
2004 :
2005 0 : OGRSpatialReference oSRS;
2006 0 : if (oSRS.importFromWkt(pszOGCWKT) != OGRERR_NONE)
2007 : {
2008 0 : return FALSE;
2009 : }
2010 0 : return GTIFSetFromOGISDefnEx(psGTIF, OGRSpatialReference::ToHandle(&oSRS),
2011 0 : GEOTIFF_KEYS_STANDARD, GEOTIFF_VERSION_1_0);
2012 : }
2013 :
2014 3388 : int GTIFSetFromOGISDefnEx(GTIF *psGTIF, OGRSpatialReferenceH hSRS,
2015 : GTIFFKeysFlavorEnum eFlavor,
2016 : GeoTIFFVersionEnum eVersion)
2017 : {
2018 6776 : std::map<geokey_t, std::string> oMapAsciiKeys;
2019 :
2020 3388 : GTIFKeySet(psGTIF, GTRasterTypeGeoKey, TYPE_SHORT, 1, RasterPixelIsArea);
2021 :
2022 3388 : const OGRSpatialReference *poSRS = OGRSpatialReference::FromHandle(hSRS);
2023 :
2024 : /* -------------------------------------------------------------------- */
2025 : /* Set version number. */
2026 : /* -------------------------------------------------------------------- */
2027 3388 : if (eVersion == GEOTIFF_VERSION_AUTO)
2028 : {
2029 6299 : if (poSRS->IsCompound() ||
2030 3143 : (poSRS->IsGeographic() && poSRS->GetAxesCount() == 3))
2031 : {
2032 15 : eVersion = GEOTIFF_VERSION_1_1;
2033 : }
2034 : else
2035 : {
2036 3141 : eVersion = GEOTIFF_VERSION_1_0;
2037 : }
2038 : }
2039 3388 : CPLAssert(eVersion == GEOTIFF_VERSION_1_0 ||
2040 : eVersion == GEOTIFF_VERSION_1_1);
2041 3388 : if (eVersion >= GEOTIFF_VERSION_1_1)
2042 : {
2043 : #if LIBGEOTIFF_VERSION >= 1600
2044 18 : GTIFSetVersionNumbers(psGTIF, GEOTIFF_SPEC_1_1_VERSION,
2045 : GEOTIFF_SPEC_1_1_KEY_REVISION,
2046 : GEOTIFF_SPEC_1_1_MINOR_REVISION);
2047 : #else
2048 : CPLError(CE_Warning, CPLE_AppDefined,
2049 : "Setting GeoTIFF 1.1 requires libgeotiff >= 1.6. Key values "
2050 : "will be written as GeoTIFF 1.1, but the version number "
2051 : "will be seen as 1.0, which might confuse GeoTIFF readers");
2052 : #endif
2053 : }
2054 :
2055 : /* -------------------------------------------------------------------- */
2056 : /* Get the ellipsoid definition. */
2057 : /* -------------------------------------------------------------------- */
2058 3388 : short nSpheroid = KvUserDefined;
2059 :
2060 4440 : if (poSRS->GetAuthorityName("PROJCS|GEOGCS|DATUM|SPHEROID") != nullptr &&
2061 1052 : EQUAL(poSRS->GetAuthorityName("PROJCS|GEOGCS|DATUM|SPHEROID"), "EPSG"))
2062 : {
2063 1052 : nSpheroid = static_cast<short>(
2064 1052 : atoi(poSRS->GetAuthorityCode("PROJCS|GEOGCS|DATUM|SPHEROID")));
2065 : }
2066 3979 : else if (poSRS->GetAuthorityName("GEOGCS|DATUM|SPHEROID") != nullptr &&
2067 1643 : EQUAL(poSRS->GetAuthorityName("GEOGCS|DATUM|SPHEROID"), "EPSG"))
2068 : {
2069 1643 : nSpheroid = static_cast<short>(
2070 1643 : atoi(poSRS->GetAuthorityCode("GEOGCS|DATUM|SPHEROID")));
2071 : }
2072 :
2073 3388 : OGRErr eErr = OGRERR_NONE;
2074 3388 : double dfSemiMajor = 0;
2075 3388 : double dfInvFlattening = 0;
2076 3388 : bool bHasEllipsoid = false;
2077 3388 : if (!poSRS->IsLocal())
2078 : {
2079 3386 : bHasEllipsoid = true;
2080 3386 : if (poSRS->IsCompound())
2081 : {
2082 26 : OGRSpatialReference oSRSTmp(*poSRS);
2083 13 : oSRSTmp.StripVertical();
2084 13 : bHasEllipsoid = CPL_TO_BOOL(!oSRSTmp.IsLocal());
2085 : }
2086 3386 : if (bHasEllipsoid)
2087 : {
2088 3385 : dfSemiMajor = poSRS->GetSemiMajor(&eErr);
2089 3385 : dfInvFlattening = poSRS->GetInvFlattening(&eErr);
2090 3385 : if (eErr != OGRERR_NONE)
2091 : {
2092 111 : dfSemiMajor = 0.0;
2093 111 : dfInvFlattening = 0.0;
2094 : }
2095 : }
2096 : }
2097 :
2098 : /* -------------------------------------------------------------------- */
2099 : /* Get the Datum so we can special case a few PCS codes. */
2100 : /* -------------------------------------------------------------------- */
2101 3388 : int nDatum = KvUserDefined;
2102 :
2103 4443 : if (poSRS->GetAuthorityName("PROJCS|GEOGCS|DATUM") != nullptr &&
2104 1055 : EQUAL(poSRS->GetAuthorityName("PROJCS|GEOGCS|DATUM"), "EPSG"))
2105 1055 : nDatum = atoi(poSRS->GetAuthorityCode("PROJCS|GEOGCS|DATUM"));
2106 3974 : else if (poSRS->GetAuthorityName("GEOGCS|DATUM") != nullptr &&
2107 1641 : EQUAL(poSRS->GetAuthorityName("GEOGCS|DATUM"), "EPSG"))
2108 1641 : nDatum = atoi(poSRS->GetAuthorityCode("GEOGCS|DATUM"));
2109 692 : else if (poSRS->GetAttrValue("DATUM") != nullptr)
2110 : nDatum =
2111 578 : OGCDatumName2EPSGDatumCode(psGTIF, poSRS->GetAttrValue("DATUM"));
2112 :
2113 : /* -------------------------------------------------------------------- */
2114 : /* Get the GCS if possible. */
2115 : /* -------------------------------------------------------------------- */
2116 3388 : int nGCS = KvUserDefined;
2117 :
2118 4417 : if (poSRS->GetAuthorityName("PROJCS|GEOGCS") != nullptr &&
2119 1029 : EQUAL(poSRS->GetAuthorityName("PROJCS|GEOGCS"), "EPSG"))
2120 1029 : nGCS = atoi(poSRS->GetAuthorityCode("PROJCS|GEOGCS"));
2121 4000 : else if (poSRS->GetAuthorityName("GEOGCS") != nullptr &&
2122 1641 : EQUAL(poSRS->GetAuthorityName("GEOGCS"), "EPSG"))
2123 1640 : nGCS = atoi(poSRS->GetAuthorityCode("GEOGCS"));
2124 :
2125 3388 : int nVerticalCSKeyValue = 0;
2126 5563 : bool hasEllipsoidHeight = !poSRS->IsCompound() && poSRS->IsGeographic() &&
2127 2175 : poSRS->GetAxesCount() == 3;
2128 3388 : if (nGCS == 4937 && eVersion >= GEOTIFF_VERSION_1_1)
2129 : {
2130 : // Workaround a bug of PROJ 6.3.0
2131 : // See https://github.com/OSGeo/PROJ/pull/1880
2132 : // EPSG:4937 = ETRS89 3D
2133 1 : hasEllipsoidHeight = true;
2134 1 : nVerticalCSKeyValue = nGCS;
2135 1 : nGCS = 4258; // ETRS89 2D
2136 : }
2137 3387 : else if (nGCS != KvUserDefined)
2138 : {
2139 5336 : OGRSpatialReference oGeogCRS;
2140 2668 : if (oGeogCRS.importFromEPSG(nGCS) == OGRERR_NONE &&
2141 2668 : oGeogCRS.IsGeographic() && oGeogCRS.GetAxesCount() == 3)
2142 : {
2143 3 : hasEllipsoidHeight = true;
2144 3 : if (eVersion >= GEOTIFF_VERSION_1_1)
2145 : {
2146 2 : const auto candidate_nVerticalCSKeyValue = nGCS;
2147 2 : nGCS = KvUserDefined;
2148 :
2149 : // In case of a geographic 3D CRS, find the corresponding
2150 : // geographic 2D CRS
2151 : auto ctx = static_cast<PJ_CONTEXT *>(
2152 2 : GTIFGetPROJContext(psGTIF, true, nullptr));
2153 2 : const auto type = PJ_TYPE_GEOGRAPHIC_2D_CRS;
2154 2 : auto list = proj_create_from_name(
2155 : ctx, "EPSG", oGeogCRS.GetName(), &type, 1,
2156 : /* approximateMatch = */ false,
2157 : 1, // result set limit size,
2158 : nullptr);
2159 2 : if (list && proj_list_get_count(list) == 1)
2160 : {
2161 2 : auto crs2D = proj_list_get(ctx, list, 0);
2162 2 : if (crs2D)
2163 : {
2164 2 : const char *pszCode = proj_get_id_code(crs2D, 0);
2165 2 : if (pszCode)
2166 : {
2167 2 : nVerticalCSKeyValue = candidate_nVerticalCSKeyValue;
2168 2 : nGCS = atoi(pszCode);
2169 : }
2170 2 : proj_destroy(crs2D);
2171 : }
2172 : }
2173 2 : proj_list_destroy(list);
2174 : }
2175 : }
2176 : }
2177 :
2178 : // Deprecated way of encoding ellipsoidal height
2179 3388 : if (hasEllipsoidHeight && nVerticalCSKeyValue == 0)
2180 : {
2181 1 : if (nGCS == 4979 || nDatum == 6326 || nSpheroid == 7030)
2182 : {
2183 0 : nVerticalCSKeyValue = 5030; // WGS_84_ellipsoid
2184 0 : if (nGCS == 4979 || nDatum == 6326)
2185 : {
2186 0 : nGCS = 4326;
2187 : }
2188 : }
2189 1 : else if (nDatum >= 6001 && nDatum <= 6033)
2190 : {
2191 0 : nVerticalCSKeyValue = nDatum - 1000;
2192 : }
2193 1 : else if (nSpheroid >= 7001 && nSpheroid <= 7033)
2194 : {
2195 1 : nVerticalCSKeyValue = nSpheroid - 2000;
2196 : }
2197 : }
2198 :
2199 3388 : if (nGCS > 32767)
2200 0 : nGCS = KvUserDefined;
2201 :
2202 : /* -------------------------------------------------------------------- */
2203 : /* Get the linear units. */
2204 : /* -------------------------------------------------------------------- */
2205 3388 : const char *pszLinearUOMNameTmp = nullptr;
2206 3388 : const double dfLinearUOM = poSRS->GetLinearUnits(&pszLinearUOMNameTmp);
2207 : const std::string osLinearUOMName(pszLinearUOMNameTmp ? pszLinearUOMNameTmp
2208 6776 : : "");
2209 3388 : int nUOMLengthCode = 9001; // Meters.
2210 :
2211 3388 : if (poSRS->GetAuthorityName("PROJCS|UNIT") != nullptr &&
2212 4471 : EQUAL(poSRS->GetAuthorityName("PROJCS|UNIT"), "EPSG") &&
2213 1083 : poSRS->GetAttrNode("PROJCS|UNIT") != poSRS->GetAttrNode("GEOGCS|UNIT"))
2214 1083 : nUOMLengthCode = atoi(poSRS->GetAuthorityCode("PROJCS|UNIT"));
2215 4610 : else if (EQUAL(osLinearUOMName.c_str(), SRS_UL_FOOT) ||
2216 2305 : fabs(dfLinearUOM - GTIFAtof(SRS_UL_FOOT_CONV)) < 0.0000001)
2217 0 : nUOMLengthCode = 9002; // International foot.
2218 4609 : else if (EQUAL(osLinearUOMName.c_str(), SRS_UL_US_FOOT) ||
2219 2304 : std::abs(dfLinearUOM - GTIFAtof(SRS_UL_US_FOOT_CONV)) < 0.0000001)
2220 2 : nUOMLengthCode = 9003; // US survey foot.
2221 2303 : else if (fabs(dfLinearUOM - 1.0) > 0.00000001)
2222 2 : nUOMLengthCode = KvUserDefined;
2223 :
2224 : /* -------------------------------------------------------------------- */
2225 : /* Get some authority values. */
2226 : /* -------------------------------------------------------------------- */
2227 3388 : int nPCS = KvUserDefined;
2228 :
2229 4411 : if (poSRS->GetAuthorityName("PROJCS") != nullptr &&
2230 1023 : EQUAL(poSRS->GetAuthorityName("PROJCS"), "EPSG"))
2231 : {
2232 1021 : nPCS = atoi(poSRS->GetAuthorityCode("PROJCS"));
2233 1021 : if (nPCS > 32767)
2234 0 : nPCS = KvUserDefined;
2235 : }
2236 :
2237 : /* -------------------------------------------------------------------- */
2238 : /* Handle the projection transformation. */
2239 : /* -------------------------------------------------------------------- */
2240 3388 : bool bWritePEString = false;
2241 3388 : bool bUnknownProjection = false;
2242 :
2243 3388 : const OGRSpatialReference *poSRSCompatibleOfWKT1 = poSRS;
2244 : #if PROJ_AT_LEAST_VERSION(6, 3, 0)
2245 6776 : OGRSpatialReference oSRSCompatibleOfWKT1;
2246 4473 : if (poSRS->IsProjected() && !poSRS->IsCompound() &&
2247 1085 : poSRS->GetAxesCount() == 3)
2248 : {
2249 0 : oSRSCompatibleOfWKT1 = *poSRS;
2250 0 : oSRSCompatibleOfWKT1.DemoteTo2D(nullptr);
2251 0 : poSRSCompatibleOfWKT1 = &oSRSCompatibleOfWKT1;
2252 0 : bWritePEString = true;
2253 : }
2254 : #endif
2255 :
2256 3388 : double dfAngUnitValue = 0;
2257 3388 : std::string osAngUnitName;
2258 3388 : if (bHasEllipsoid)
2259 : {
2260 3385 : const char *angUnitNameTmp = "";
2261 3385 : dfAngUnitValue = poSRS->GetAngularUnits(&angUnitNameTmp);
2262 3385 : osAngUnitName = angUnitNameTmp;
2263 : }
2264 :
2265 : // Convert angular projection parameters from its normalized value in degree
2266 : // to the units of GeogAngularUnitsGeoKey.
2267 : // Note: for GDAL <= 3.9.0, we always have written them in degrees !
2268 : // We can set GTIFF_WRITE_ANGULAR_PARAMS_IN_DEGREE=YES to get that
2269 : // non-conformant behavior...
2270 266 : const auto ConvertAngularParam = [dfAngUnitValue](double dfValInDeg)
2271 : {
2272 132 : constexpr double DEG_TO_RAD = M_PI / 180.0;
2273 132 : return dfAngUnitValue != 0 &&
2274 132 : std::fabs(dfAngUnitValue - DEG_TO_RAD) > 1e-10 &&
2275 4 : !CPLTestBool(CPLGetConfigOption(
2276 : "GTIFF_WRITE_ANGULAR_PARAMS_IN_DEGREE", "NO"))
2277 264 : ? dfValInDeg * DEG_TO_RAD / dfAngUnitValue
2278 132 : : dfValInDeg;
2279 3388 : };
2280 :
2281 : const char *pszProjection =
2282 3388 : poSRSCompatibleOfWKT1->GetAttrValue("PROJECTION");
2283 3388 : if (nPCS != KvUserDefined)
2284 : {
2285 : // If ESRI_PE flavor is explicitly required, then for EPSG:3857
2286 : // we will have to write a completely non-standard definition
2287 : // that requires not setting GTModelTypeGeoKey to ProjectedCSTypeGeoKey.
2288 1021 : if (eFlavor == GEOTIFF_KEYS_ESRI_PE && nPCS == 3857)
2289 : {
2290 1 : bWritePEString = true;
2291 : }
2292 : else
2293 : {
2294 1020 : GTIFKeySet(psGTIF, GTModelTypeGeoKey, TYPE_SHORT, 1,
2295 : ModelTypeProjected);
2296 1020 : GTIFKeySet(psGTIF, ProjectedCSTypeGeoKey, TYPE_SHORT, 1, nPCS);
2297 : }
2298 : }
2299 2367 : else if (poSRSCompatibleOfWKT1->IsGeocentric())
2300 : {
2301 2 : GTIFKeySet(psGTIF, GTModelTypeGeoKey, TYPE_SHORT, 1,
2302 : ModelTypeGeocentric);
2303 : }
2304 2365 : else if (pszProjection == nullptr)
2305 : {
2306 2294 : if (poSRSCompatibleOfWKT1->IsGeographic())
2307 2180 : GTIFKeySet(psGTIF, GTModelTypeGeoKey, TYPE_SHORT, 1,
2308 : ModelTypeGeographic);
2309 : // Otherwise, presumably something like LOCAL_CS.
2310 : }
2311 71 : else if (EQUAL(pszProjection, SRS_PT_ALBERS_CONIC_EQUAL_AREA))
2312 : {
2313 3 : GTIFKeySet(psGTIF, GTModelTypeGeoKey, TYPE_SHORT, 1,
2314 : ModelTypeProjected);
2315 3 : GTIFKeySet(psGTIF, ProjectedCSTypeGeoKey, TYPE_SHORT, 1, KvUserDefined);
2316 3 : GTIFKeySet(psGTIF, ProjectionGeoKey, TYPE_SHORT, 1, KvUserDefined);
2317 :
2318 3 : GTIFKeySet(psGTIF, ProjCoordTransGeoKey, TYPE_SHORT, 1,
2319 : CT_AlbersEqualArea);
2320 :
2321 3 : GTIFKeySet(psGTIF, ProjStdParallelGeoKey, TYPE_DOUBLE, 1,
2322 : ConvertAngularParam(poSRSCompatibleOfWKT1->GetNormProjParm(
2323 : SRS_PP_STANDARD_PARALLEL_1, 0.0)));
2324 :
2325 3 : GTIFKeySet(psGTIF, ProjStdParallel2GeoKey, TYPE_DOUBLE, 1,
2326 : ConvertAngularParam(poSRSCompatibleOfWKT1->GetNormProjParm(
2327 : SRS_PP_STANDARD_PARALLEL_2, 0.0)));
2328 :
2329 3 : GTIFKeySet(psGTIF, ProjNatOriginLatGeoKey, TYPE_DOUBLE, 1,
2330 : ConvertAngularParam(poSRSCompatibleOfWKT1->GetNormProjParm(
2331 : SRS_PP_LATITUDE_OF_CENTER, 0.0)));
2332 :
2333 3 : GTIFKeySet(psGTIF, ProjNatOriginLongGeoKey, TYPE_DOUBLE, 1,
2334 : ConvertAngularParam(poSRSCompatibleOfWKT1->GetNormProjParm(
2335 : SRS_PP_LONGITUDE_OF_CENTER, 0.0)));
2336 :
2337 3 : GTIFKeySet(
2338 : psGTIF, ProjFalseEastingGeoKey, TYPE_DOUBLE, 1,
2339 : poSRSCompatibleOfWKT1->GetProjParm(SRS_PP_FALSE_EASTING, 0.0));
2340 :
2341 3 : GTIFKeySet(
2342 : psGTIF, ProjFalseNorthingGeoKey, TYPE_DOUBLE, 1,
2343 : poSRSCompatibleOfWKT1->GetProjParm(SRS_PP_FALSE_NORTHING, 0.0));
2344 : }
2345 :
2346 68 : else if (poSRSCompatibleOfWKT1->GetUTMZone() != 0)
2347 : {
2348 10 : GTIFKeySet(psGTIF, GTModelTypeGeoKey, TYPE_SHORT, 1,
2349 : ModelTypeProjected);
2350 :
2351 10 : int bNorth = 0;
2352 10 : const int nZone = poSRSCompatibleOfWKT1->GetUTMZone(&bNorth);
2353 :
2354 10 : if (nDatum == Datum_North_American_Datum_1983 && nZone >= 3 &&
2355 1 : nZone <= 22 && bNorth && nUOMLengthCode == 9001)
2356 : {
2357 1 : nPCS = 26900 + nZone;
2358 :
2359 1 : GTIFKeySet(psGTIF, ProjectedCSTypeGeoKey, TYPE_SHORT, 1, nPCS);
2360 : }
2361 9 : else if (nDatum == Datum_North_American_Datum_1927 && nZone >= 3 &&
2362 1 : nZone <= 22 && bNorth && nUOMLengthCode == 9001)
2363 : {
2364 1 : nPCS = 26700 + nZone;
2365 :
2366 1 : GTIFKeySet(psGTIF, ProjectedCSTypeGeoKey, TYPE_SHORT, 1, nPCS);
2367 : }
2368 8 : else if (nDatum == Datum_WGS84 && nUOMLengthCode == 9001)
2369 : {
2370 2 : if (bNorth)
2371 1 : nPCS = 32600 + nZone;
2372 : else
2373 1 : nPCS = 32700 + nZone;
2374 :
2375 2 : GTIFKeySet(psGTIF, ProjectedCSTypeGeoKey, TYPE_SHORT, 1, nPCS);
2376 : }
2377 : else
2378 : {
2379 6 : const int nProjection = nZone + (bNorth ? 16000 : 16100);
2380 6 : GTIFKeySet(psGTIF, ProjectedCSTypeGeoKey, TYPE_SHORT, 1,
2381 : KvUserDefined);
2382 :
2383 6 : GTIFKeySet(psGTIF, ProjectionGeoKey, TYPE_SHORT, 1, nProjection);
2384 : }
2385 : }
2386 :
2387 58 : else if (EQUAL(pszProjection, SRS_PT_TRANSVERSE_MERCATOR))
2388 : {
2389 9 : GTIFKeySet(psGTIF, GTModelTypeGeoKey, TYPE_SHORT, 1,
2390 : ModelTypeProjected);
2391 9 : GTIFKeySet(psGTIF, ProjectedCSTypeGeoKey, TYPE_SHORT, 1, KvUserDefined);
2392 9 : GTIFKeySet(psGTIF, ProjectionGeoKey, TYPE_SHORT, 1, KvUserDefined);
2393 :
2394 9 : GTIFKeySet(psGTIF, ProjCoordTransGeoKey, TYPE_SHORT, 1,
2395 : CT_TransverseMercator);
2396 :
2397 9 : GTIFKeySet(psGTIF, ProjNatOriginLatGeoKey, TYPE_DOUBLE, 1,
2398 : ConvertAngularParam(poSRSCompatibleOfWKT1->GetNormProjParm(
2399 : SRS_PP_LATITUDE_OF_ORIGIN, 0.0)));
2400 :
2401 9 : GTIFKeySet(psGTIF, ProjNatOriginLongGeoKey, TYPE_DOUBLE, 1,
2402 : ConvertAngularParam(poSRSCompatibleOfWKT1->GetNormProjParm(
2403 : SRS_PP_CENTRAL_MERIDIAN, 0.0)));
2404 :
2405 9 : GTIFKeySet(
2406 : psGTIF, ProjScaleAtNatOriginGeoKey, TYPE_DOUBLE, 1,
2407 : poSRSCompatibleOfWKT1->GetNormProjParm(SRS_PP_SCALE_FACTOR, 1.0));
2408 :
2409 9 : GTIFKeySet(
2410 : psGTIF, ProjFalseEastingGeoKey, TYPE_DOUBLE, 1,
2411 : poSRSCompatibleOfWKT1->GetProjParm(SRS_PP_FALSE_EASTING, 0.0));
2412 :
2413 9 : GTIFKeySet(
2414 : psGTIF, ProjFalseNorthingGeoKey, TYPE_DOUBLE, 1,
2415 : poSRSCompatibleOfWKT1->GetProjParm(SRS_PP_FALSE_NORTHING, 0.0));
2416 : }
2417 :
2418 49 : else if (EQUAL(pszProjection, SRS_PT_TRANSVERSE_MERCATOR_SOUTH_ORIENTED))
2419 : {
2420 0 : GTIFKeySet(psGTIF, GTModelTypeGeoKey, TYPE_SHORT, 1,
2421 : ModelTypeProjected);
2422 0 : GTIFKeySet(psGTIF, ProjectedCSTypeGeoKey, TYPE_SHORT, 1, KvUserDefined);
2423 0 : GTIFKeySet(psGTIF, ProjectionGeoKey, TYPE_SHORT, 1, KvUserDefined);
2424 :
2425 0 : GTIFKeySet(psGTIF, ProjCoordTransGeoKey, TYPE_SHORT, 1,
2426 : CT_TransvMercator_SouthOriented);
2427 :
2428 0 : GTIFKeySet(psGTIF, ProjNatOriginLatGeoKey, TYPE_DOUBLE, 1,
2429 : ConvertAngularParam(poSRSCompatibleOfWKT1->GetNormProjParm(
2430 : SRS_PP_LATITUDE_OF_ORIGIN, 0.0)));
2431 :
2432 0 : GTIFKeySet(psGTIF, ProjNatOriginLongGeoKey, TYPE_DOUBLE, 1,
2433 : ConvertAngularParam(poSRSCompatibleOfWKT1->GetNormProjParm(
2434 : SRS_PP_CENTRAL_MERIDIAN, 0.0)));
2435 :
2436 0 : GTIFKeySet(
2437 : psGTIF, ProjScaleAtNatOriginGeoKey, TYPE_DOUBLE, 1,
2438 : poSRSCompatibleOfWKT1->GetNormProjParm(SRS_PP_SCALE_FACTOR, 1.0));
2439 :
2440 0 : GTIFKeySet(
2441 : psGTIF, ProjFalseEastingGeoKey, TYPE_DOUBLE, 1,
2442 : poSRSCompatibleOfWKT1->GetProjParm(SRS_PP_FALSE_EASTING, 0.0));
2443 :
2444 0 : GTIFKeySet(
2445 : psGTIF, ProjFalseNorthingGeoKey, TYPE_DOUBLE, 1,
2446 : poSRSCompatibleOfWKT1->GetProjParm(SRS_PP_FALSE_NORTHING, 0.0));
2447 : }
2448 :
2449 49 : else if (EQUAL(pszProjection, SRS_PT_MERCATOR_2SP) ||
2450 48 : EQUAL(pszProjection, SRS_PT_MERCATOR_1SP))
2451 :
2452 : {
2453 5 : GTIFKeySet(psGTIF, GTModelTypeGeoKey, TYPE_SHORT, 1,
2454 : ModelTypeProjected);
2455 5 : GTIFKeySet(psGTIF, ProjectedCSTypeGeoKey, TYPE_SHORT, 1, KvUserDefined);
2456 5 : GTIFKeySet(psGTIF, ProjectionGeoKey, TYPE_SHORT, 1, KvUserDefined);
2457 :
2458 5 : GTIFKeySet(psGTIF, ProjCoordTransGeoKey, TYPE_SHORT, 1, CT_Mercator);
2459 :
2460 5 : GTIFKeySet(psGTIF, ProjNatOriginLatGeoKey, TYPE_DOUBLE, 1,
2461 : ConvertAngularParam(poSRSCompatibleOfWKT1->GetNormProjParm(
2462 : SRS_PP_LATITUDE_OF_ORIGIN, 0.0)));
2463 :
2464 5 : GTIFKeySet(psGTIF, ProjNatOriginLongGeoKey, TYPE_DOUBLE, 1,
2465 : ConvertAngularParam(poSRSCompatibleOfWKT1->GetNormProjParm(
2466 : SRS_PP_CENTRAL_MERIDIAN, 0.0)));
2467 :
2468 5 : if (EQUAL(pszProjection, SRS_PT_MERCATOR_2SP))
2469 1 : GTIFKeySet(
2470 : psGTIF, ProjStdParallel1GeoKey, TYPE_DOUBLE, 1,
2471 : ConvertAngularParam(poSRSCompatibleOfWKT1->GetNormProjParm(
2472 : SRS_PP_STANDARD_PARALLEL_1, 0.0)));
2473 : else
2474 4 : GTIFKeySet(psGTIF, ProjScaleAtNatOriginGeoKey, TYPE_DOUBLE, 1,
2475 : poSRSCompatibleOfWKT1->GetNormProjParm(
2476 : SRS_PP_SCALE_FACTOR, 1.0));
2477 :
2478 5 : GTIFKeySet(
2479 : psGTIF, ProjFalseEastingGeoKey, TYPE_DOUBLE, 1,
2480 : poSRSCompatibleOfWKT1->GetProjParm(SRS_PP_FALSE_EASTING, 0.0));
2481 :
2482 5 : GTIFKeySet(
2483 : psGTIF, ProjFalseNorthingGeoKey, TYPE_DOUBLE, 1,
2484 : poSRSCompatibleOfWKT1->GetProjParm(SRS_PP_FALSE_NORTHING, 0.0));
2485 : }
2486 :
2487 44 : else if (EQUAL(pszProjection, SRS_PT_OBLIQUE_STEREOGRAPHIC))
2488 : {
2489 1 : GTIFKeySet(psGTIF, GTModelTypeGeoKey, TYPE_SHORT, 1,
2490 : ModelTypeProjected);
2491 1 : GTIFKeySet(psGTIF, ProjectedCSTypeGeoKey, TYPE_SHORT, 1, KvUserDefined);
2492 1 : GTIFKeySet(psGTIF, ProjectionGeoKey, TYPE_SHORT, 1, KvUserDefined);
2493 :
2494 1 : GTIFKeySet(psGTIF, ProjCoordTransGeoKey, TYPE_SHORT, 1,
2495 : CT_ObliqueStereographic);
2496 :
2497 1 : GTIFKeySet(psGTIF, ProjNatOriginLatGeoKey, TYPE_DOUBLE, 1,
2498 : ConvertAngularParam(poSRSCompatibleOfWKT1->GetNormProjParm(
2499 : SRS_PP_LATITUDE_OF_ORIGIN, 0.0)));
2500 :
2501 1 : GTIFKeySet(psGTIF, ProjNatOriginLongGeoKey, TYPE_DOUBLE, 1,
2502 : ConvertAngularParam(poSRSCompatibleOfWKT1->GetNormProjParm(
2503 : SRS_PP_CENTRAL_MERIDIAN, 0.0)));
2504 :
2505 1 : GTIFKeySet(
2506 : psGTIF, ProjScaleAtNatOriginGeoKey, TYPE_DOUBLE, 1,
2507 : poSRSCompatibleOfWKT1->GetNormProjParm(SRS_PP_SCALE_FACTOR, 1.0));
2508 :
2509 1 : GTIFKeySet(
2510 : psGTIF, ProjFalseEastingGeoKey, TYPE_DOUBLE, 1,
2511 : poSRSCompatibleOfWKT1->GetProjParm(SRS_PP_FALSE_EASTING, 0.0));
2512 :
2513 1 : GTIFKeySet(
2514 : psGTIF, ProjFalseNorthingGeoKey, TYPE_DOUBLE, 1,
2515 : poSRSCompatibleOfWKT1->GetProjParm(SRS_PP_FALSE_NORTHING, 0.0));
2516 : }
2517 :
2518 43 : else if (EQUAL(pszProjection, SRS_PT_STEREOGRAPHIC))
2519 : {
2520 0 : GTIFKeySet(psGTIF, GTModelTypeGeoKey, TYPE_SHORT, 1,
2521 : ModelTypeProjected);
2522 0 : GTIFKeySet(psGTIF, ProjectedCSTypeGeoKey, TYPE_SHORT, 1, KvUserDefined);
2523 0 : GTIFKeySet(psGTIF, ProjectionGeoKey, TYPE_SHORT, 1, KvUserDefined);
2524 :
2525 0 : GTIFKeySet(psGTIF, ProjCoordTransGeoKey, TYPE_SHORT, 1,
2526 : CT_Stereographic);
2527 :
2528 0 : GTIFKeySet(psGTIF, ProjCenterLatGeoKey, TYPE_DOUBLE, 1,
2529 : ConvertAngularParam(poSRSCompatibleOfWKT1->GetNormProjParm(
2530 : SRS_PP_LATITUDE_OF_ORIGIN, 0.0)));
2531 :
2532 0 : GTIFKeySet(psGTIF, ProjCenterLongGeoKey, TYPE_DOUBLE, 1,
2533 : ConvertAngularParam(poSRSCompatibleOfWKT1->GetNormProjParm(
2534 : SRS_PP_CENTRAL_MERIDIAN, 0.0)));
2535 :
2536 0 : GTIFKeySet(
2537 : psGTIF, ProjScaleAtNatOriginGeoKey, TYPE_DOUBLE, 1,
2538 : poSRSCompatibleOfWKT1->GetNormProjParm(SRS_PP_SCALE_FACTOR, 1.0));
2539 :
2540 0 : GTIFKeySet(
2541 : psGTIF, ProjFalseEastingGeoKey, TYPE_DOUBLE, 1,
2542 : poSRSCompatibleOfWKT1->GetProjParm(SRS_PP_FALSE_EASTING, 0.0));
2543 :
2544 0 : GTIFKeySet(
2545 : psGTIF, ProjFalseNorthingGeoKey, TYPE_DOUBLE, 1,
2546 : poSRSCompatibleOfWKT1->GetProjParm(SRS_PP_FALSE_NORTHING, 0.0));
2547 : }
2548 :
2549 43 : else if (EQUAL(pszProjection, SRS_PT_POLAR_STEREOGRAPHIC))
2550 : {
2551 4 : GTIFKeySet(psGTIF, GTModelTypeGeoKey, TYPE_SHORT, 1,
2552 : ModelTypeProjected);
2553 4 : GTIFKeySet(psGTIF, ProjectedCSTypeGeoKey, TYPE_SHORT, 1, KvUserDefined);
2554 4 : GTIFKeySet(psGTIF, ProjectionGeoKey, TYPE_SHORT, 1, KvUserDefined);
2555 :
2556 4 : GTIFKeySet(psGTIF, ProjCoordTransGeoKey, TYPE_SHORT, 1,
2557 : CT_PolarStereographic);
2558 :
2559 4 : GTIFKeySet(psGTIF, ProjNatOriginLatGeoKey, TYPE_DOUBLE, 1,
2560 : ConvertAngularParam(poSRSCompatibleOfWKT1->GetNormProjParm(
2561 : SRS_PP_LATITUDE_OF_ORIGIN, 0.0)));
2562 :
2563 4 : GTIFKeySet(psGTIF, ProjStraightVertPoleLongGeoKey, TYPE_DOUBLE, 1,
2564 : ConvertAngularParam(poSRSCompatibleOfWKT1->GetNormProjParm(
2565 : SRS_PP_CENTRAL_MERIDIAN, 0.0)));
2566 :
2567 4 : GTIFKeySet(
2568 : psGTIF, ProjScaleAtNatOriginGeoKey, TYPE_DOUBLE, 1,
2569 : poSRSCompatibleOfWKT1->GetNormProjParm(SRS_PP_SCALE_FACTOR, 1.0));
2570 :
2571 4 : GTIFKeySet(
2572 : psGTIF, ProjFalseEastingGeoKey, TYPE_DOUBLE, 1,
2573 : poSRSCompatibleOfWKT1->GetProjParm(SRS_PP_FALSE_EASTING, 0.0));
2574 :
2575 4 : GTIFKeySet(
2576 : psGTIF, ProjFalseNorthingGeoKey, TYPE_DOUBLE, 1,
2577 : poSRSCompatibleOfWKT1->GetProjParm(SRS_PP_FALSE_NORTHING, 0.0));
2578 : }
2579 :
2580 39 : else if (EQUAL(pszProjection, SRS_PT_HOTINE_OBLIQUE_MERCATOR))
2581 : {
2582 1 : GTIFKeySet(psGTIF, GTModelTypeGeoKey, TYPE_SHORT, 1,
2583 : ModelTypeProjected);
2584 1 : GTIFKeySet(psGTIF, ProjectedCSTypeGeoKey, TYPE_SHORT, 1, KvUserDefined);
2585 1 : GTIFKeySet(psGTIF, ProjectionGeoKey, TYPE_SHORT, 1, KvUserDefined);
2586 :
2587 1 : GTIFKeySet(psGTIF, ProjCoordTransGeoKey, TYPE_SHORT, 1,
2588 : CT_ObliqueMercator);
2589 :
2590 1 : GTIFKeySet(psGTIF, ProjCenterLatGeoKey, TYPE_DOUBLE, 1,
2591 : ConvertAngularParam(poSRSCompatibleOfWKT1->GetNormProjParm(
2592 : SRS_PP_LATITUDE_OF_CENTER, 0.0)));
2593 :
2594 1 : GTIFKeySet(psGTIF, ProjCenterLongGeoKey, TYPE_DOUBLE, 1,
2595 : ConvertAngularParam(poSRSCompatibleOfWKT1->GetNormProjParm(
2596 : SRS_PP_LONGITUDE_OF_CENTER, 0.0)));
2597 :
2598 1 : GTIFKeySet(psGTIF, ProjAzimuthAngleGeoKey, TYPE_DOUBLE, 1,
2599 : poSRSCompatibleOfWKT1->GetNormProjParm(SRS_PP_AZIMUTH, 0.0));
2600 :
2601 1 : GTIFKeySet(psGTIF, ProjRectifiedGridAngleGeoKey, TYPE_DOUBLE, 1,
2602 : ConvertAngularParam(poSRSCompatibleOfWKT1->GetNormProjParm(
2603 : SRS_PP_RECTIFIED_GRID_ANGLE, 0.0)));
2604 :
2605 1 : GTIFKeySet(
2606 : psGTIF, ProjScaleAtCenterGeoKey, TYPE_DOUBLE, 1,
2607 : poSRSCompatibleOfWKT1->GetNormProjParm(SRS_PP_SCALE_FACTOR, 1.0));
2608 :
2609 1 : GTIFKeySet(
2610 : psGTIF, ProjFalseEastingGeoKey, TYPE_DOUBLE, 1,
2611 : poSRSCompatibleOfWKT1->GetProjParm(SRS_PP_FALSE_EASTING, 0.0));
2612 :
2613 1 : GTIFKeySet(
2614 : psGTIF, ProjFalseNorthingGeoKey, TYPE_DOUBLE, 1,
2615 : poSRSCompatibleOfWKT1->GetProjParm(SRS_PP_FALSE_NORTHING, 0.0));
2616 : }
2617 :
2618 38 : else if (EQUAL(pszProjection,
2619 : SRS_PT_HOTINE_OBLIQUE_MERCATOR_AZIMUTH_CENTER))
2620 : {
2621 2 : GTIFKeySet(psGTIF, GTModelTypeGeoKey, TYPE_SHORT, 1,
2622 : ModelTypeProjected);
2623 2 : GTIFKeySet(psGTIF, ProjectedCSTypeGeoKey, TYPE_SHORT, 1, KvUserDefined);
2624 2 : GTIFKeySet(psGTIF, ProjectionGeoKey, TYPE_SHORT, 1, KvUserDefined);
2625 :
2626 2 : GTIFKeySet(psGTIF, ProjCoordTransGeoKey, TYPE_SHORT, 1,
2627 : CT_HotineObliqueMercatorAzimuthCenter);
2628 :
2629 2 : GTIFKeySet(psGTIF, ProjCenterLatGeoKey, TYPE_DOUBLE, 1,
2630 : ConvertAngularParam(poSRSCompatibleOfWKT1->GetNormProjParm(
2631 : SRS_PP_LATITUDE_OF_CENTER, 0.0)));
2632 :
2633 2 : GTIFKeySet(psGTIF, ProjCenterLongGeoKey, TYPE_DOUBLE, 1,
2634 : ConvertAngularParam(poSRSCompatibleOfWKT1->GetNormProjParm(
2635 : SRS_PP_LONGITUDE_OF_CENTER, 0.0)));
2636 :
2637 2 : GTIFKeySet(psGTIF, ProjAzimuthAngleGeoKey, TYPE_DOUBLE, 1,
2638 : poSRSCompatibleOfWKT1->GetNormProjParm(SRS_PP_AZIMUTH, 0.0));
2639 :
2640 2 : GTIFKeySet(psGTIF, ProjRectifiedGridAngleGeoKey, TYPE_DOUBLE, 1,
2641 : ConvertAngularParam(poSRSCompatibleOfWKT1->GetNormProjParm(
2642 : SRS_PP_RECTIFIED_GRID_ANGLE, 0.0)));
2643 :
2644 2 : GTIFKeySet(
2645 : psGTIF, ProjScaleAtCenterGeoKey, TYPE_DOUBLE, 1,
2646 : poSRSCompatibleOfWKT1->GetNormProjParm(SRS_PP_SCALE_FACTOR, 1.0));
2647 :
2648 2 : GTIFKeySet(
2649 : psGTIF, ProjFalseEastingGeoKey, TYPE_DOUBLE, 1,
2650 : poSRSCompatibleOfWKT1->GetProjParm(SRS_PP_FALSE_EASTING, 0.0));
2651 :
2652 2 : GTIFKeySet(
2653 : psGTIF, ProjFalseNorthingGeoKey, TYPE_DOUBLE, 1,
2654 : poSRSCompatibleOfWKT1->GetProjParm(SRS_PP_FALSE_NORTHING, 0.0));
2655 : }
2656 :
2657 36 : else if (EQUAL(pszProjection, "Laborde_Oblique_Mercator"))
2658 : {
2659 0 : GTIFKeySet(psGTIF, GTModelTypeGeoKey, TYPE_SHORT, 1,
2660 : ModelTypeProjected);
2661 0 : GTIFKeySet(psGTIF, ProjectedCSTypeGeoKey, TYPE_SHORT, 1, KvUserDefined);
2662 0 : GTIFKeySet(psGTIF, ProjectionGeoKey, TYPE_SHORT, 1, KvUserDefined);
2663 :
2664 0 : GTIFKeySet(psGTIF, ProjCoordTransGeoKey, TYPE_SHORT, 1,
2665 : CT_ObliqueMercator_Laborde);
2666 :
2667 0 : GTIFKeySet(psGTIF, ProjCenterLatGeoKey, TYPE_DOUBLE, 1,
2668 : ConvertAngularParam(poSRSCompatibleOfWKT1->GetNormProjParm(
2669 : SRS_PP_LATITUDE_OF_CENTER, 0.0)));
2670 :
2671 0 : GTIFKeySet(psGTIF, ProjCenterLongGeoKey, TYPE_DOUBLE, 1,
2672 : ConvertAngularParam(poSRSCompatibleOfWKT1->GetNormProjParm(
2673 : SRS_PP_LONGITUDE_OF_CENTER, 0.0)));
2674 :
2675 0 : GTIFKeySet(psGTIF, ProjAzimuthAngleGeoKey, TYPE_DOUBLE, 1,
2676 : poSRSCompatibleOfWKT1->GetNormProjParm(SRS_PP_AZIMUTH, 0.0));
2677 :
2678 0 : GTIFKeySet(
2679 : psGTIF, ProjScaleAtCenterGeoKey, TYPE_DOUBLE, 1,
2680 : poSRSCompatibleOfWKT1->GetNormProjParm(SRS_PP_SCALE_FACTOR, 1.0));
2681 :
2682 0 : GTIFKeySet(
2683 : psGTIF, ProjFalseEastingGeoKey, TYPE_DOUBLE, 1,
2684 : poSRSCompatibleOfWKT1->GetProjParm(SRS_PP_FALSE_EASTING, 0.0));
2685 :
2686 0 : GTIFKeySet(
2687 : psGTIF, ProjFalseNorthingGeoKey, TYPE_DOUBLE, 1,
2688 : poSRSCompatibleOfWKT1->GetProjParm(SRS_PP_FALSE_NORTHING, 0.0));
2689 : }
2690 :
2691 36 : else if (EQUAL(pszProjection, SRS_PT_CASSINI_SOLDNER))
2692 : {
2693 1 : GTIFKeySet(psGTIF, GTModelTypeGeoKey, TYPE_SHORT, 1,
2694 : ModelTypeProjected);
2695 1 : GTIFKeySet(psGTIF, ProjectedCSTypeGeoKey, TYPE_SHORT, 1, KvUserDefined);
2696 1 : GTIFKeySet(psGTIF, ProjectionGeoKey, TYPE_SHORT, 1, KvUserDefined);
2697 :
2698 1 : GTIFKeySet(psGTIF, ProjCoordTransGeoKey, TYPE_SHORT, 1,
2699 : CT_CassiniSoldner);
2700 :
2701 1 : GTIFKeySet(psGTIF, ProjNatOriginLatGeoKey, TYPE_DOUBLE, 1,
2702 : ConvertAngularParam(poSRSCompatibleOfWKT1->GetNormProjParm(
2703 : SRS_PP_LATITUDE_OF_ORIGIN, 0.0)));
2704 :
2705 1 : GTIFKeySet(psGTIF, ProjNatOriginLongGeoKey, TYPE_DOUBLE, 1,
2706 : ConvertAngularParam(poSRSCompatibleOfWKT1->GetNormProjParm(
2707 : SRS_PP_CENTRAL_MERIDIAN, 0.0)));
2708 :
2709 1 : GTIFKeySet(
2710 : psGTIF, ProjFalseEastingGeoKey, TYPE_DOUBLE, 1,
2711 : poSRSCompatibleOfWKT1->GetProjParm(SRS_PP_FALSE_EASTING, 0.0));
2712 :
2713 1 : GTIFKeySet(
2714 : psGTIF, ProjFalseNorthingGeoKey, TYPE_DOUBLE, 1,
2715 : poSRSCompatibleOfWKT1->GetProjParm(SRS_PP_FALSE_NORTHING, 0.0));
2716 : }
2717 :
2718 35 : else if (EQUAL(pszProjection, SRS_PT_EQUIDISTANT_CONIC))
2719 : {
2720 1 : GTIFKeySet(psGTIF, GTModelTypeGeoKey, TYPE_SHORT, 1,
2721 : ModelTypeProjected);
2722 1 : GTIFKeySet(psGTIF, ProjectedCSTypeGeoKey, TYPE_SHORT, 1, KvUserDefined);
2723 1 : GTIFKeySet(psGTIF, ProjectionGeoKey, TYPE_SHORT, 1, KvUserDefined);
2724 :
2725 1 : GTIFKeySet(psGTIF, ProjCoordTransGeoKey, TYPE_SHORT, 1,
2726 : CT_EquidistantConic);
2727 :
2728 1 : GTIFKeySet(psGTIF, ProjStdParallel1GeoKey, TYPE_DOUBLE, 1,
2729 : ConvertAngularParam(poSRSCompatibleOfWKT1->GetNormProjParm(
2730 : SRS_PP_STANDARD_PARALLEL_1, 0.0)));
2731 :
2732 1 : GTIFKeySet(psGTIF, ProjStdParallel2GeoKey, TYPE_DOUBLE, 1,
2733 : ConvertAngularParam(poSRSCompatibleOfWKT1->GetNormProjParm(
2734 : SRS_PP_STANDARD_PARALLEL_2, 0.0)));
2735 :
2736 1 : GTIFKeySet(psGTIF, ProjNatOriginLatGeoKey, TYPE_DOUBLE, 1,
2737 : ConvertAngularParam(poSRSCompatibleOfWKT1->GetNormProjParm(
2738 : SRS_PP_LATITUDE_OF_CENTER, 0.0)));
2739 :
2740 1 : GTIFKeySet(psGTIF, ProjNatOriginLongGeoKey, TYPE_DOUBLE, 1,
2741 : ConvertAngularParam(poSRSCompatibleOfWKT1->GetNormProjParm(
2742 : SRS_PP_LONGITUDE_OF_CENTER, 0.0)));
2743 :
2744 1 : GTIFKeySet(
2745 : psGTIF, ProjFalseEastingGeoKey, TYPE_DOUBLE, 1,
2746 : poSRSCompatibleOfWKT1->GetProjParm(SRS_PP_FALSE_EASTING, 0.0));
2747 :
2748 1 : GTIFKeySet(
2749 : psGTIF, ProjFalseNorthingGeoKey, TYPE_DOUBLE, 1,
2750 : poSRSCompatibleOfWKT1->GetProjParm(SRS_PP_FALSE_NORTHING, 0.0));
2751 : }
2752 :
2753 34 : else if (EQUAL(pszProjection, SRS_PT_POLYCONIC))
2754 : {
2755 1 : GTIFKeySet(psGTIF, GTModelTypeGeoKey, TYPE_SHORT, 1,
2756 : ModelTypeProjected);
2757 1 : GTIFKeySet(psGTIF, ProjectedCSTypeGeoKey, TYPE_SHORT, 1, KvUserDefined);
2758 1 : GTIFKeySet(psGTIF, ProjectionGeoKey, TYPE_SHORT, 1, KvUserDefined);
2759 :
2760 1 : GTIFKeySet(psGTIF, ProjCoordTransGeoKey, TYPE_SHORT, 1, CT_Polyconic);
2761 :
2762 1 : GTIFKeySet(psGTIF, ProjNatOriginLatGeoKey, TYPE_DOUBLE, 1,
2763 : poSRSCompatibleOfWKT1->GetNormProjParm(
2764 : SRS_PP_LATITUDE_OF_ORIGIN, 0.0));
2765 :
2766 1 : GTIFKeySet(psGTIF, ProjNatOriginLongGeoKey, TYPE_DOUBLE, 1,
2767 : poSRSCompatibleOfWKT1->GetNormProjParm(
2768 : SRS_PP_CENTRAL_MERIDIAN, 0.0));
2769 :
2770 1 : GTIFKeySet(
2771 : psGTIF, ProjScaleAtNatOriginGeoKey, TYPE_DOUBLE, 1,
2772 : poSRSCompatibleOfWKT1->GetNormProjParm(SRS_PP_SCALE_FACTOR, 1.0));
2773 :
2774 1 : GTIFKeySet(
2775 : psGTIF, ProjFalseEastingGeoKey, TYPE_DOUBLE, 1,
2776 : poSRSCompatibleOfWKT1->GetProjParm(SRS_PP_FALSE_EASTING, 0.0));
2777 :
2778 1 : GTIFKeySet(
2779 : psGTIF, ProjFalseNorthingGeoKey, TYPE_DOUBLE, 1,
2780 : poSRSCompatibleOfWKT1->GetProjParm(SRS_PP_FALSE_NORTHING, 0.0));
2781 : }
2782 :
2783 33 : else if (EQUAL(pszProjection, SRS_PT_AZIMUTHAL_EQUIDISTANT))
2784 : {
2785 3 : GTIFKeySet(psGTIF, GTModelTypeGeoKey, TYPE_SHORT, 1,
2786 : ModelTypeProjected);
2787 3 : GTIFKeySet(psGTIF, ProjectedCSTypeGeoKey, TYPE_SHORT, 1, KvUserDefined);
2788 3 : GTIFKeySet(psGTIF, ProjectionGeoKey, TYPE_SHORT, 1, KvUserDefined);
2789 :
2790 3 : GTIFKeySet(psGTIF, ProjCoordTransGeoKey, TYPE_SHORT, 1,
2791 : CT_AzimuthalEquidistant);
2792 :
2793 3 : GTIFKeySet(psGTIF, ProjCenterLatGeoKey, TYPE_DOUBLE, 1,
2794 : ConvertAngularParam(poSRSCompatibleOfWKT1->GetNormProjParm(
2795 : SRS_PP_LATITUDE_OF_CENTER, 0.0)));
2796 :
2797 3 : GTIFKeySet(psGTIF, ProjCenterLongGeoKey, TYPE_DOUBLE, 1,
2798 : ConvertAngularParam(poSRSCompatibleOfWKT1->GetNormProjParm(
2799 : SRS_PP_LONGITUDE_OF_CENTER, 0.0)));
2800 :
2801 3 : GTIFKeySet(
2802 : psGTIF, ProjFalseEastingGeoKey, TYPE_DOUBLE, 1,
2803 : poSRSCompatibleOfWKT1->GetProjParm(SRS_PP_FALSE_EASTING, 0.0));
2804 :
2805 3 : GTIFKeySet(
2806 : psGTIF, ProjFalseNorthingGeoKey, TYPE_DOUBLE, 1,
2807 : poSRSCompatibleOfWKT1->GetProjParm(SRS_PP_FALSE_NORTHING, 0.0));
2808 : }
2809 :
2810 30 : else if (EQUAL(pszProjection, SRS_PT_MILLER_CYLINDRICAL))
2811 : {
2812 1 : GTIFKeySet(psGTIF, GTModelTypeGeoKey, TYPE_SHORT, 1,
2813 : ModelTypeProjected);
2814 1 : GTIFKeySet(psGTIF, ProjectedCSTypeGeoKey, TYPE_SHORT, 1, KvUserDefined);
2815 1 : GTIFKeySet(psGTIF, ProjectionGeoKey, TYPE_SHORT, 1, KvUserDefined);
2816 :
2817 1 : GTIFKeySet(psGTIF, ProjCoordTransGeoKey, TYPE_SHORT, 1,
2818 : CT_MillerCylindrical);
2819 :
2820 1 : GTIFKeySet(psGTIF, ProjCenterLatGeoKey, TYPE_DOUBLE, 1,
2821 : ConvertAngularParam(poSRSCompatibleOfWKT1->GetNormProjParm(
2822 : SRS_PP_LATITUDE_OF_CENTER, 0.0)));
2823 :
2824 1 : GTIFKeySet(psGTIF, ProjCenterLongGeoKey, TYPE_DOUBLE, 1,
2825 : ConvertAngularParam(poSRSCompatibleOfWKT1->GetNormProjParm(
2826 : SRS_PP_LONGITUDE_OF_CENTER, 0.0)));
2827 :
2828 1 : GTIFKeySet(
2829 : psGTIF, ProjFalseEastingGeoKey, TYPE_DOUBLE, 1,
2830 : poSRSCompatibleOfWKT1->GetProjParm(SRS_PP_FALSE_EASTING, 0.0));
2831 :
2832 1 : GTIFKeySet(
2833 : psGTIF, ProjFalseNorthingGeoKey, TYPE_DOUBLE, 1,
2834 : poSRSCompatibleOfWKT1->GetProjParm(SRS_PP_FALSE_NORTHING, 0.0));
2835 : }
2836 :
2837 29 : else if (EQUAL(pszProjection, SRS_PT_EQUIRECTANGULAR))
2838 : {
2839 9 : GTIFKeySet(psGTIF, GTModelTypeGeoKey, TYPE_SHORT, 1,
2840 : ModelTypeProjected);
2841 9 : GTIFKeySet(psGTIF, ProjectedCSTypeGeoKey, TYPE_SHORT, 1, KvUserDefined);
2842 9 : GTIFKeySet(psGTIF, ProjectionGeoKey, TYPE_SHORT, 1, KvUserDefined);
2843 :
2844 9 : GTIFKeySet(psGTIF, ProjCoordTransGeoKey, TYPE_SHORT, 1,
2845 : CT_Equirectangular);
2846 :
2847 9 : GTIFKeySet(psGTIF, ProjCenterLatGeoKey, TYPE_DOUBLE, 1,
2848 : ConvertAngularParam(poSRSCompatibleOfWKT1->GetNormProjParm(
2849 : SRS_PP_LATITUDE_OF_ORIGIN, 0.0)));
2850 :
2851 9 : GTIFKeySet(psGTIF, ProjCenterLongGeoKey, TYPE_DOUBLE, 1,
2852 : ConvertAngularParam(poSRSCompatibleOfWKT1->GetNormProjParm(
2853 : SRS_PP_CENTRAL_MERIDIAN, 0.0)));
2854 :
2855 9 : GTIFKeySet(psGTIF, ProjStdParallel1GeoKey, TYPE_DOUBLE, 1,
2856 : ConvertAngularParam(poSRSCompatibleOfWKT1->GetNormProjParm(
2857 : SRS_PP_STANDARD_PARALLEL_1, 0.0)));
2858 :
2859 9 : GTIFKeySet(
2860 : psGTIF, ProjFalseEastingGeoKey, TYPE_DOUBLE, 1,
2861 : poSRSCompatibleOfWKT1->GetProjParm(SRS_PP_FALSE_EASTING, 0.0));
2862 :
2863 9 : GTIFKeySet(
2864 : psGTIF, ProjFalseNorthingGeoKey, TYPE_DOUBLE, 1,
2865 : poSRSCompatibleOfWKT1->GetProjParm(SRS_PP_FALSE_NORTHING, 0.0));
2866 : }
2867 :
2868 20 : else if (EQUAL(pszProjection, SRS_PT_GNOMONIC))
2869 : {
2870 1 : GTIFKeySet(psGTIF, GTModelTypeGeoKey, TYPE_SHORT, 1,
2871 : ModelTypeProjected);
2872 1 : GTIFKeySet(psGTIF, ProjectedCSTypeGeoKey, TYPE_SHORT, 1, KvUserDefined);
2873 1 : GTIFKeySet(psGTIF, ProjectionGeoKey, TYPE_SHORT, 1, KvUserDefined);
2874 :
2875 1 : GTIFKeySet(psGTIF, ProjCoordTransGeoKey, TYPE_SHORT, 1, CT_Gnomonic);
2876 :
2877 1 : GTIFKeySet(psGTIF, ProjCenterLatGeoKey, TYPE_DOUBLE, 1,
2878 : ConvertAngularParam(poSRSCompatibleOfWKT1->GetNormProjParm(
2879 : SRS_PP_LATITUDE_OF_ORIGIN, 0.0)));
2880 :
2881 1 : GTIFKeySet(psGTIF, ProjCenterLongGeoKey, TYPE_DOUBLE, 1,
2882 : ConvertAngularParam(poSRSCompatibleOfWKT1->GetNormProjParm(
2883 : SRS_PP_CENTRAL_MERIDIAN, 0.0)));
2884 :
2885 1 : GTIFKeySet(
2886 : psGTIF, ProjFalseEastingGeoKey, TYPE_DOUBLE, 1,
2887 : poSRSCompatibleOfWKT1->GetProjParm(SRS_PP_FALSE_EASTING, 0.0));
2888 :
2889 1 : GTIFKeySet(
2890 : psGTIF, ProjFalseNorthingGeoKey, TYPE_DOUBLE, 1,
2891 : poSRSCompatibleOfWKT1->GetProjParm(SRS_PP_FALSE_NORTHING, 0.0));
2892 : }
2893 :
2894 19 : else if (EQUAL(pszProjection, SRS_PT_LAMBERT_AZIMUTHAL_EQUAL_AREA))
2895 : {
2896 2 : GTIFKeySet(psGTIF, GTModelTypeGeoKey, TYPE_SHORT, 1,
2897 : ModelTypeProjected);
2898 2 : GTIFKeySet(psGTIF, ProjectedCSTypeGeoKey, TYPE_SHORT, 1, KvUserDefined);
2899 2 : GTIFKeySet(psGTIF, ProjectionGeoKey, TYPE_SHORT, 1, KvUserDefined);
2900 :
2901 2 : GTIFKeySet(psGTIF, ProjCoordTransGeoKey, TYPE_SHORT, 1,
2902 : CT_LambertAzimEqualArea);
2903 :
2904 2 : GTIFKeySet(psGTIF, ProjCenterLatGeoKey, TYPE_DOUBLE, 1,
2905 : ConvertAngularParam(poSRSCompatibleOfWKT1->GetNormProjParm(
2906 : SRS_PP_LATITUDE_OF_CENTER, 0.0)));
2907 :
2908 2 : GTIFKeySet(psGTIF, ProjCenterLongGeoKey, TYPE_DOUBLE, 1,
2909 : ConvertAngularParam(poSRSCompatibleOfWKT1->GetNormProjParm(
2910 : SRS_PP_LONGITUDE_OF_CENTER, 0.0)));
2911 :
2912 2 : GTIFKeySet(
2913 : psGTIF, ProjFalseEastingGeoKey, TYPE_DOUBLE, 1,
2914 : poSRSCompatibleOfWKT1->GetProjParm(SRS_PP_FALSE_EASTING, 0.0));
2915 :
2916 2 : GTIFKeySet(
2917 : psGTIF, ProjFalseNorthingGeoKey, TYPE_DOUBLE, 1,
2918 : poSRSCompatibleOfWKT1->GetProjParm(SRS_PP_FALSE_NORTHING, 0.0));
2919 : }
2920 :
2921 17 : else if (EQUAL(pszProjection, SRS_PT_ORTHOGRAPHIC))
2922 : {
2923 0 : GTIFKeySet(psGTIF, GTModelTypeGeoKey, TYPE_SHORT, 1,
2924 : ModelTypeProjected);
2925 0 : GTIFKeySet(psGTIF, ProjectedCSTypeGeoKey, TYPE_SHORT, 1, KvUserDefined);
2926 0 : GTIFKeySet(psGTIF, ProjectionGeoKey, TYPE_SHORT, 1, KvUserDefined);
2927 :
2928 0 : GTIFKeySet(psGTIF, ProjCoordTransGeoKey, TYPE_SHORT, 1,
2929 : CT_Orthographic);
2930 :
2931 0 : GTIFKeySet(psGTIF, ProjCenterLatGeoKey, TYPE_DOUBLE, 1,
2932 : ConvertAngularParam(poSRSCompatibleOfWKT1->GetNormProjParm(
2933 : SRS_PP_LATITUDE_OF_ORIGIN, 0.0)));
2934 :
2935 0 : GTIFKeySet(psGTIF, ProjCenterLongGeoKey, TYPE_DOUBLE, 1,
2936 : ConvertAngularParam(poSRSCompatibleOfWKT1->GetNormProjParm(
2937 : SRS_PP_CENTRAL_MERIDIAN, 0.0)));
2938 :
2939 0 : GTIFKeySet(
2940 : psGTIF, ProjFalseEastingGeoKey, TYPE_DOUBLE, 1,
2941 : poSRSCompatibleOfWKT1->GetProjParm(SRS_PP_FALSE_EASTING, 0.0));
2942 :
2943 0 : GTIFKeySet(
2944 : psGTIF, ProjFalseNorthingGeoKey, TYPE_DOUBLE, 1,
2945 : poSRSCompatibleOfWKT1->GetProjParm(SRS_PP_FALSE_NORTHING, 0.0));
2946 : }
2947 :
2948 17 : else if (EQUAL(pszProjection, SRS_PT_NEW_ZEALAND_MAP_GRID))
2949 : {
2950 1 : GTIFKeySet(psGTIF, GTModelTypeGeoKey, TYPE_SHORT, 1,
2951 : ModelTypeProjected);
2952 1 : GTIFKeySet(psGTIF, ProjectedCSTypeGeoKey, TYPE_SHORT, 1, KvUserDefined);
2953 1 : GTIFKeySet(psGTIF, ProjectionGeoKey, TYPE_SHORT, 1, KvUserDefined);
2954 :
2955 1 : GTIFKeySet(psGTIF, ProjCoordTransGeoKey, TYPE_SHORT, 1,
2956 : CT_NewZealandMapGrid);
2957 :
2958 1 : GTIFKeySet(psGTIF, ProjNatOriginLatGeoKey, TYPE_DOUBLE, 1,
2959 : ConvertAngularParam(poSRSCompatibleOfWKT1->GetNormProjParm(
2960 : SRS_PP_LATITUDE_OF_ORIGIN, 0.0)));
2961 :
2962 1 : GTIFKeySet(psGTIF, ProjNatOriginLongGeoKey, TYPE_DOUBLE, 1,
2963 : ConvertAngularParam(poSRSCompatibleOfWKT1->GetNormProjParm(
2964 : SRS_PP_CENTRAL_MERIDIAN, 0.0)));
2965 :
2966 1 : GTIFKeySet(
2967 : psGTIF, ProjFalseEastingGeoKey, TYPE_DOUBLE, 1,
2968 : poSRSCompatibleOfWKT1->GetProjParm(SRS_PP_FALSE_EASTING, 0.0));
2969 :
2970 1 : GTIFKeySet(
2971 : psGTIF, ProjFalseNorthingGeoKey, TYPE_DOUBLE, 1,
2972 : poSRSCompatibleOfWKT1->GetProjParm(SRS_PP_FALSE_NORTHING, 0.0));
2973 : }
2974 :
2975 16 : else if (EQUAL(pszProjection, SRS_PT_ROBINSON))
2976 : {
2977 1 : GTIFKeySet(psGTIF, GTModelTypeGeoKey, TYPE_SHORT, 1,
2978 : ModelTypeProjected);
2979 1 : GTIFKeySet(psGTIF, ProjectedCSTypeGeoKey, TYPE_SHORT, 1, KvUserDefined);
2980 1 : GTIFKeySet(psGTIF, ProjectionGeoKey, TYPE_SHORT, 1, KvUserDefined);
2981 :
2982 1 : GTIFKeySet(psGTIF, ProjCoordTransGeoKey, TYPE_SHORT, 1, CT_Robinson);
2983 :
2984 1 : GTIFKeySet(psGTIF, ProjCenterLongGeoKey, TYPE_DOUBLE, 1,
2985 : ConvertAngularParam(poSRSCompatibleOfWKT1->GetNormProjParm(
2986 : SRS_PP_LONGITUDE_OF_CENTER, 0.0)));
2987 :
2988 1 : GTIFKeySet(
2989 : psGTIF, ProjFalseEastingGeoKey, TYPE_DOUBLE, 1,
2990 : poSRSCompatibleOfWKT1->GetProjParm(SRS_PP_FALSE_EASTING, 0.0));
2991 :
2992 1 : GTIFKeySet(
2993 : psGTIF, ProjFalseNorthingGeoKey, TYPE_DOUBLE, 1,
2994 : poSRSCompatibleOfWKT1->GetProjParm(SRS_PP_FALSE_NORTHING, 0.0));
2995 : }
2996 :
2997 15 : else if (EQUAL(pszProjection, SRS_PT_SINUSOIDAL))
2998 : {
2999 4 : GTIFKeySet(psGTIF, GTModelTypeGeoKey, TYPE_SHORT, 1,
3000 : ModelTypeProjected);
3001 4 : GTIFKeySet(psGTIF, ProjectedCSTypeGeoKey, TYPE_SHORT, 1, KvUserDefined);
3002 4 : GTIFKeySet(psGTIF, ProjectionGeoKey, TYPE_SHORT, 1, KvUserDefined);
3003 :
3004 4 : GTIFKeySet(psGTIF, ProjCoordTransGeoKey, TYPE_SHORT, 1, CT_Sinusoidal);
3005 :
3006 4 : GTIFKeySet(psGTIF, ProjCenterLongGeoKey, TYPE_DOUBLE, 1,
3007 : ConvertAngularParam(poSRSCompatibleOfWKT1->GetNormProjParm(
3008 : SRS_PP_LONGITUDE_OF_CENTER, 0.0)));
3009 :
3010 4 : GTIFKeySet(
3011 : psGTIF, ProjFalseEastingGeoKey, TYPE_DOUBLE, 1,
3012 : poSRSCompatibleOfWKT1->GetProjParm(SRS_PP_FALSE_EASTING, 0.0));
3013 :
3014 4 : GTIFKeySet(
3015 : psGTIF, ProjFalseNorthingGeoKey, TYPE_DOUBLE, 1,
3016 : poSRSCompatibleOfWKT1->GetProjParm(SRS_PP_FALSE_NORTHING, 0.0));
3017 : }
3018 :
3019 11 : else if (EQUAL(pszProjection, SRS_PT_VANDERGRINTEN))
3020 : {
3021 2 : GTIFKeySet(psGTIF, GTModelTypeGeoKey, TYPE_SHORT, 1,
3022 : ModelTypeProjected);
3023 2 : GTIFKeySet(psGTIF, ProjectedCSTypeGeoKey, TYPE_SHORT, 1, KvUserDefined);
3024 2 : GTIFKeySet(psGTIF, ProjectionGeoKey, TYPE_SHORT, 1, KvUserDefined);
3025 :
3026 2 : GTIFKeySet(psGTIF, ProjCoordTransGeoKey, TYPE_SHORT, 1,
3027 : CT_VanDerGrinten);
3028 :
3029 2 : GTIFKeySet(psGTIF, ProjCenterLongGeoKey, TYPE_DOUBLE, 1,
3030 : ConvertAngularParam(poSRSCompatibleOfWKT1->GetNormProjParm(
3031 : SRS_PP_CENTRAL_MERIDIAN, 0.0)));
3032 :
3033 2 : GTIFKeySet(
3034 : psGTIF, ProjFalseEastingGeoKey, TYPE_DOUBLE, 1,
3035 : poSRSCompatibleOfWKT1->GetProjParm(SRS_PP_FALSE_EASTING, 0.0));
3036 :
3037 2 : GTIFKeySet(
3038 : psGTIF, ProjFalseNorthingGeoKey, TYPE_DOUBLE, 1,
3039 : poSRSCompatibleOfWKT1->GetProjParm(SRS_PP_FALSE_NORTHING, 0.0));
3040 : }
3041 :
3042 9 : else if (EQUAL(pszProjection, SRS_PT_LAMBERT_CONFORMAL_CONIC_2SP))
3043 : {
3044 2 : GTIFKeySet(psGTIF, GTModelTypeGeoKey, TYPE_SHORT, 1,
3045 : ModelTypeProjected);
3046 2 : GTIFKeySet(psGTIF, ProjectedCSTypeGeoKey, TYPE_SHORT, 1, KvUserDefined);
3047 2 : GTIFKeySet(psGTIF, ProjectionGeoKey, TYPE_SHORT, 1, KvUserDefined);
3048 :
3049 2 : GTIFKeySet(psGTIF, ProjCoordTransGeoKey, TYPE_SHORT, 1,
3050 : CT_LambertConfConic_2SP);
3051 :
3052 2 : GTIFKeySet(psGTIF, ProjFalseOriginLatGeoKey, TYPE_DOUBLE, 1,
3053 : ConvertAngularParam(poSRSCompatibleOfWKT1->GetNormProjParm(
3054 : SRS_PP_LATITUDE_OF_ORIGIN, 0.0)));
3055 :
3056 2 : GTIFKeySet(psGTIF, ProjFalseOriginLongGeoKey, TYPE_DOUBLE, 1,
3057 : ConvertAngularParam(poSRSCompatibleOfWKT1->GetNormProjParm(
3058 : SRS_PP_CENTRAL_MERIDIAN, 0.0)));
3059 :
3060 2 : GTIFKeySet(psGTIF, ProjStdParallel1GeoKey, TYPE_DOUBLE, 1,
3061 : ConvertAngularParam(poSRSCompatibleOfWKT1->GetNormProjParm(
3062 : SRS_PP_STANDARD_PARALLEL_1, 0.0)));
3063 :
3064 2 : GTIFKeySet(psGTIF, ProjStdParallel2GeoKey, TYPE_DOUBLE, 1,
3065 : ConvertAngularParam(poSRSCompatibleOfWKT1->GetNormProjParm(
3066 : SRS_PP_STANDARD_PARALLEL_2, 0.0)));
3067 :
3068 2 : GTIFKeySet(
3069 : psGTIF, ProjFalseOriginEastingGeoKey, TYPE_DOUBLE, 1,
3070 : poSRSCompatibleOfWKT1->GetProjParm(SRS_PP_FALSE_EASTING, 0.0));
3071 :
3072 2 : GTIFKeySet(
3073 : psGTIF, ProjFalseOriginNorthingGeoKey, TYPE_DOUBLE, 1,
3074 : poSRSCompatibleOfWKT1->GetProjParm(SRS_PP_FALSE_NORTHING, 0.0));
3075 : }
3076 :
3077 7 : else if (EQUAL(pszProjection, SRS_PT_LAMBERT_CONFORMAL_CONIC_1SP))
3078 : {
3079 3 : GTIFKeySet(psGTIF, GTModelTypeGeoKey, TYPE_SHORT, 1,
3080 : ModelTypeProjected);
3081 3 : GTIFKeySet(psGTIF, ProjectedCSTypeGeoKey, TYPE_SHORT, 1, KvUserDefined);
3082 3 : GTIFKeySet(psGTIF, ProjectionGeoKey, TYPE_SHORT, 1, KvUserDefined);
3083 :
3084 3 : GTIFKeySet(psGTIF, ProjCoordTransGeoKey, TYPE_SHORT, 1,
3085 : CT_LambertConfConic_1SP);
3086 :
3087 3 : GTIFKeySet(psGTIF, ProjNatOriginLatGeoKey, TYPE_DOUBLE, 1,
3088 : ConvertAngularParam(poSRSCompatibleOfWKT1->GetNormProjParm(
3089 : SRS_PP_LATITUDE_OF_ORIGIN, 0.0)));
3090 :
3091 3 : GTIFKeySet(psGTIF, ProjNatOriginLongGeoKey, TYPE_DOUBLE, 1,
3092 : ConvertAngularParam(poSRSCompatibleOfWKT1->GetNormProjParm(
3093 : SRS_PP_CENTRAL_MERIDIAN, 0.0)));
3094 :
3095 3 : GTIFKeySet(
3096 : psGTIF, ProjScaleAtNatOriginGeoKey, TYPE_DOUBLE, 1,
3097 : poSRSCompatibleOfWKT1->GetNormProjParm(SRS_PP_SCALE_FACTOR, 1.0));
3098 :
3099 3 : GTIFKeySet(
3100 : psGTIF, ProjFalseEastingGeoKey, TYPE_DOUBLE, 1,
3101 : poSRSCompatibleOfWKT1->GetProjParm(SRS_PP_FALSE_EASTING, 0.0));
3102 :
3103 3 : GTIFKeySet(
3104 : psGTIF, ProjFalseNorthingGeoKey, TYPE_DOUBLE, 1,
3105 : poSRSCompatibleOfWKT1->GetProjParm(SRS_PP_FALSE_NORTHING, 0.0));
3106 : }
3107 :
3108 4 : else if (EQUAL(pszProjection, SRS_PT_CYLINDRICAL_EQUAL_AREA))
3109 : {
3110 1 : GTIFKeySet(psGTIF, GTModelTypeGeoKey, TYPE_SHORT, 1,
3111 : ModelTypeProjected);
3112 1 : GTIFKeySet(psGTIF, ProjectedCSTypeGeoKey, TYPE_SHORT, 1, KvUserDefined);
3113 1 : GTIFKeySet(psGTIF, ProjectionGeoKey, TYPE_SHORT, 1, KvUserDefined);
3114 :
3115 1 : GTIFKeySet(psGTIF, ProjCoordTransGeoKey, TYPE_SHORT, 1,
3116 : CT_CylindricalEqualArea);
3117 :
3118 1 : GTIFKeySet(psGTIF, ProjNatOriginLongGeoKey, TYPE_DOUBLE, 1,
3119 : ConvertAngularParam(poSRSCompatibleOfWKT1->GetNormProjParm(
3120 : SRS_PP_CENTRAL_MERIDIAN, 0.0)));
3121 :
3122 1 : GTIFKeySet(psGTIF, ProjStdParallel1GeoKey, TYPE_DOUBLE, 1,
3123 : ConvertAngularParam(poSRSCompatibleOfWKT1->GetNormProjParm(
3124 : SRS_PP_STANDARD_PARALLEL_1, 0.0)));
3125 :
3126 1 : GTIFKeySet(
3127 : psGTIF, ProjFalseEastingGeoKey, TYPE_DOUBLE, 1,
3128 : poSRSCompatibleOfWKT1->GetProjParm(SRS_PP_FALSE_EASTING, 0.0));
3129 :
3130 1 : GTIFKeySet(
3131 : psGTIF, ProjFalseNorthingGeoKey, TYPE_DOUBLE, 1,
3132 : poSRSCompatibleOfWKT1->GetProjParm(SRS_PP_FALSE_NORTHING, 0.0));
3133 : }
3134 : else
3135 : {
3136 3 : bWritePEString = true;
3137 3 : bUnknownProjection = true;
3138 : }
3139 :
3140 : // Note that VERTCS is an ESRI "spelling" of VERT_CS so we assume if
3141 : // we find it that we should try to treat this as a PE string.
3142 6775 : if (eFlavor == GEOTIFF_KEYS_ESRI_PE ||
3143 3387 : poSRS->GetAttrValue("VERTCS") != nullptr)
3144 : {
3145 1 : bWritePEString = true;
3146 : }
3147 :
3148 3388 : if (nPCS == KvUserDefined)
3149 : {
3150 : const char *pszPROJ4Ext =
3151 2363 : poSRS->GetExtension("PROJCS", "PROJ4", nullptr);
3152 2363 : if (pszPROJ4Ext &&
3153 4 : strstr(pszPROJ4Ext, "+proj=merc +a=6378137 +b=6378137"))
3154 : {
3155 3 : bWritePEString = true;
3156 : }
3157 : }
3158 :
3159 3388 : bWritePEString &=
3160 3388 : CPLTestBool(CPLGetConfigOption("GTIFF_ESRI_CITATION", "YES"));
3161 :
3162 3388 : bool peStrStored = false;
3163 :
3164 3388 : if (bWritePEString)
3165 : {
3166 : // Anything we can't map, store as an ESRI PE string with a citation
3167 : // key.
3168 7 : char *pszPEString = nullptr;
3169 : // We cheat a bit, but if we have a custom_proj4, do not morph to ESRI
3170 : // so as to keep the EXTENSION PROJ4 node
3171 7 : const char *const apszOptionsDefault[] = {nullptr};
3172 7 : const char *const apszOptionsEsri[] = {"FORMAT=WKT1_ESRI", nullptr};
3173 7 : const char *const *papszOptions = apszOptionsDefault;
3174 17 : if (!(bUnknownProjection &&
3175 14 : poSRS->GetExtension("PROJCS", "PROJ4", nullptr) != nullptr) &&
3176 7 : !(poSRS->IsProjected() && !poSRS->IsCompound() &&
3177 7 : poSRS->GetAxesCount() == 3))
3178 : {
3179 7 : papszOptions = apszOptionsEsri;
3180 : }
3181 7 : poSRS->exportToWkt(&pszPEString, papszOptions);
3182 7 : const int peStrLen = static_cast<int>(strlen(pszPEString));
3183 7 : if (peStrLen > 0)
3184 : {
3185 : char *outPeStr = static_cast<char *>(
3186 7 : CPLMalloc(peStrLen + strlen("ESRI PE String = ") + 1));
3187 7 : strcpy(outPeStr, "ESRI PE String = ");
3188 7 : strcat(outPeStr, pszPEString);
3189 7 : oMapAsciiKeys[PCSCitationGeoKey] = outPeStr;
3190 7 : peStrStored = true;
3191 7 : CPLFree(outPeStr);
3192 : }
3193 7 : CPLFree(pszPEString);
3194 7 : GTIFKeySet(psGTIF, GTModelTypeGeoKey, TYPE_SHORT, 1, KvUserDefined);
3195 :
3196 : // Not completely sure we really need to imitate ArcGIS to that point
3197 : // but that cannot hurt.
3198 7 : if (nPCS == 3857)
3199 : {
3200 1 : oMapAsciiKeys[GTCitationGeoKey] =
3201 1 : "PCS Name = WGS_1984_Web_Mercator_Auxiliary_Sphere";
3202 1 : GTIFKeySet(psGTIF, GeographicTypeGeoKey, TYPE_SHORT, 1, GCS_WGS_84);
3203 1 : GTIFKeySet(psGTIF, GeogSemiMajorAxisGeoKey, TYPE_DOUBLE, 1,
3204 : 6378137.0);
3205 1 : GTIFKeySet(psGTIF, GeogInvFlatteningGeoKey, TYPE_DOUBLE, 1,
3206 : 298.257223563);
3207 : }
3208 : }
3209 :
3210 : /* -------------------------------------------------------------------- */
3211 : /* Is there a false easting/northing set? If so, write out a */
3212 : /* special geokey tag to indicate that GDAL has written these */
3213 : /* with the proper interpretation of the linear units. */
3214 : /* -------------------------------------------------------------------- */
3215 3388 : double dfFE = 0.0;
3216 3388 : double dfFN = 0.0;
3217 :
3218 6758 : if (eVersion == GEOTIFF_VERSION_1_0 &&
3219 3370 : (GDALGTIFKeyGetDOUBLE(psGTIF, ProjFalseEastingGeoKey, &dfFE, 0, 1) ||
3220 3314 : GDALGTIFKeyGetDOUBLE(psGTIF, ProjFalseNorthingGeoKey, &dfFN, 0, 1) ||
3221 3314 : GDALGTIFKeyGetDOUBLE(psGTIF, ProjFalseOriginEastingGeoKey, &dfFE, 0,
3222 3312 : 1) ||
3223 3312 : GDALGTIFKeyGetDOUBLE(psGTIF, ProjFalseOriginNorthingGeoKey, &dfFN, 0,
3224 58 : 1)) &&
3225 6758 : (dfFE != 0.0 || dfFN != 0.0) && nUOMLengthCode != 9001)
3226 : {
3227 2 : GTIFKeySet(psGTIF, ProjLinearUnitsInterpCorrectGeoKey, TYPE_SHORT, 1,
3228 : static_cast<short>(1));
3229 : }
3230 :
3231 : /* -------------------------------------------------------------------- */
3232 : /* Write linear units information. */
3233 : /* -------------------------------------------------------------------- */
3234 3388 : if (poSRS->IsGeocentric())
3235 : {
3236 2 : GTIFKeySet(psGTIF, GeogLinearUnitsGeoKey, TYPE_SHORT, 1,
3237 : nUOMLengthCode);
3238 2 : if (nUOMLengthCode == KvUserDefined)
3239 0 : GTIFKeySet(psGTIF, GeogLinearUnitSizeGeoKey, TYPE_DOUBLE, 1,
3240 : dfLinearUOM);
3241 : }
3242 4411 : else if (!poSRS->IsGeographic() &&
3243 1025 : (nPCS == KvUserDefined || eVersion == GEOTIFF_VERSION_1_0))
3244 : {
3245 1198 : GTIFKeySet(psGTIF, ProjLinearUnitsGeoKey, TYPE_SHORT, 1,
3246 : nUOMLengthCode);
3247 1198 : if (nUOMLengthCode == KvUserDefined)
3248 2 : GTIFKeySet(psGTIF, ProjLinearUnitSizeGeoKey, TYPE_DOUBLE, 1,
3249 : dfLinearUOM);
3250 :
3251 : // If linear units name is available and user defined, store it as
3252 : // citation.
3253 1191 : if (!peStrStored && nUOMLengthCode == KvUserDefined &&
3254 2391 : !osLinearUOMName.empty() &&
3255 2 : CPLTestBool(CPLGetConfigOption("GTIFF_ESRI_CITATION", "YES")))
3256 : {
3257 2 : SetLinearUnitCitation(oMapAsciiKeys, osLinearUOMName.c_str());
3258 : }
3259 : }
3260 :
3261 : /* -------------------------------------------------------------------- */
3262 : /* Write angular units. */
3263 : /* -------------------------------------------------------------------- */
3264 :
3265 3388 : if (bHasEllipsoid &&
3266 2669 : (nGCS == KvUserDefined || eVersion == GEOTIFF_VERSION_1_0))
3267 : {
3268 3368 : if (EQUAL(osAngUnitName.c_str(), "Degree"))
3269 3357 : GTIFKeySet(psGTIF, GeogAngularUnitsGeoKey, TYPE_SHORT, 1,
3270 : Angular_Degree);
3271 11 : else if (EQUAL(osAngUnitName.c_str(), "arc-second"))
3272 1 : GTIFKeySet(psGTIF, GeogAngularUnitsGeoKey, TYPE_SHORT, 1,
3273 : Angular_Arc_Second);
3274 10 : else if (EQUAL(osAngUnitName.c_str(), "arc-minute"))
3275 1 : GTIFKeySet(psGTIF, GeogAngularUnitsGeoKey, TYPE_SHORT, 1,
3276 : Angular_Arc_Minute);
3277 9 : else if (EQUAL(osAngUnitName.c_str(), "grad"))
3278 4 : GTIFKeySet(psGTIF, GeogAngularUnitsGeoKey, TYPE_SHORT, 1,
3279 : Angular_Grad);
3280 5 : else if (EQUAL(osAngUnitName.c_str(), "gon"))
3281 1 : GTIFKeySet(psGTIF, GeogAngularUnitsGeoKey, TYPE_SHORT, 1,
3282 : Angular_Gon);
3283 4 : else if (EQUAL(osAngUnitName.c_str(), "radian"))
3284 1 : GTIFKeySet(psGTIF, GeogAngularUnitsGeoKey, TYPE_SHORT, 1,
3285 : Angular_Radian);
3286 : // else if (EQUAL(osAngUnitName.c_str(), "microradian"))
3287 : // GTIFKeySet(psGTIF, GeogAngularUnitsGeoKey, TYPE_SHORT, 1,
3288 : // 9109);
3289 : else
3290 : {
3291 : // GeogCitationGeoKey may be rewritten if the gcs is user defined.
3292 3 : oMapAsciiKeys[GeogCitationGeoKey] = osAngUnitName;
3293 3 : GTIFKeySet(psGTIF, GeogAngularUnitSizeGeoKey, TYPE_DOUBLE, 1,
3294 : dfAngUnitValue);
3295 : }
3296 : }
3297 :
3298 : /* -------------------------------------------------------------------- */
3299 : /* Try to write a citation from the main coordinate system */
3300 : /* name. */
3301 : /* -------------------------------------------------------------------- */
3302 6665 : if (poSRS->GetName() != nullptr &&
3303 3277 : ((poSRS->IsProjected() &&
3304 1025 : (nPCS == KvUserDefined || eVersion == GEOTIFF_VERSION_1_0)) ||
3305 2193 : poSRS->IsCompound() || poSRS->IsLocal() ||
3306 2178 : (poSRS->IsGeocentric() &&
3307 0 : (nGCS == KvUserDefined || eVersion == GEOTIFF_VERSION_1_0))))
3308 : {
3309 1101 : if (!(bWritePEString && nPCS == 3857))
3310 : {
3311 1100 : oMapAsciiKeys[GTCitationGeoKey] = poSRS->GetName();
3312 : }
3313 : }
3314 :
3315 : /* -------------------------------------------------------------------- */
3316 : /* Try to write a GCS citation. */
3317 : /* -------------------------------------------------------------------- */
3318 3388 : if (nGCS == KvUserDefined || eVersion == GEOTIFF_VERSION_1_0)
3319 : {
3320 3371 : const OGR_SRSNode *poGCS = poSRS->GetAttrNode("GEOGCS");
3321 :
3322 3371 : if (poGCS != nullptr && poGCS->GetChild(0) != nullptr)
3323 : {
3324 3254 : oMapAsciiKeys[GeogCitationGeoKey] = poGCS->GetChild(0)->GetValue();
3325 : }
3326 : }
3327 :
3328 : /* -------------------------------------------------------------------- */
3329 : /* Try to identify the GCS/datum, scanning the EPSG datum file for */
3330 : /* a match. */
3331 : /* -------------------------------------------------------------------- */
3332 3388 : if (nPCS == KvUserDefined)
3333 : {
3334 2903 : if (nGCS == KvUserDefined && poSRS->IsGeographic() &&
3335 540 : std::fabs(poSRS->GetAngularUnits() - CPLAtof(SRS_UA_DEGREE_CONV)) <
3336 : 1e-9)
3337 : {
3338 533 : if (nDatum == Datum_North_American_Datum_1927)
3339 0 : nGCS = GCS_NAD27;
3340 533 : else if (nDatum == Datum_North_American_Datum_1983)
3341 0 : nGCS = GCS_NAD83;
3342 533 : else if (nDatum == Datum_WGS84 || nDatum == DatumE_WGS84)
3343 527 : nGCS = GCS_WGS_84;
3344 : }
3345 :
3346 2363 : if (nGCS != KvUserDefined)
3347 : {
3348 2176 : GTIFKeySet(psGTIF, GeographicTypeGeoKey, TYPE_SHORT, 1, nGCS);
3349 : }
3350 187 : else if (nDatum != KvUserDefined)
3351 : {
3352 29 : GTIFKeySet(psGTIF, GeographicTypeGeoKey, TYPE_SHORT, 1,
3353 : KvUserDefined);
3354 29 : GTIFKeySet(psGTIF, GeogGeodeticDatumGeoKey, TYPE_SHORT, 1, nDatum);
3355 : }
3356 158 : else if (nSpheroid != KvUserDefined)
3357 : {
3358 3 : GTIFKeySet(psGTIF, GeographicTypeGeoKey, TYPE_SHORT, 1,
3359 : KvUserDefined);
3360 3 : GTIFKeySet(psGTIF, GeogGeodeticDatumGeoKey, TYPE_SHORT, 1,
3361 : KvUserDefined);
3362 3 : GTIFKeySet(psGTIF, GeogEllipsoidGeoKey, TYPE_SHORT, 1, nSpheroid);
3363 : }
3364 155 : else if (dfSemiMajor != 0.0)
3365 : {
3366 41 : GTIFKeySet(psGTIF, GeographicTypeGeoKey, TYPE_SHORT, 1,
3367 : KvUserDefined);
3368 41 : GTIFKeySet(psGTIF, GeogGeodeticDatumGeoKey, TYPE_SHORT, 1,
3369 : KvUserDefined);
3370 41 : GTIFKeySet(psGTIF, GeogEllipsoidGeoKey, TYPE_SHORT, 1,
3371 : KvUserDefined);
3372 41 : GTIFKeySet(psGTIF, GeogSemiMajorAxisGeoKey, TYPE_DOUBLE, 1,
3373 : dfSemiMajor);
3374 41 : if (dfInvFlattening == 0.0)
3375 14 : GTIFKeySet(psGTIF, GeogSemiMinorAxisGeoKey, TYPE_DOUBLE, 1,
3376 : dfSemiMajor);
3377 : else
3378 27 : GTIFKeySet(psGTIF, GeogInvFlatteningGeoKey, TYPE_DOUBLE, 1,
3379 : dfInvFlattening);
3380 : }
3381 114 : else if (poSRS->GetAttrValue("DATUM") != nullptr &&
3382 114 : strstr(poSRS->GetAttrValue("DATUM"), "unknown") == nullptr &&
3383 0 : strstr(poSRS->GetAttrValue("DATUM"), "unnamed") == nullptr)
3384 :
3385 : {
3386 0 : CPLError(CE_Warning, CPLE_AppDefined,
3387 : "Couldn't translate `%s' to a GeoTIFF datum.",
3388 : poSRS->GetAttrValue("DATUM"));
3389 : }
3390 :
3391 2363 : if (nGCS == KvUserDefined || eVersion == GEOTIFF_VERSION_1_0)
3392 : {
3393 : // Always set InvFlattening if it is available.
3394 : // So that it doesn't need to calculate from SemiMinor.
3395 2354 : if (dfInvFlattening != 0.0)
3396 2222 : GTIFKeySet(psGTIF, GeogInvFlatteningGeoKey, TYPE_DOUBLE, 1,
3397 : dfInvFlattening);
3398 : // Always set SemiMajor to keep the precision and in case of
3399 : // editing.
3400 2354 : if (dfSemiMajor != 0.0)
3401 2240 : GTIFKeySet(psGTIF, GeogSemiMajorAxisGeoKey, TYPE_DOUBLE, 1,
3402 : dfSemiMajor);
3403 :
3404 2541 : if (nGCS == KvUserDefined &&
3405 187 : CPLTestBool(CPLGetConfigOption("GTIFF_ESRI_CITATION", "YES")))
3406 : {
3407 187 : SetGeogCSCitation(psGTIF, oMapAsciiKeys, poSRS,
3408 : osAngUnitName.c_str(), nDatum, nSpheroid);
3409 : }
3410 : }
3411 : }
3412 :
3413 : /* -------------------------------------------------------------------- */
3414 : /* Do we have TOWGS84 parameters? */
3415 : /* -------------------------------------------------------------------- */
3416 : #if !defined(GEO_NORMALIZE_DISABLE_TOWGS84)
3417 3388 : double adfTOWGS84[7] = {0.0};
3418 :
3419 6759 : if ((nGCS == KvUserDefined || eVersion == GEOTIFF_VERSION_1_0) &&
3420 3371 : poSRS->GetTOWGS84(adfTOWGS84) == OGRERR_NONE)
3421 : {
3422 : // If we are writing a SRS with a EPSG code, and that the EPSG code
3423 : // of the current SRS object and the one coming from the EPSG code
3424 : // are the same, then by default, do not write them.
3425 23 : bool bUseReferenceTOWGS84 = false;
3426 23 : const char *pszAuthName = poSRS->GetAuthorityName(nullptr);
3427 23 : const char *pszAuthCode = poSRS->GetAuthorityCode(nullptr);
3428 23 : if (pszAuthName && EQUAL(pszAuthName, "EPSG") && pszAuthCode)
3429 : {
3430 24 : CPLErrorHandlerPusher oErrorHandler(CPLQuietErrorHandler);
3431 24 : OGRSpatialReference oRefSRS;
3432 12 : double adfRefTOWGS84[7] = {0.0};
3433 12 : if (oRefSRS.importFromEPSG(atoi(pszAuthCode)) == OGRERR_NONE)
3434 : {
3435 12 : oRefSRS.AddGuessedTOWGS84();
3436 15 : if (oRefSRS.GetTOWGS84(adfRefTOWGS84) == OGRERR_NONE &&
3437 3 : memcmp(adfRefTOWGS84, adfTOWGS84, sizeof(adfTOWGS84)) == 0)
3438 : {
3439 2 : bUseReferenceTOWGS84 = true;
3440 : }
3441 : }
3442 : }
3443 : const char *pszWriteTOWGS84 =
3444 23 : CPLGetConfigOption("GTIFF_WRITE_TOWGS84", "AUTO");
3445 23 : if ((EQUAL(pszWriteTOWGS84, "YES") || EQUAL(pszWriteTOWGS84, "TRUE") ||
3446 22 : EQUAL(pszWriteTOWGS84, "ON")) ||
3447 22 : (!bUseReferenceTOWGS84 && EQUAL(pszWriteTOWGS84, "AUTO")))
3448 : {
3449 22 : if (adfTOWGS84[3] == 0.0 && adfTOWGS84[4] == 0.0 &&
3450 18 : adfTOWGS84[5] == 0.0 && adfTOWGS84[6] == 0.0)
3451 : {
3452 15 : if (nGCS == GCS_WGS_84 && adfTOWGS84[0] == 0.0 &&
3453 9 : adfTOWGS84[1] == 0.0 && adfTOWGS84[2] == 0.0)
3454 : {
3455 : ; // Do nothing.
3456 : }
3457 : else
3458 6 : GTIFKeySet(psGTIF, GeogTOWGS84GeoKey, TYPE_DOUBLE, 3,
3459 : adfTOWGS84);
3460 : }
3461 : else
3462 7 : GTIFKeySet(psGTIF, GeogTOWGS84GeoKey, TYPE_DOUBLE, 7,
3463 : adfTOWGS84);
3464 : }
3465 : }
3466 : #endif
3467 :
3468 : /* -------------------------------------------------------------------- */
3469 : /* Do we have vertical information to set? */
3470 : /* -------------------------------------------------------------------- */
3471 3388 : if (poSRS->GetAttrValue("COMPD_CS|VERT_CS") != nullptr)
3472 : {
3473 13 : bool bGotVertCSCode = false;
3474 13 : const char *pszVertCSCode = poSRS->GetAuthorityCode("COMPD_CS|VERT_CS");
3475 : const char *pszVertCSAuthName =
3476 13 : poSRS->GetAuthorityName("COMPD_CS|VERT_CS");
3477 13 : if (pszVertCSCode && pszVertCSAuthName && atoi(pszVertCSCode) &&
3478 11 : EQUAL(pszVertCSAuthName, "EPSG"))
3479 : {
3480 10 : bGotVertCSCode = true;
3481 10 : GTIFKeySet(psGTIF, VerticalCSTypeGeoKey, TYPE_SHORT, 1,
3482 : atoi(pszVertCSCode));
3483 : }
3484 3 : else if (eVersion >= GEOTIFF_VERSION_1_1)
3485 : {
3486 3 : GTIFKeySet(psGTIF, VerticalCSTypeGeoKey, TYPE_SHORT, 1,
3487 : KvUserDefined);
3488 : }
3489 :
3490 13 : if (eVersion == GEOTIFF_VERSION_1_0 || !bGotVertCSCode)
3491 : {
3492 0 : oMapAsciiKeys[VerticalCitationGeoKey] =
3493 3 : poSRS->GetAttrValue("COMPD_CS|VERT_CS");
3494 :
3495 : const char *pszVertDatumCode =
3496 3 : poSRS->GetAuthorityCode("COMPD_CS|VERT_CS|VERT_DATUM");
3497 : const char *pszVertDatumAuthName =
3498 3 : poSRS->GetAuthorityName("COMPD_CS|VERT_CS|VERT_DATUM");
3499 3 : if (pszVertDatumCode && pszVertDatumAuthName &&
3500 2 : atoi(pszVertDatumCode) && EQUAL(pszVertDatumAuthName, "EPSG"))
3501 : {
3502 1 : GTIFKeySet(psGTIF, VerticalDatumGeoKey, TYPE_SHORT, 1,
3503 : atoi(pszVertDatumCode));
3504 : }
3505 2 : else if (eVersion >= GEOTIFF_VERSION_1_1)
3506 : {
3507 2 : GTIFKeySet(psGTIF, VerticalDatumGeoKey, TYPE_SHORT, 1,
3508 : KvUserDefined);
3509 : }
3510 :
3511 : const char *pszVertUnitCode =
3512 3 : poSRS->GetAuthorityCode("COMPD_CS|VERT_CS|UNIT");
3513 : const char *pszVertUnitAuthName =
3514 3 : poSRS->GetAuthorityName("COMPD_CS|VERT_CS|UNIT");
3515 3 : if (pszVertUnitCode && pszVertUnitAuthName &&
3516 3 : atoi(pszVertUnitCode) && EQUAL(pszVertUnitAuthName, "EPSG"))
3517 : {
3518 3 : GTIFKeySet(psGTIF, VerticalUnitsGeoKey, TYPE_SHORT, 1,
3519 : atoi(pszVertUnitCode));
3520 : }
3521 0 : else if (eVersion >= GEOTIFF_VERSION_1_1)
3522 : {
3523 0 : GTIFKeySet(psGTIF, VerticalUnitsGeoKey, TYPE_SHORT, 1,
3524 : KvUserDefined);
3525 : }
3526 : }
3527 : }
3528 3375 : else if (eVersion >= GEOTIFF_VERSION_1_1 && nVerticalCSKeyValue != 0)
3529 : {
3530 3 : GTIFKeySet(psGTIF, VerticalCSTypeGeoKey, TYPE_SHORT, 1,
3531 : nVerticalCSKeyValue);
3532 : }
3533 :
3534 3388 : const double dfCoordinateEpoch = poSRS->GetCoordinateEpoch();
3535 3388 : if (dfCoordinateEpoch > 0)
3536 : {
3537 2 : GTIFKeySet(psGTIF, CoordinateEpochGeoKey, TYPE_DOUBLE, 1,
3538 : dfCoordinateEpoch);
3539 : }
3540 :
3541 : /* -------------------------------------------------------------------- */
3542 : /* Write all ascii keys */
3543 : /* -------------------------------------------------------------------- */
3544 7755 : for (const auto &oIter : oMapAsciiKeys)
3545 : {
3546 4367 : GTIFKeySet(psGTIF, oIter.first, TYPE_ASCII, 0, oIter.second.c_str());
3547 : }
3548 :
3549 6776 : return TRUE;
3550 : }
3551 :
3552 : /************************************************************************/
3553 : /* GTIFWktFromMemBuf() */
3554 : /************************************************************************/
3555 :
3556 0 : CPLErr GTIFWktFromMemBuf(int nSize, unsigned char *pabyBuffer, char **ppszWKT,
3557 : double *padfGeoTransform, int *pnGCPCount,
3558 : GDAL_GCP **ppasGCPList)
3559 : {
3560 0 : OGRSpatialReferenceH hSRS = nullptr;
3561 0 : if (ppszWKT)
3562 0 : *ppszWKT = nullptr;
3563 : CPLErr eErr =
3564 0 : GTIFWktFromMemBufEx(nSize, pabyBuffer, &hSRS, padfGeoTransform,
3565 : pnGCPCount, ppasGCPList, nullptr, nullptr);
3566 0 : if (eErr == CE_None)
3567 : {
3568 0 : if (hSRS && ppszWKT)
3569 : {
3570 0 : OSRExportToWkt(hSRS, ppszWKT);
3571 : }
3572 : }
3573 0 : OSRDestroySpatialReference(hSRS);
3574 0 : return eErr;
3575 : }
3576 :
3577 581 : CPLErr GTIFWktFromMemBufEx(int nSize, unsigned char *pabyBuffer,
3578 : OGRSpatialReferenceH *phSRS,
3579 : double *padfGeoTransform, int *pnGCPCount,
3580 : GDAL_GCP **ppasGCPList, int *pbPixelIsPoint,
3581 : char ***ppapszRPCMD)
3582 :
3583 : {
3584 : const std::string osFilename(
3585 1162 : VSIMemGenerateHiddenFilename("wkt_from_mem_buf.tif"));
3586 :
3587 : /* -------------------------------------------------------------------- */
3588 : /* Initialization of libtiff and libgeotiff. */
3589 : /* -------------------------------------------------------------------- */
3590 581 : GTiffOneTimeInit(); // For RPC tag.
3591 581 : LibgeotiffOneTimeInit();
3592 :
3593 : /* -------------------------------------------------------------------- */
3594 : /* Create a memory file from the buffer. */
3595 : /* -------------------------------------------------------------------- */
3596 : VSILFILE *fp =
3597 581 : VSIFileFromMemBuffer(osFilename.c_str(), pabyBuffer, nSize, FALSE);
3598 581 : if (fp == nullptr)
3599 0 : return CE_Failure;
3600 :
3601 : /* -------------------------------------------------------------------- */
3602 : /* Initialize access to the memory geotiff structure. */
3603 : /* -------------------------------------------------------------------- */
3604 581 : TIFF *hTIFF = VSI_TIFFOpen(osFilename.c_str(), "rc", fp);
3605 :
3606 581 : if (hTIFF == nullptr)
3607 : {
3608 0 : CPLError(CE_Failure, CPLE_AppDefined,
3609 : "TIFF/GeoTIFF structure is corrupt.");
3610 0 : VSIUnlink(osFilename.c_str());
3611 0 : CPL_IGNORE_RET_VAL(VSIFCloseL(fp));
3612 0 : return CE_Failure;
3613 : }
3614 :
3615 : /* -------------------------------------------------------------------- */
3616 : /* Get the projection definition. */
3617 : /* -------------------------------------------------------------------- */
3618 581 : bool bPixelIsPoint = false;
3619 581 : bool bPointGeoIgnore = false;
3620 581 : unsigned short nRasterType = 0;
3621 :
3622 581 : GTIF *hGTIF = GTIFNew(hTIFF);
3623 581 : if (hGTIF)
3624 581 : GTIFAttachPROJContext(hGTIF, OSRGetProjTLSContext());
3625 :
3626 1162 : if (hGTIF != nullptr &&
3627 581 : GDALGTIFKeyGetSHORT(hGTIF, GTRasterTypeGeoKey, &nRasterType, 0, 1) ==
3628 1162 : 1 &&
3629 580 : nRasterType == static_cast<unsigned short>(RasterPixelIsPoint))
3630 : {
3631 5 : bPixelIsPoint = true;
3632 : bPointGeoIgnore =
3633 5 : CPLTestBool(CPLGetConfigOption("GTIFF_POINT_GEO_IGNORE", "FALSE"));
3634 : }
3635 581 : if (pbPixelIsPoint)
3636 581 : *pbPixelIsPoint = bPixelIsPoint;
3637 581 : if (ppapszRPCMD)
3638 581 : *ppapszRPCMD = nullptr;
3639 :
3640 581 : if (phSRS)
3641 : {
3642 581 : *phSRS = nullptr;
3643 581 : if (hGTIF != nullptr)
3644 : {
3645 581 : GTIFDefn *psGTIFDefn = GTIFAllocDefn();
3646 581 : if (GTIFGetDefn(hGTIF, psGTIFDefn))
3647 : {
3648 580 : *phSRS = GTIFGetOGISDefnAsOSR(hGTIF, psGTIFDefn);
3649 : }
3650 581 : GTIFFreeDefn(psGTIFDefn);
3651 : }
3652 : }
3653 581 : if (hGTIF)
3654 581 : GTIFFree(hGTIF);
3655 :
3656 : /* -------------------------------------------------------------------- */
3657 : /* Get geotransform or tiepoints. */
3658 : /* -------------------------------------------------------------------- */
3659 581 : double *padfTiePoints = nullptr;
3660 581 : double *padfScale = nullptr;
3661 581 : double *padfMatrix = nullptr;
3662 581 : int16_t nCount = 0;
3663 :
3664 581 : padfGeoTransform[0] = 0.0;
3665 581 : padfGeoTransform[1] = 1.0;
3666 581 : padfGeoTransform[2] = 0.0;
3667 581 : padfGeoTransform[3] = 0.0;
3668 581 : padfGeoTransform[4] = 0.0;
3669 581 : padfGeoTransform[5] = 1.0;
3670 :
3671 581 : *pnGCPCount = 0;
3672 581 : *ppasGCPList = nullptr;
3673 :
3674 1123 : if (TIFFGetField(hTIFF, TIFFTAG_GEOPIXELSCALE, &nCount, &padfScale) &&
3675 542 : nCount >= 2)
3676 : {
3677 542 : padfGeoTransform[1] = padfScale[0];
3678 542 : padfGeoTransform[5] = -std::abs(padfScale[1]);
3679 :
3680 542 : if (TIFFGetField(hTIFF, TIFFTAG_GEOTIEPOINTS, &nCount,
3681 1084 : &padfTiePoints) &&
3682 542 : nCount >= 6)
3683 : {
3684 542 : padfGeoTransform[0] =
3685 542 : padfTiePoints[3] - padfTiePoints[0] * padfGeoTransform[1];
3686 542 : padfGeoTransform[3] =
3687 542 : padfTiePoints[4] - padfTiePoints[1] * padfGeoTransform[5];
3688 :
3689 : // Adjust for pixel is point in transform.
3690 542 : if (bPixelIsPoint && !bPointGeoIgnore)
3691 : {
3692 4 : padfGeoTransform[0] -=
3693 4 : padfGeoTransform[1] * 0.5 + padfGeoTransform[2] * 0.5;
3694 4 : padfGeoTransform[3] -=
3695 4 : padfGeoTransform[4] * 0.5 + padfGeoTransform[5] * 0.5;
3696 : }
3697 : }
3698 : }
3699 39 : else if (TIFFGetField(hTIFF, TIFFTAG_GEOTIEPOINTS, &nCount,
3700 52 : &padfTiePoints) &&
3701 13 : nCount >= 6)
3702 : {
3703 13 : *pnGCPCount = nCount / 6;
3704 13 : *ppasGCPList =
3705 13 : static_cast<GDAL_GCP *>(CPLCalloc(sizeof(GDAL_GCP), *pnGCPCount));
3706 :
3707 48 : for (int iGCP = 0; iGCP < *pnGCPCount; iGCP++)
3708 : {
3709 35 : char szID[32] = {};
3710 35 : GDAL_GCP *psGCP = *ppasGCPList + iGCP;
3711 :
3712 35 : snprintf(szID, sizeof(szID), "%d", iGCP + 1);
3713 35 : psGCP->pszId = CPLStrdup(szID);
3714 35 : psGCP->pszInfo = CPLStrdup("");
3715 35 : psGCP->dfGCPPixel = padfTiePoints[iGCP * 6 + 0];
3716 35 : psGCP->dfGCPLine = padfTiePoints[iGCP * 6 + 1];
3717 35 : psGCP->dfGCPX = padfTiePoints[iGCP * 6 + 3];
3718 35 : psGCP->dfGCPY = padfTiePoints[iGCP * 6 + 4];
3719 35 : psGCP->dfGCPZ = padfTiePoints[iGCP * 6 + 5];
3720 : }
3721 : }
3722 26 : else if (TIFFGetField(hTIFF, TIFFTAG_GEOTRANSMATRIX, &nCount,
3723 33 : &padfMatrix) &&
3724 7 : nCount == 16)
3725 : {
3726 7 : padfGeoTransform[0] = padfMatrix[3];
3727 7 : padfGeoTransform[1] = padfMatrix[0];
3728 7 : padfGeoTransform[2] = padfMatrix[1];
3729 7 : padfGeoTransform[3] = padfMatrix[7];
3730 7 : padfGeoTransform[4] = padfMatrix[4];
3731 7 : padfGeoTransform[5] = padfMatrix[5];
3732 : }
3733 :
3734 : /* -------------------------------------------------------------------- */
3735 : /* Read RPC */
3736 : /* -------------------------------------------------------------------- */
3737 581 : if (ppapszRPCMD != nullptr)
3738 : {
3739 581 : *ppapszRPCMD = GTiffDatasetReadRPCTag(hTIFF);
3740 : }
3741 :
3742 : /* -------------------------------------------------------------------- */
3743 : /* Cleanup. */
3744 : /* -------------------------------------------------------------------- */
3745 581 : XTIFFClose(hTIFF);
3746 581 : CPL_IGNORE_RET_VAL(VSIFCloseL(fp));
3747 :
3748 581 : VSIUnlink(osFilename.c_str());
3749 :
3750 581 : if (phSRS && *phSRS == nullptr)
3751 1 : return CE_Failure;
3752 :
3753 580 : return CE_None;
3754 : }
3755 :
3756 : /************************************************************************/
3757 : /* GTIFMemBufFromWkt() */
3758 : /************************************************************************/
3759 :
3760 0 : CPLErr GTIFMemBufFromWkt(const char *pszWKT, const double *padfGeoTransform,
3761 : int nGCPCount, const GDAL_GCP *pasGCPList, int *pnSize,
3762 : unsigned char **ppabyBuffer)
3763 : {
3764 0 : OGRSpatialReference oSRS;
3765 0 : if (pszWKT != nullptr)
3766 0 : oSRS.importFromWkt(pszWKT);
3767 0 : return GTIFMemBufFromSRS(OGRSpatialReference::ToHandle(&oSRS),
3768 : padfGeoTransform, nGCPCount, pasGCPList, pnSize,
3769 0 : ppabyBuffer, FALSE, nullptr);
3770 : }
3771 :
3772 226 : CPLErr GTIFMemBufFromSRS(OGRSpatialReferenceH hSRS,
3773 : const double *padfGeoTransform, int nGCPCount,
3774 : const GDAL_GCP *pasGCPList, int *pnSize,
3775 : unsigned char **ppabyBuffer, int bPixelIsPoint,
3776 : char **papszRPCMD)
3777 :
3778 : {
3779 : const std::string osFilename(
3780 452 : VSIMemGenerateHiddenFilename("wkt_from_mem_buf.tif"));
3781 :
3782 : /* -------------------------------------------------------------------- */
3783 : /* Initialization of libtiff and libgeotiff. */
3784 : /* -------------------------------------------------------------------- */
3785 226 : GTiffOneTimeInit(); // For RPC tag.
3786 226 : LibgeotiffOneTimeInit();
3787 :
3788 : /* -------------------------------------------------------------------- */
3789 : /* Initialize access to the memory geotiff structure. */
3790 : /* -------------------------------------------------------------------- */
3791 226 : VSILFILE *fpL = VSIFOpenL(osFilename.c_str(), "w");
3792 226 : if (fpL == nullptr)
3793 0 : return CE_Failure;
3794 :
3795 226 : TIFF *hTIFF = VSI_TIFFOpen(osFilename.c_str(), "w", fpL);
3796 :
3797 226 : if (hTIFF == nullptr)
3798 : {
3799 0 : CPLError(CE_Failure, CPLE_AppDefined,
3800 : "TIFF/GeoTIFF structure is corrupt.");
3801 0 : CPL_IGNORE_RET_VAL(VSIFCloseL(fpL));
3802 0 : VSIUnlink(osFilename.c_str());
3803 0 : return CE_Failure;
3804 : }
3805 :
3806 : /* -------------------------------------------------------------------- */
3807 : /* Write some minimal set of image parameters. */
3808 : /* -------------------------------------------------------------------- */
3809 226 : TIFFSetField(hTIFF, TIFFTAG_IMAGEWIDTH, 1);
3810 226 : TIFFSetField(hTIFF, TIFFTAG_IMAGELENGTH, 1);
3811 226 : TIFFSetField(hTIFF, TIFFTAG_BITSPERSAMPLE, 8);
3812 226 : TIFFSetField(hTIFF, TIFFTAG_SAMPLESPERPIXEL, 1);
3813 226 : TIFFSetField(hTIFF, TIFFTAG_ROWSPERSTRIP, 1);
3814 226 : TIFFSetField(hTIFF, TIFFTAG_PLANARCONFIG, PLANARCONFIG_CONTIG);
3815 226 : TIFFSetField(hTIFF, TIFFTAG_PHOTOMETRIC, PHOTOMETRIC_MINISBLACK);
3816 :
3817 : /* -------------------------------------------------------------------- */
3818 : /* Get the projection definition. */
3819 : /* -------------------------------------------------------------------- */
3820 :
3821 226 : bool bPointGeoIgnore = false;
3822 226 : if (bPixelIsPoint)
3823 : {
3824 : bPointGeoIgnore =
3825 1 : CPLTestBool(CPLGetConfigOption("GTIFF_POINT_GEO_IGNORE", "FALSE"));
3826 : }
3827 :
3828 226 : GTIF *hGTIF = nullptr;
3829 226 : if (hSRS != nullptr || bPixelIsPoint)
3830 : {
3831 226 : hGTIF = GTIFNew(hTIFF);
3832 226 : if (hGTIF)
3833 : {
3834 226 : GTIFAttachPROJContext(hGTIF, OSRGetProjTLSContext());
3835 :
3836 226 : if (hSRS != nullptr)
3837 226 : GTIFSetFromOGISDefnEx(hGTIF, hSRS, GEOTIFF_KEYS_STANDARD,
3838 : GEOTIFF_VERSION_1_0);
3839 :
3840 226 : if (bPixelIsPoint)
3841 : {
3842 1 : GTIFKeySet(hGTIF, GTRasterTypeGeoKey, TYPE_SHORT, 1,
3843 : RasterPixelIsPoint);
3844 : }
3845 :
3846 226 : GTIFWriteKeys(hGTIF);
3847 226 : GTIFFree(hGTIF);
3848 : }
3849 : }
3850 :
3851 : /* -------------------------------------------------------------------- */
3852 : /* Set the geotransform, or GCPs. */
3853 : /* -------------------------------------------------------------------- */
3854 :
3855 28 : if (padfGeoTransform[0] != 0.0 || padfGeoTransform[1] != 1.0 ||
3856 28 : padfGeoTransform[2] != 0.0 || padfGeoTransform[3] != 0.0 ||
3857 254 : padfGeoTransform[4] != 0.0 || std::abs(padfGeoTransform[5]) != 1.0)
3858 : {
3859 :
3860 209 : if (padfGeoTransform[2] == 0.0 && padfGeoTransform[4] == 0.0)
3861 : {
3862 206 : double adfPixelScale[3] = {padfGeoTransform[1],
3863 206 : fabs(padfGeoTransform[5]), 0.0};
3864 :
3865 206 : TIFFSetField(hTIFF, TIFFTAG_GEOPIXELSCALE, 3, adfPixelScale);
3866 :
3867 206 : double adfTiePoints[6] = {
3868 206 : 0.0, 0.0, 0.0, padfGeoTransform[0], padfGeoTransform[3], 0.0};
3869 :
3870 206 : if (bPixelIsPoint && !bPointGeoIgnore)
3871 : {
3872 1 : adfTiePoints[3] +=
3873 1 : padfGeoTransform[1] * 0.5 + padfGeoTransform[2] * 0.5;
3874 1 : adfTiePoints[4] +=
3875 1 : padfGeoTransform[4] * 0.5 + padfGeoTransform[5] * 0.5;
3876 : }
3877 :
3878 206 : TIFFSetField(hTIFF, TIFFTAG_GEOTIEPOINTS, 6, adfTiePoints);
3879 : }
3880 : else
3881 : {
3882 3 : double adfMatrix[16] = {0.0};
3883 :
3884 3 : adfMatrix[0] = padfGeoTransform[1];
3885 3 : adfMatrix[1] = padfGeoTransform[2];
3886 3 : adfMatrix[3] = padfGeoTransform[0];
3887 3 : adfMatrix[4] = padfGeoTransform[4];
3888 3 : adfMatrix[5] = padfGeoTransform[5];
3889 3 : adfMatrix[7] = padfGeoTransform[3];
3890 3 : adfMatrix[15] = 1.0;
3891 :
3892 3 : if (bPixelIsPoint && !bPointGeoIgnore)
3893 : {
3894 0 : adfMatrix[3] +=
3895 0 : padfGeoTransform[1] * 0.5 + padfGeoTransform[2] * 0.5;
3896 0 : adfMatrix[7] +=
3897 0 : padfGeoTransform[4] * 0.5 + padfGeoTransform[5] * 0.5;
3898 : }
3899 :
3900 3 : TIFFSetField(hTIFF, TIFFTAG_GEOTRANSMATRIX, 16, adfMatrix);
3901 : }
3902 : }
3903 :
3904 : /* -------------------------------------------------------------------- */
3905 : /* Otherwise write tiepoints if they are available. */
3906 : /* -------------------------------------------------------------------- */
3907 17 : else if (nGCPCount > 0)
3908 : {
3909 : double *padfTiePoints =
3910 5 : static_cast<double *>(CPLMalloc(6 * sizeof(double) * nGCPCount));
3911 :
3912 16 : for (int iGCP = 0; iGCP < nGCPCount; iGCP++)
3913 : {
3914 :
3915 11 : padfTiePoints[iGCP * 6 + 0] = pasGCPList[iGCP].dfGCPPixel;
3916 11 : padfTiePoints[iGCP * 6 + 1] = pasGCPList[iGCP].dfGCPLine;
3917 11 : padfTiePoints[iGCP * 6 + 2] = 0;
3918 11 : padfTiePoints[iGCP * 6 + 3] = pasGCPList[iGCP].dfGCPX;
3919 11 : padfTiePoints[iGCP * 6 + 4] = pasGCPList[iGCP].dfGCPY;
3920 11 : padfTiePoints[iGCP * 6 + 5] = pasGCPList[iGCP].dfGCPZ;
3921 : }
3922 :
3923 5 : TIFFSetField(hTIFF, TIFFTAG_GEOTIEPOINTS, 6 * nGCPCount, padfTiePoints);
3924 5 : CPLFree(padfTiePoints);
3925 : }
3926 :
3927 : /* -------------------------------------------------------------------- */
3928 : /* Write RPC */
3929 : /* -------------------------------------------------------------------- */
3930 226 : if (papszRPCMD != nullptr)
3931 : {
3932 1 : GTiffDatasetWriteRPCTag(hTIFF, papszRPCMD);
3933 : }
3934 :
3935 : /* -------------------------------------------------------------------- */
3936 : /* Cleanup and return the created memory buffer. */
3937 : /* -------------------------------------------------------------------- */
3938 226 : GByte bySmallImage = 0;
3939 :
3940 226 : TIFFWriteEncodedStrip(hTIFF, 0, reinterpret_cast<char *>(&bySmallImage), 1);
3941 226 : TIFFWriteCheck(hTIFF, TIFFIsTiled(hTIFF), "GTIFMemBufFromWkt");
3942 226 : TIFFWriteDirectory(hTIFF);
3943 :
3944 226 : XTIFFClose(hTIFF);
3945 226 : CPL_IGNORE_RET_VAL(VSIFCloseL(fpL));
3946 :
3947 : /* -------------------------------------------------------------------- */
3948 : /* Read back from the memory buffer. It would be preferable */
3949 : /* to be able to "steal" the memory buffer, but there isn't */
3950 : /* currently any support for this. */
3951 : /* -------------------------------------------------------------------- */
3952 226 : GUIntBig nBigLength = 0;
3953 :
3954 226 : *ppabyBuffer = VSIGetMemFileBuffer(osFilename.c_str(), &nBigLength, TRUE);
3955 226 : *pnSize = static_cast<int>(nBigLength);
3956 :
3957 226 : return CE_None;
3958 : }
|
