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