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