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