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