"Fossies" - the Fresh Open Source Software Archive

Member "libgeotiff-1.6.0/geo_normalize.c" (22 Apr 2020, 111876 Bytes) of package /linux/privat/libgeotiff-1.6.0.tar.gz:


As a special service "Fossies" has tried to format the requested source page into HTML format using (guessed) C and C++ source code syntax highlighting (style: standard) with prefixed line numbers and code folding option. Alternatively you can here view or download the uninterpreted source code file. For more information about "geo_normalize.c" see the Fossies "Dox" file reference documentation and the latest Fossies "Diffs" side-by-side code changes report: 1.5.1_vs_1.6.0.

    1 /******************************************************************************
    2  * $Id$
    3  *
    4  * Project:  libgeotiff
    5  * Purpose:  Code to normalize PCS and other composite codes in a GeoTIFF file.
    6  * Author:   Frank Warmerdam, warmerdam@pobox.com
    7  *
    8  ******************************************************************************
    9  * Copyright (c) 1999, Frank Warmerdam
   10  * Copyright (c) 2018, Even Rouault <even.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 <assert.h>
   32 
   33 #include "cpl_serv.h"
   34 #include "geo_tiffp.h"
   35 #include "geovalues.h"
   36 #include "geo_normalize.h"
   37 #include "geo_keyp.h"
   38 
   39 #include "proj.h"
   40 
   41 #ifndef KvUserDefined
   42 #  define KvUserDefined 32767
   43 #endif
   44 
   45 #ifndef M_PI
   46 #  define M_PI 3.14159265358979323846
   47 #endif
   48 
   49 /* EPSG Codes for projection parameters.  Unfortunately, these bear no
   50    relationship to the GeoTIFF codes even though the names are so similar. */
   51 
   52 #define EPSGNatOriginLat         8801
   53 #define EPSGNatOriginLong        8802
   54 #define EPSGNatOriginScaleFactor 8805
   55 #define EPSGFalseEasting         8806
   56 #define EPSGFalseNorthing        8807
   57 #define EPSGProjCenterLat        8811
   58 #define EPSGProjCenterLong       8812
   59 #define EPSGAzimuth              8813
   60 #define EPSGAngleRectifiedToSkewedGrid 8814
   61 #define EPSGInitialLineScaleFactor 8815
   62 #define EPSGProjCenterEasting    8816
   63 #define EPSGProjCenterNorthing   8817
   64 #define EPSGPseudoStdParallelLat 8818
   65 #define EPSGPseudoStdParallelScaleFactor 8819
   66 #define EPSGFalseOriginLat       8821
   67 #define EPSGFalseOriginLong      8822
   68 #define EPSGStdParallel1Lat      8823
   69 #define EPSGStdParallel2Lat      8824
   70 #define EPSGFalseOriginEasting   8826
   71 #define EPSGFalseOriginNorthing  8827
   72 #define EPSGSphericalOriginLat   8828
   73 #define EPSGSphericalOriginLong  8829
   74 #define EPSGInitialLongitude     8830
   75 #define EPSGZoneWidth            8831
   76 #define EPSGLatOfStdParallel     8832
   77 #define EPSGOriginLong           8833
   78 #define EPSGTopocentricOriginLat 8834
   79 #define EPSGTopocentricOriginLong 8835
   80 #define EPSGTopocentricOriginHeight 8836
   81 
   82 #define CT_Ext_Mercator_2SP     -CT_Mercator
   83 
   84 #ifndef CPL_INLINE
   85 #  if (defined(__GNUC__) && !defined(__NO_INLINE__)) || defined(_MSC_VER)
   86 #    define HAS_CPL_INLINE  1
   87 #    define CPL_INLINE __inline
   88 #  elif defined(__SUNPRO_CC)
   89 #    define HAS_CPL_INLINE  1
   90 #    define CPL_INLINE inline
   91 #  else
   92 #    define CPL_INLINE
   93 #  endif
   94 #endif
   95 
   96 #ifndef CPL_UNUSED
   97 #if defined(__GNUC__) && __GNUC__ >= 4
   98 #  define CPL_UNUSED __attribute((__unused__))
   99 #else
  100 /* TODO: add cases for other compilers */
  101 #  define CPL_UNUSED
  102 #endif
  103 #endif
  104 
  105 CPL_INLINE static void CPL_IGNORE_RET_VAL_INT(CPL_UNUSED int unused) {}
  106 
  107 /************************************************************************/
  108 /*                         GTIFKeyGetSSHORT()                           */
  109 /************************************************************************/
  110 
  111 // Geotiff SHORT keys are supposed to be unsigned, but geo_normalize interface
  112 // uses signed short...
  113 static int GTIFKeyGetSSHORT( GTIF *gtif, geokey_t key, short* pnVal )
  114 {
  115     unsigned short sVal;
  116     if( GTIFKeyGetSHORT(gtif, key, &sVal, 0, 1) == 1 )
  117     {
  118         memcpy(pnVal, &sVal, 2);
  119         return 1;
  120     }
  121     return 0;
  122 }
  123 
  124 /************************************************************************/
  125 /*                           GTIFGetPCSInfo()                           */
  126 /************************************************************************/
  127 
  128 int GTIFGetPCSInfoEx( void* ctxIn,
  129                       int nPCSCode, char **ppszEPSGName,
  130                       short *pnProjOp, short *pnUOMLengthCode,
  131                       short *pnGeogCS )
  132 
  133 {
  134     PJ_CONTEXT* ctx = (PJ_CONTEXT*)ctxIn;
  135     int         nDatum;
  136     int         nZone;
  137 
  138     /* Deal with a few well known CRS */
  139     int Proj = GTIFPCSToMapSys( nPCSCode, &nDatum, &nZone );
  140     if ((Proj == MapSys_UTM_North || Proj == MapSys_UTM_South) &&
  141         nDatum != KvUserDefined)
  142     {
  143         const char* pszDatumName = NULL;
  144         switch (nDatum)
  145         {
  146             case GCS_NAD27: pszDatumName = "NAD27"; break;
  147             case GCS_NAD83: pszDatumName = "NAD83"; break;
  148             case GCS_WGS_72: pszDatumName = "WGS 72"; break;
  149             case GCS_WGS_72BE: pszDatumName = "WGS 72BE"; break;
  150             case GCS_WGS_84: pszDatumName = "WGS 84"; break;
  151             default: break;
  152         }
  153 
  154         if (pszDatumName)
  155         {
  156             if (ppszEPSGName)
  157             {
  158                 char szEPSGName[64];
  159                 sprintf(szEPSGName, "%s / UTM zone %d%c",
  160                         pszDatumName, nZone, (Proj == MapSys_UTM_North) ? 'N' : 'S');
  161                 *ppszEPSGName = CPLStrdup(szEPSGName);
  162             }
  163 
  164             if (pnProjOp)
  165                 *pnProjOp = (short) (((Proj == MapSys_UTM_North) ? Proj_UTM_zone_1N - 1 : Proj_UTM_zone_1S - 1) + nZone);
  166 
  167             if (pnUOMLengthCode)
  168                 *pnUOMLengthCode = 9001; /* Linear_Meter */
  169 
  170             if (pnGeogCS)
  171                 *pnGeogCS = (short) nDatum;
  172 
  173             return TRUE;
  174         }
  175     }
  176 
  177     if( nPCSCode == KvUserDefined )
  178         return FALSE;
  179 
  180     {
  181         char szCode[12];
  182         PJ* proj_crs;
  183 
  184         sprintf(szCode, "%d", nPCSCode);
  185         proj_crs = proj_create_from_database(
  186             ctx, "EPSG", szCode, PJ_CATEGORY_CRS, 0, NULL);
  187         if( !proj_crs )
  188         {
  189             return FALSE;
  190         }
  191 
  192         if( proj_get_type(proj_crs) != PJ_TYPE_PROJECTED_CRS )
  193         {
  194             proj_destroy(proj_crs);
  195             return FALSE;
  196         }
  197 
  198         if( ppszEPSGName )
  199         {
  200             const char* pszName = proj_get_name(proj_crs);
  201             if( !pszName )
  202             {
  203                 // shouldn't happen
  204                 proj_destroy(proj_crs);
  205                 return FALSE;
  206             }
  207             *ppszEPSGName = CPLStrdup(pszName);
  208         }
  209 
  210         if( pnProjOp )
  211         {
  212             PJ* conversion = proj_crs_get_coordoperation(
  213                 ctx, proj_crs);
  214             if( !conversion )
  215             {
  216                 // shouldn't happen except out of memory
  217                 proj_destroy(proj_crs);
  218                 return FALSE;
  219             }
  220 
  221             {
  222                 const char* pszConvCode = proj_get_id_code(conversion, 0);
  223                 assert( pszConvCode );
  224                 *pnProjOp = (short) atoi(pszConvCode);
  225             }
  226 
  227             proj_destroy(conversion);
  228         }
  229 
  230         if( pnUOMLengthCode )
  231         {
  232             PJ* coordSys = proj_crs_get_coordinate_system(
  233                 ctx, proj_crs);
  234             if( !coordSys )
  235             {
  236                 // shouldn't happen except out of memory
  237                 proj_destroy(proj_crs);
  238                 return FALSE;
  239             }
  240 
  241             {
  242                 const char* pszUnitCode = NULL;
  243                 if( !proj_cs_get_axis_info(
  244                     ctx, coordSys, 0,
  245                     NULL, /* name */
  246                     NULL, /* abbreviation*/
  247                     NULL, /* direction */
  248                     NULL, /* conversion factor */
  249                     NULL, /* unit name */
  250                     NULL, /* unit auth name (should be EPSG) */
  251                     &pszUnitCode) || pszUnitCode == NULL )
  252                 {
  253                     proj_destroy(coordSys);
  254                     return FALSE;
  255                 }
  256                 *pnUOMLengthCode = (short) atoi(pszUnitCode);
  257                 proj_destroy(coordSys);
  258             }
  259         }
  260 
  261         if( pnGeogCS )
  262         {
  263             PJ* geod_crs = proj_crs_get_geodetic_crs(ctx, proj_crs);
  264             if( !geod_crs )
  265             {
  266                 // shouldn't happen except out of memory
  267                 proj_destroy(proj_crs);
  268                 return FALSE;
  269             }
  270 
  271             {
  272                 const char* pszGeodCode = proj_get_id_code(geod_crs, 0);
  273                 assert( pszGeodCode );
  274                 *pnGeogCS = (short) atoi(pszGeodCode);
  275             }
  276 
  277             proj_destroy(geod_crs);
  278         }
  279 
  280 
  281         proj_destroy(proj_crs);
  282         return TRUE;
  283     }
  284 }
  285 
  286 
  287 int GTIFGetPCSInfo( int nPCSCode, char **ppszEPSGName,
  288                     short *pnProjOp, short *pnUOMLengthCode,
  289                     short *pnGeogCS )
  290 
  291 {
  292     PJ_CONTEXT* ctx = proj_context_create();
  293     int ret = GTIFGetPCSInfoEx(ctx, nPCSCode, ppszEPSGName, pnProjOp,
  294                                pnUOMLengthCode, pnGeogCS);
  295     proj_context_destroy(ctx);
  296     return ret;
  297 }
  298 
  299 /************************************************************************/
  300 /*                           GTIFAngleToDD()                            */
  301 /*                                                                      */
  302 /*      Convert a numeric angle to decimal degrees.                     */
  303 /************************************************************************/
  304 
  305 double GTIFAngleToDD( double dfAngle, int nUOMAngle )
  306 
  307 {
  308     if( nUOMAngle == 9110 )     /* DDD.MMSSsss */
  309     {
  310         if( dfAngle > -999.9 && dfAngle < 999.9 )
  311         {
  312             char    szAngleString[32];
  313 
  314             sprintf( szAngleString, "%12.7f", dfAngle );
  315             dfAngle = GTIFAngleStringToDD( szAngleString, nUOMAngle );
  316         }
  317     }
  318     else if ( nUOMAngle != KvUserDefined )
  319     {
  320         double      dfInDegrees = 1.0;
  321 
  322         GTIFGetUOMAngleInfo( nUOMAngle, NULL, &dfInDegrees );
  323         dfAngle = dfAngle * dfInDegrees;
  324     }
  325 
  326     return dfAngle;
  327 }
  328 
  329 /************************************************************************/
  330 /*                        GTIFAngleStringToDD()                         */
  331 /*                                                                      */
  332 /*      Convert an angle in the specified units to decimal degrees.     */
  333 /************************************************************************/
  334 
  335 double GTIFAngleStringToDD( const char * pszAngle, int nUOMAngle )
  336 
  337 {
  338     double  dfAngle;
  339 
  340     if( nUOMAngle == 9110 )     /* DDD.MMSSsss */
  341     {
  342         char    *pszDecimal;
  343 
  344         dfAngle = ABS(atoi(pszAngle));
  345         pszDecimal = strchr(pszAngle,'.');
  346         if( pszDecimal != NULL && strlen(pszDecimal) > 1 )
  347         {
  348             char    szMinutes[3];
  349             char    szSeconds[64];
  350 
  351             szMinutes[0] = pszDecimal[1];
  352             if( pszDecimal[2] >= '0' && pszDecimal[2] <= '9' )
  353                 szMinutes[1] = pszDecimal[2];
  354             else
  355                 szMinutes[1] = '0';
  356 
  357             szMinutes[2] = '\0';
  358             dfAngle += atoi(szMinutes) / 60.0;
  359 
  360             if( strlen(pszDecimal) > 3 )
  361             {
  362                 szSeconds[0] = pszDecimal[3];
  363                 if( pszDecimal[4] >= '0' && pszDecimal[4] <= '9' )
  364                 {
  365                     szSeconds[1] = pszDecimal[4];
  366                     szSeconds[2] = '.';
  367                     strncpy( szSeconds+3, pszDecimal + 5, sizeof(szSeconds) - 3 );
  368                     szSeconds[sizeof(szSeconds) - 1] = 0;
  369                 }
  370                 else
  371                 {
  372                     szSeconds[1] = '0';
  373                     szSeconds[2] = '\0';
  374                 }
  375                 dfAngle += GTIFAtof(szSeconds) / 3600.0;
  376             }
  377         }
  378 
  379         if( pszAngle[0] == '-' )
  380             dfAngle *= -1;
  381     }
  382     else if( nUOMAngle == 9105 || nUOMAngle == 9106 )   /* grad */
  383     {
  384         dfAngle = 180 * (GTIFAtof(pszAngle ) / 200);
  385     }
  386     else if( nUOMAngle == 9101 )            /* radians */
  387     {
  388         dfAngle = 180 * (GTIFAtof(pszAngle ) / M_PI);
  389     }
  390     else if( nUOMAngle == 9103 )            /* arc-minute */
  391     {
  392         dfAngle = GTIFAtof(pszAngle) / 60;
  393     }
  394     else if( nUOMAngle == 9104 )            /* arc-second */
  395     {
  396         dfAngle = GTIFAtof(pszAngle) / 3600;
  397     }
  398     else /* decimal degrees ... some cases missing but seemingly never used */
  399     {
  400         CPLAssert( nUOMAngle == 9102 || nUOMAngle == KvUserDefined
  401                    || nUOMAngle == 0 );
  402 
  403         dfAngle = GTIFAtof(pszAngle );
  404     }
  405 
  406     return dfAngle;
  407 }
  408 
  409 /************************************************************************/
  410 /*                           GTIFGetGCSInfo()                           */
  411 /*                                                                      */
  412 /*      Fetch the datum, and prime meridian related to a particular     */
  413 /*      GCS.                                                            */
  414 /************************************************************************/
  415 
  416 int GTIFGetGCSInfoEx( void* ctxIn,
  417                       int nGCSCode, char ** ppszName,
  418                       short * pnDatum, short * pnPM, short *pnUOMAngle )
  419 
  420 {
  421     PJ_CONTEXT* ctx = (PJ_CONTEXT*)ctxIn;
  422     int     nDatum=0, nPM, nUOMAngle;
  423 
  424 /* -------------------------------------------------------------------- */
  425 /*      Handle some "well known" GCS codes directly                     */
  426 /* -------------------------------------------------------------------- */
  427     const char * pszName = NULL;
  428     nPM = PM_Greenwich;
  429     nUOMAngle = Angular_DMS_Hemisphere;
  430     if( nGCSCode == GCS_NAD27 )
  431     {
  432         nDatum = Datum_North_American_Datum_1927;
  433         pszName = "NAD27";
  434     }
  435     else if( nGCSCode == GCS_NAD83 )
  436     {
  437         nDatum = Datum_North_American_Datum_1983;
  438         pszName = "NAD83";
  439     }
  440     else if( nGCSCode == GCS_WGS_84 )
  441     {
  442         nDatum = Datum_WGS84;
  443         pszName = "WGS 84";
  444     }
  445     else if( nGCSCode == GCS_WGS_72 )
  446     {
  447         nDatum = Datum_WGS72;
  448         pszName = "WGS 72";
  449     }
  450     else if ( nGCSCode == KvUserDefined )
  451     {
  452         return FALSE;
  453     }
  454 
  455     if (pszName != NULL)
  456     {
  457         if( ppszName != NULL )
  458             *ppszName = CPLStrdup( pszName );
  459         if( pnDatum != NULL )
  460             *pnDatum = (short) nDatum;
  461         if( pnPM != NULL )
  462             *pnPM = (short) nPM;
  463         if( pnUOMAngle != NULL )
  464             *pnUOMAngle = (short) nUOMAngle;
  465 
  466         return TRUE;
  467     }
  468 
  469 
  470     if( nGCSCode == KvUserDefined )
  471         return FALSE;
  472 
  473 /* -------------------------------------------------------------------- */
  474 /*      Search the database.                                            */
  475 /* -------------------------------------------------------------------- */
  476 
  477     {
  478         char szCode[12];
  479         PJ* geod_crs;
  480 
  481         sprintf(szCode, "%d", nGCSCode);
  482         geod_crs = proj_create_from_database(
  483             ctx, "EPSG", szCode, PJ_CATEGORY_CRS, 0, NULL);
  484         if( !geod_crs )
  485         {
  486             return FALSE;
  487         }
  488 
  489         {
  490             int objType = proj_get_type(geod_crs);
  491             if( objType != PJ_TYPE_GEODETIC_CRS &&
  492                 objType != PJ_TYPE_GEOCENTRIC_CRS &&
  493                 objType != PJ_TYPE_GEOGRAPHIC_2D_CRS &&
  494                 objType != PJ_TYPE_GEOGRAPHIC_3D_CRS )
  495             {
  496                 proj_destroy(geod_crs);
  497                 return FALSE;
  498             }
  499         }
  500 
  501         if( ppszName )
  502         {
  503             pszName = proj_get_name(geod_crs);
  504             if( !pszName )
  505             {
  506                 // shouldn't happen
  507                 proj_destroy(geod_crs);
  508                 return FALSE;
  509             }
  510             *ppszName = CPLStrdup(pszName);
  511         }
  512 
  513         if( pnDatum )
  514         {
  515             PJ* datum = proj_crs_get_datum(ctx, geod_crs);
  516             if( !datum )
  517             {
  518                 proj_destroy(geod_crs);
  519                 return FALSE;
  520             }
  521 
  522             {
  523                 const char* pszDatumCode = proj_get_id_code(datum, 0);
  524                 assert( pszDatumCode );
  525                 *pnDatum = (short) atoi(pszDatumCode);
  526             }
  527 
  528             proj_destroy(datum);
  529         }
  530 
  531         if( pnPM )
  532         {
  533             PJ* pm = proj_get_prime_meridian(ctx, geod_crs);
  534             if( !pm )
  535             {
  536                 proj_destroy(geod_crs);
  537                 return FALSE;
  538             }
  539 
  540             {
  541                 const char* pszPMCode = proj_get_id_code(pm, 0);
  542                 assert( pszPMCode );
  543                 *pnPM = (short) atoi(pszPMCode);
  544             }
  545 
  546             proj_destroy(pm);
  547         }
  548 
  549         if( pnUOMAngle )
  550         {
  551             PJ* coordSys = proj_crs_get_coordinate_system(
  552                 ctx, geod_crs);
  553             if( !coordSys )
  554             {
  555                 // shouldn't happen except out of memory
  556                 proj_destroy(geod_crs);
  557                 return FALSE;
  558             }
  559 
  560             {
  561                 const char* pszUnitCode = NULL;
  562                 if( !proj_cs_get_axis_info(
  563                     ctx, coordSys, 0,
  564                     NULL, /* name */
  565                     NULL, /* abbreviation*/
  566                     NULL, /* direction */
  567                     NULL, /* conversion factor */
  568                     NULL, /* unit name */
  569                     NULL, /* unit auth name (should be EPSG) */
  570                     &pszUnitCode) || pszUnitCode == NULL )
  571                 {
  572                     proj_destroy(coordSys);
  573                     return FALSE;
  574                 }
  575                 *pnUOMAngle = (short) atoi(pszUnitCode);
  576                 proj_destroy(coordSys);
  577             }
  578         }
  579 
  580         proj_destroy(geod_crs);
  581         return TRUE;
  582     }
  583 }
  584 
  585 int GTIFGetGCSInfo( int nGCSCode, char ** ppszName,
  586                     short * pnDatum, short * pnPM, short *pnUOMAngle )
  587 
  588 {
  589     PJ_CONTEXT* ctx = proj_context_create();
  590     int ret = GTIFGetGCSInfoEx(ctx, nGCSCode, ppszName, pnDatum,
  591                                pnPM, pnUOMAngle);
  592     proj_context_destroy(ctx);
  593     return ret;
  594 }
  595 
  596 /************************************************************************/
  597 /*                        GTIFGetEllipsoidInfo()                        */
  598 /*                                                                      */
  599 /*      Fetch info about an ellipsoid.  Axes are always returned in     */
  600 /*      meters.  SemiMajor computed based on inverse flattening         */
  601 /*      where that is provided.                                         */
  602 /************************************************************************/
  603 
  604 int GTIFGetEllipsoidInfoEx( void* ctxIn,
  605                             int nEllipseCode, char ** ppszName,
  606                             double * pdfSemiMajor, double * pdfSemiMinor )
  607 
  608 {
  609     PJ_CONTEXT* ctx = (PJ_CONTEXT*)ctxIn;
  610 /* -------------------------------------------------------------------- */
  611 /*      Try some well known ellipsoids.                                 */
  612 /* -------------------------------------------------------------------- */
  613     double  dfSemiMajor=0.0;
  614     double     dfInvFlattening=0.0, dfSemiMinor=0.0;
  615     const char *pszName = NULL;
  616 
  617     if( nEllipseCode == Ellipse_Clarke_1866 )
  618     {
  619         pszName = "Clarke 1866";
  620         dfSemiMajor = 6378206.4;
  621         dfSemiMinor = 6356583.8;
  622         dfInvFlattening = 0.0;
  623     }
  624     else if( nEllipseCode == Ellipse_GRS_1980 )
  625     {
  626         pszName = "GRS 1980";
  627         dfSemiMajor = 6378137.0;
  628         dfSemiMinor = 0.0;
  629         dfInvFlattening = 298.257222101;
  630     }
  631     else if( nEllipseCode == Ellipse_WGS_84 )
  632     {
  633         pszName = "WGS 84";
  634         dfSemiMajor = 6378137.0;
  635         dfSemiMinor = 0.0;
  636         dfInvFlattening = 298.257223563;
  637     }
  638     else if( nEllipseCode == 7043 )
  639     {
  640         pszName = "WGS 72";
  641         dfSemiMajor = 6378135.0;
  642         dfSemiMinor = 0.0;
  643         dfInvFlattening = 298.26;
  644     }
  645 
  646     if (pszName != NULL)
  647     {
  648         if( dfSemiMinor == 0.0 )
  649             dfSemiMinor = dfSemiMajor * (1 - 1.0/dfInvFlattening);
  650 
  651         if( pdfSemiMinor != NULL )
  652             *pdfSemiMinor = dfSemiMinor;
  653         if( pdfSemiMajor != NULL )
  654             *pdfSemiMajor = dfSemiMajor;
  655         if( ppszName != NULL )
  656             *ppszName = CPLStrdup( pszName );
  657 
  658         return TRUE;
  659     }
  660 
  661     if( nEllipseCode == KvUserDefined )
  662         return FALSE;
  663 
  664 /* -------------------------------------------------------------------- */
  665 /*      Search the database.                                            */
  666 /* -------------------------------------------------------------------- */
  667     {
  668         char szCode[12];
  669         PJ* ellipsoid;
  670 
  671         sprintf(szCode, "%d", nEllipseCode);
  672         ellipsoid = proj_create_from_database(
  673             ctx, "EPSG", szCode, PJ_CATEGORY_ELLIPSOID, 0, NULL);
  674         if( !ellipsoid )
  675         {
  676             return FALSE;
  677         }
  678 
  679         if( ppszName )
  680         {
  681             pszName = proj_get_name(ellipsoid);
  682             if( !pszName )
  683             {
  684                 // shouldn't happen
  685                 proj_destroy(ellipsoid);
  686                 return FALSE;
  687             }
  688             *ppszName = CPLStrdup(pszName);
  689         }
  690 
  691         proj_ellipsoid_get_parameters(
  692             ctx, ellipsoid, pdfSemiMajor, pdfSemiMinor, NULL, NULL);
  693 
  694         proj_destroy(ellipsoid);
  695 
  696         return TRUE;
  697     }
  698 }
  699 
  700 int GTIFGetEllipsoidInfo( int nEllipseCode, char ** ppszName,
  701                           double * pdfSemiMajor, double * pdfSemiMinor )
  702 
  703 {
  704     PJ_CONTEXT* ctx = proj_context_create();
  705     int ret = GTIFGetEllipsoidInfoEx(ctx, nEllipseCode, ppszName, pdfSemiMajor,
  706                                      pdfSemiMinor);
  707     proj_context_destroy(ctx);
  708     return ret;
  709 }
  710 
  711 /************************************************************************/
  712 /*                           GTIFGetPMInfo()                            */
  713 /*                                                                      */
  714 /*      Get the offset between a given prime meridian and Greenwich     */
  715 /*      in degrees.                                                     */
  716 /************************************************************************/
  717 
  718 int GTIFGetPMInfoEx( void* ctxIn,
  719                      int nPMCode, char ** ppszName, double *pdfOffset )
  720 
  721 {
  722     PJ_CONTEXT* ctx = (PJ_CONTEXT*)ctxIn;
  723 
  724 /* -------------------------------------------------------------------- */
  725 /*      Use a special short cut for Greenwich, since it is so common.   */
  726 /* -------------------------------------------------------------------- */
  727     if( nPMCode == PM_Greenwich )
  728     {
  729         if( pdfOffset != NULL )
  730             *pdfOffset = 0.0;
  731         if( ppszName != NULL )
  732             *ppszName = CPLStrdup( "Greenwich" );
  733         return TRUE;
  734     }
  735 
  736 
  737     if( nPMCode == KvUserDefined )
  738         return FALSE;
  739 
  740 /* -------------------------------------------------------------------- */
  741 /*      Search the database.                                            */
  742 /* -------------------------------------------------------------------- */
  743     {
  744         char szCode[12];
  745         PJ* pm;
  746 
  747         sprintf(szCode, "%d", nPMCode);
  748         pm = proj_create_from_database(
  749             ctx, "EPSG", szCode, PJ_CATEGORY_PRIME_MERIDIAN, 0, NULL);
  750         if( !pm )
  751         {
  752             return FALSE;
  753         }
  754 
  755         if( ppszName )
  756         {
  757             const char* pszName = proj_get_name(pm);
  758             if( !pszName )
  759             {
  760                 // shouldn't happen
  761                 proj_destroy(pm);
  762                 return FALSE;
  763             }
  764             *ppszName = CPLStrdup(pszName);
  765         }
  766 
  767         if( pdfOffset )
  768         {
  769             double conv_factor = 0;
  770             proj_prime_meridian_get_parameters(
  771                 ctx, pm, pdfOffset, &conv_factor, NULL);
  772             *pdfOffset *= conv_factor * 180.0 / M_PI;
  773         }
  774 
  775         proj_destroy(pm);
  776 
  777         return TRUE;
  778     }
  779 }
  780 
  781 int GTIFGetPMInfo( int nPMCode, char ** ppszName, double *pdfOffset )
  782 
  783 {
  784     PJ_CONTEXT* ctx = proj_context_create();
  785     int ret = GTIFGetPMInfoEx(ctx, nPMCode, ppszName, pdfOffset);
  786     proj_context_destroy(ctx);
  787     return ret;
  788 }
  789 
  790 /************************************************************************/
  791 /*                          GTIFGetDatumInfo()                          */
  792 /*                                                                      */
  793 /*      Fetch the ellipsoid, and name for a datum.                      */
  794 /************************************************************************/
  795 
  796 int GTIFGetDatumInfoEx( void* ctxIn,
  797                         int nDatumCode, char ** ppszName, short * pnEllipsoid )
  798 
  799 {
  800     PJ_CONTEXT* ctx = (PJ_CONTEXT*)ctxIn;
  801     const char* pszName = NULL;
  802     int     nEllipsoid = 0;
  803 
  804 /* -------------------------------------------------------------------- */
  805 /*      Handle a few built-in datums.                                   */
  806 /* -------------------------------------------------------------------- */
  807     if( nDatumCode == Datum_North_American_Datum_1927 )
  808     {
  809         nEllipsoid = Ellipse_Clarke_1866;
  810         pszName = "North American Datum 1927";
  811     }
  812     else if( nDatumCode == Datum_North_American_Datum_1983 )
  813     {
  814         nEllipsoid = Ellipse_GRS_1980;
  815         pszName = "North American Datum 1983";
  816     }
  817     else if( nDatumCode == Datum_WGS84 )
  818     {
  819         nEllipsoid = Ellipse_WGS_84;
  820         pszName = "World Geodetic System 1984";
  821     }
  822     else if( nDatumCode == Datum_WGS72 )
  823     {
  824         nEllipsoid = 7043; /* WGS72 */
  825         pszName = "World Geodetic System 1972";
  826     }
  827 
  828     if (pszName != NULL)
  829     {
  830         if( pnEllipsoid != NULL )
  831             *pnEllipsoid = (short) nEllipsoid;
  832 
  833         if( ppszName != NULL )
  834             *ppszName = CPLStrdup( pszName );
  835 
  836         return TRUE;
  837     }
  838 
  839     if( nDatumCode == KvUserDefined )
  840         return FALSE;
  841 
  842 /* -------------------------------------------------------------------- */
  843 /*      Search the database.                                            */
  844 /* -------------------------------------------------------------------- */
  845     {
  846         char szCode[12];
  847         PJ* datum;
  848 
  849         sprintf(szCode, "%d", nDatumCode);
  850         datum = proj_create_from_database(
  851             ctx, "EPSG", szCode, PJ_CATEGORY_DATUM, 0, NULL);
  852         if( !datum )
  853         {
  854             return FALSE;
  855         }
  856 
  857         if( proj_get_type(datum) != PJ_TYPE_GEODETIC_REFERENCE_FRAME )
  858         {
  859             proj_destroy(datum);
  860             return FALSE;
  861         }
  862 
  863         if( ppszName )
  864         {
  865             pszName = proj_get_name(datum);
  866             if( !pszName )
  867             {
  868                 // shouldn't happen
  869                 proj_destroy(datum);
  870                 return FALSE;
  871             }
  872             *ppszName = CPLStrdup(pszName);
  873         }
  874 
  875         if( pnEllipsoid )
  876         {
  877             PJ* ellipsoid = proj_get_ellipsoid(ctx, datum);
  878             if( !ellipsoid )
  879             {
  880                 proj_destroy(datum);
  881                 return FALSE;
  882             }
  883 
  884             {
  885                 const char* pszEllipsoidCode = proj_get_id_code(
  886                     ellipsoid, 0);
  887                 assert( pszEllipsoidCode );
  888                 *pnEllipsoid = (short) atoi(pszEllipsoidCode);
  889             }
  890 
  891             proj_destroy(ellipsoid);
  892         }
  893 
  894         proj_destroy(datum);
  895 
  896         return TRUE;
  897     }
  898 }
  899 
  900 int GTIFGetDatumInfo( int nDatumCode, char ** ppszName, short * pnEllipsoid )
  901 
  902 {
  903     PJ_CONTEXT* ctx = proj_context_create();
  904     int ret = GTIFGetDatumInfoEx(ctx, nDatumCode, ppszName, pnEllipsoid);
  905     proj_context_destroy(ctx);
  906     return ret;
  907 }
  908 
  909 /************************************************************************/
  910 /*                        GTIFGetUOMLengthInfo()                        */
  911 /*                                                                      */
  912 /*      Note: This function should eventually also know how to          */
  913 /*      lookup length aliases in the UOM_LE_ALIAS table.                */
  914 /************************************************************************/
  915 
  916 int GTIFGetUOMLengthInfoEx( void* ctxIn,
  917                             int nUOMLengthCode,
  918                             char **ppszUOMName,
  919                             double * pdfInMeters )
  920 
  921 {
  922     PJ_CONTEXT* ctx = (PJ_CONTEXT*)ctxIn;
  923 /* -------------------------------------------------------------------- */
  924 /*      We short cut meter to save work and avoid failure for missing   */
  925 /*      in the most common cases.                       */
  926 /* -------------------------------------------------------------------- */
  927     if( nUOMLengthCode == 9001 )
  928     {
  929         if( ppszUOMName != NULL )
  930             *ppszUOMName = CPLStrdup( "metre" );
  931         if( pdfInMeters != NULL )
  932             *pdfInMeters = 1.0;
  933 
  934         return TRUE;
  935     }
  936 
  937     if( nUOMLengthCode == 9002 )
  938     {
  939         if( ppszUOMName != NULL )
  940             *ppszUOMName = CPLStrdup( "foot" );
  941         if( pdfInMeters != NULL )
  942             *pdfInMeters = 0.3048;
  943 
  944         return TRUE;
  945     }
  946 
  947     if( nUOMLengthCode == 9003 )
  948     {
  949         if( ppszUOMName != NULL )
  950             *ppszUOMName = CPLStrdup( "US survey foot" );
  951         if( pdfInMeters != NULL )
  952             *pdfInMeters = 12.0 / 39.37;
  953 
  954         return TRUE;
  955     }
  956 
  957     if( nUOMLengthCode == KvUserDefined )
  958         return FALSE;
  959 
  960 /* -------------------------------------------------------------------- */
  961 /*      Search the units database for this unit.  If we don't find      */
  962 /*      it return failure.                                              */
  963 /* -------------------------------------------------------------------- */
  964     {
  965         char szCode[12];
  966         const char* pszName = NULL;
  967 
  968         sprintf(szCode, "%d", nUOMLengthCode);
  969         if( !proj_uom_get_info_from_database(
  970             ctx, "EPSG", szCode, &pszName, pdfInMeters,  NULL) )
  971         {
  972             return FALSE;
  973         }
  974         if( ppszUOMName )
  975         {
  976             *ppszUOMName = CPLStrdup(pszName);
  977         }
  978         return TRUE;
  979     }
  980 }
  981 
  982 int GTIFGetUOMLengthInfo( int nUOMLengthCode,
  983                           char **ppszUOMName,
  984                           double * pdfInMeters )
  985 
  986 {
  987     PJ_CONTEXT* ctx = proj_context_create();
  988     int ret = GTIFGetUOMLengthInfoEx(
  989         ctx, nUOMLengthCode, ppszUOMName, pdfInMeters);
  990     proj_context_destroy(ctx);
  991     return ret;
  992 }
  993 
  994 /************************************************************************/
  995 /*                        GTIFGetUOMAngleInfo()                         */
  996 /************************************************************************/
  997 
  998 int GTIFGetUOMAngleInfoEx( void* ctxIn,
  999                            int nUOMAngleCode,
 1000                            char **ppszUOMName,
 1001                            double * pdfInDegrees )
 1002 
 1003 {
 1004     PJ_CONTEXT* ctx = (PJ_CONTEXT*)ctxIn;
 1005     const char  *pszUOMName = NULL;
 1006     double  dfInDegrees = 1.0;
 1007 
 1008     switch( nUOMAngleCode )
 1009     {
 1010       case 9101:
 1011         pszUOMName = "radian";
 1012         dfInDegrees = 180.0 / M_PI;
 1013         break;
 1014 
 1015       case 9102:
 1016       case 9107:
 1017       case 9108:
 1018       case 9110:
 1019       case 9122:
 1020         pszUOMName = "degree";
 1021         dfInDegrees = 1.0;
 1022         break;
 1023 
 1024       case 9103:
 1025         pszUOMName = "arc-minute";
 1026         dfInDegrees = 1 / 60.0;
 1027         break;
 1028 
 1029       case 9104:
 1030         pszUOMName = "arc-second";
 1031         dfInDegrees = 1 / 3600.0;
 1032         break;
 1033 
 1034       case 9105:
 1035         pszUOMName = "grad";
 1036         dfInDegrees = 180.0 / 200.0;
 1037         break;
 1038 
 1039       case 9106:
 1040         pszUOMName = "gon";
 1041         dfInDegrees = 180.0 / 200.0;
 1042         break;
 1043 
 1044       case 9109:
 1045         pszUOMName = "microradian";
 1046         dfInDegrees = 180.0 / (M_PI * 1000000.0);
 1047         break;
 1048 
 1049       default:
 1050         break;
 1051     }
 1052 
 1053     if (pszUOMName)
 1054     {
 1055         if( ppszUOMName != NULL )
 1056         {
 1057             *ppszUOMName = CPLStrdup( pszUOMName );
 1058         }
 1059 
 1060         if( pdfInDegrees != NULL )
 1061             *pdfInDegrees = dfInDegrees;
 1062 
 1063         return TRUE;
 1064     }
 1065 
 1066 /* -------------------------------------------------------------------- */
 1067 /*      Search the units database for this unit.  If we don't find      */
 1068 /*      it return failure.                                              */
 1069 /* -------------------------------------------------------------------- */
 1070     {
 1071         char szCode[12];
 1072         const char* pszName = NULL;
 1073         double dfConvFactorToRadians = 0;
 1074 
 1075         sprintf(szCode, "%d", nUOMAngleCode);
 1076         if( !proj_uom_get_info_from_database(
 1077             ctx, "EPSG", szCode, &pszName, &dfConvFactorToRadians, NULL) )
 1078         {
 1079             return FALSE;
 1080         }
 1081         if( ppszUOMName )
 1082         {
 1083             *ppszUOMName = CPLStrdup(pszName);
 1084         }
 1085         if( pdfInDegrees )
 1086         {
 1087             *pdfInDegrees = dfConvFactorToRadians * 180.0 / M_PI;
 1088         }
 1089         return TRUE;
 1090     }
 1091 }
 1092 
 1093 int GTIFGetUOMAngleInfo( int nUOMAngleCode,
 1094                          char **ppszUOMName,
 1095                          double * pdfInDegrees )
 1096 
 1097 {
 1098     PJ_CONTEXT* ctx = proj_context_create();
 1099     int ret = GTIFGetUOMAngleInfoEx(
 1100         ctx, nUOMAngleCode, ppszUOMName, pdfInDegrees);
 1101     proj_context_destroy(ctx);
 1102     return ret;
 1103 }
 1104 
 1105 /************************************************************************/
 1106 /*                    EPSGProjMethodToCTProjMethod()                    */
 1107 /*                                                                      */
 1108 /*      Convert between the EPSG enumeration for projection methods,    */
 1109 /*      and the GeoTIFF CT codes.                                       */
 1110 /************************************************************************/
 1111 
 1112 static int EPSGProjMethodToCTProjMethod( int nEPSG, int bReturnExtendedCTCode )
 1113 
 1114 {
 1115     switch( nEPSG )
 1116     {
 1117       case 9801:
 1118         return CT_LambertConfConic_1SP;
 1119 
 1120       case 9802:
 1121         return CT_LambertConfConic_2SP;
 1122 
 1123       case 9803:
 1124         return CT_LambertConfConic_2SP; /* Belgian variant not supported */
 1125 
 1126       case 9804:
 1127         return CT_Mercator;  /* 1SP and 2SP not differentiated */
 1128 
 1129       case 9805:
 1130         if( bReturnExtendedCTCode )
 1131             return CT_Ext_Mercator_2SP;
 1132         else
 1133             return CT_Mercator;  /* 1SP and 2SP not differentiated */
 1134 
 1135       /* Mercator 1SP (Spherical) For EPSG:3785 */
 1136       case 9841:
 1137         return CT_Mercator;  /* 1SP and 2SP not differentiated */
 1138 
 1139       /* Google Mercator For EPSG:3857 */
 1140       case 1024:
 1141         return CT_Mercator;  /* 1SP and 2SP not differentiated */
 1142 
 1143       case 9806:
 1144         return CT_CassiniSoldner;
 1145 
 1146       case 9807:
 1147         return CT_TransverseMercator;
 1148 
 1149       case 9808:
 1150         return CT_TransvMercator_SouthOriented;
 1151 
 1152       case 9809:
 1153         return CT_ObliqueStereographic;
 1154 
 1155       case 9810:
 1156       case 9829: /* variant B not quite the same - not sure how to handle */
 1157         return CT_PolarStereographic;
 1158 
 1159       case 9811:
 1160         return CT_NewZealandMapGrid;
 1161 
 1162       case 9812:
 1163         return CT_ObliqueMercator; /* is hotine actually different? */
 1164 
 1165       case 9813:
 1166         return CT_ObliqueMercator_Laborde;
 1167 
 1168       case 9814:
 1169         return CT_ObliqueMercator_Rosenmund; /* swiss  */
 1170 
 1171       case 9815:
 1172         return CT_HotineObliqueMercatorAzimuthCenter;
 1173 
 1174       case 9816: /* tunesia mining grid has no counterpart */
 1175         return KvUserDefined;
 1176 
 1177       case 9818:
 1178         return CT_Polyconic;
 1179 
 1180       case 9820:
 1181       case 1027:
 1182         return CT_LambertAzimEqualArea;
 1183 
 1184       case 9822:
 1185         return CT_AlbersEqualArea;
 1186 
 1187       case 9834:
 1188         return CT_CylindricalEqualArea;
 1189 
 1190       case 1028:
 1191       case 1029:
 1192       case 9823: /* spherical */
 1193       case 9842: /* elliptical */
 1194         return CT_Equirectangular;
 1195 
 1196       default: /* use the EPSG code for other methods */
 1197         return nEPSG;
 1198     }
 1199 }
 1200 
 1201 /************************************************************************/
 1202 /*                            SetGTParmIds()                            */
 1203 /*                                                                      */
 1204 /*      This is hardcoded logic to set the GeoTIFF parameter            */
 1205 /*      identifiers for all the EPSG supported projections.  As new     */
 1206 /*      projection methods are added, this code will need to be updated */
 1207 /************************************************************************/
 1208 
 1209 static int SetGTParmIds( int nCTProjection,
 1210                          int nEPSGProjMethod,
 1211                          int *panProjParmId,
 1212                          int *panEPSGCodes )
 1213 
 1214 {
 1215     int anWorkingDummy[7];
 1216 
 1217     if( panEPSGCodes == NULL )
 1218         panEPSGCodes = anWorkingDummy;
 1219     if( panProjParmId == NULL )
 1220         panProjParmId = anWorkingDummy;
 1221 
 1222     memset( panEPSGCodes, 0, sizeof(int) * 7 );
 1223 
 1224     /* psDefn->nParms = 7; */
 1225 
 1226     switch( nCTProjection )
 1227     {
 1228       case CT_CassiniSoldner:
 1229       case CT_NewZealandMapGrid:
 1230       case CT_Polyconic:
 1231         panProjParmId[0] = ProjNatOriginLatGeoKey;
 1232         panProjParmId[1] = ProjNatOriginLongGeoKey;
 1233         panProjParmId[5] = ProjFalseEastingGeoKey;
 1234         panProjParmId[6] = ProjFalseNorthingGeoKey;
 1235 
 1236         panEPSGCodes[0] = EPSGNatOriginLat;
 1237         panEPSGCodes[1] = EPSGNatOriginLong;
 1238         panEPSGCodes[5] = EPSGFalseEasting;
 1239         panEPSGCodes[6] = EPSGFalseNorthing;
 1240         return TRUE;
 1241 
 1242       case CT_ObliqueMercator:
 1243       case CT_HotineObliqueMercatorAzimuthCenter:
 1244         panProjParmId[0] = ProjCenterLatGeoKey;
 1245         panProjParmId[1] = ProjCenterLongGeoKey;
 1246         panProjParmId[2] = ProjAzimuthAngleGeoKey;
 1247         panProjParmId[3] = ProjRectifiedGridAngleGeoKey;
 1248         panProjParmId[4] = ProjScaleAtCenterGeoKey;
 1249         panProjParmId[5] = ProjFalseEastingGeoKey;
 1250         panProjParmId[6] = ProjFalseNorthingGeoKey;
 1251 
 1252         panEPSGCodes[0] = EPSGProjCenterLat;
 1253         panEPSGCodes[1] = EPSGProjCenterLong;
 1254         panEPSGCodes[2] = EPSGAzimuth;
 1255         panEPSGCodes[3] = EPSGAngleRectifiedToSkewedGrid;
 1256         panEPSGCodes[4] = EPSGInitialLineScaleFactor;
 1257         panEPSGCodes[5] = EPSGProjCenterEasting; /* EPSG proj method 9812 uses EPSGFalseEasting, but 9815 uses EPSGProjCenterEasting */
 1258         panEPSGCodes[6] = EPSGProjCenterNorthing; /* EPSG proj method 9812 uses EPSGFalseNorthing, but 9815 uses EPSGProjCenterNorthing */
 1259         return TRUE;
 1260 
 1261       case CT_ObliqueMercator_Laborde:
 1262         panProjParmId[0] = ProjCenterLatGeoKey;
 1263         panProjParmId[1] = ProjCenterLongGeoKey;
 1264         panProjParmId[2] = ProjAzimuthAngleGeoKey;
 1265         panProjParmId[4] = ProjScaleAtCenterGeoKey;
 1266         panProjParmId[5] = ProjFalseEastingGeoKey;
 1267         panProjParmId[6] = ProjFalseNorthingGeoKey;
 1268 
 1269         panEPSGCodes[0] = EPSGProjCenterLat;
 1270         panEPSGCodes[1] = EPSGProjCenterLong;
 1271         panEPSGCodes[2] = EPSGAzimuth;
 1272         panEPSGCodes[4] = EPSGInitialLineScaleFactor;
 1273         panEPSGCodes[5] = EPSGFalseEasting;
 1274         panEPSGCodes[6] = EPSGFalseNorthing;
 1275         return TRUE;
 1276 
 1277       case CT_LambertConfConic_1SP:
 1278       case CT_Mercator:
 1279       case CT_ObliqueStereographic:
 1280       case CT_PolarStereographic:
 1281       case CT_TransverseMercator:
 1282       case CT_TransvMercator_SouthOriented:
 1283         panProjParmId[0] = ProjNatOriginLatGeoKey;
 1284         if( nCTProjection == CT_PolarStereographic )
 1285         {
 1286             panProjParmId[1] = ProjStraightVertPoleLongGeoKey;
 1287         }
 1288         else
 1289         {
 1290             panProjParmId[1] = ProjNatOriginLongGeoKey;
 1291         }
 1292         if( nEPSGProjMethod == 9805 ) /* Mercator_2SP */
 1293         {
 1294             panProjParmId[2] = ProjStdParallel1GeoKey;
 1295         }
 1296         panProjParmId[4] = ProjScaleAtNatOriginGeoKey;
 1297         panProjParmId[5] = ProjFalseEastingGeoKey;
 1298         panProjParmId[6] = ProjFalseNorthingGeoKey;
 1299 
 1300         panEPSGCodes[0] = EPSGNatOriginLat;
 1301         panEPSGCodes[1] = EPSGNatOriginLong;
 1302         if( nEPSGProjMethod == 9805 ) /* Mercator_2SP */
 1303         {
 1304             panEPSGCodes[2] = EPSGStdParallel1Lat;
 1305         }
 1306         panEPSGCodes[4] = EPSGNatOriginScaleFactor;
 1307         panEPSGCodes[5] = EPSGFalseEasting;
 1308         panEPSGCodes[6] = EPSGFalseNorthing;
 1309         return TRUE;
 1310 
 1311       case CT_LambertConfConic_2SP:
 1312         panProjParmId[0] = ProjFalseOriginLatGeoKey;
 1313         panProjParmId[1] = ProjFalseOriginLongGeoKey;
 1314         panProjParmId[2] = ProjStdParallel1GeoKey;
 1315         panProjParmId[3] = ProjStdParallel2GeoKey;
 1316         panProjParmId[5] = ProjFalseEastingGeoKey;
 1317         panProjParmId[6] = ProjFalseNorthingGeoKey;
 1318 
 1319         panEPSGCodes[0] = EPSGFalseOriginLat;
 1320         panEPSGCodes[1] = EPSGFalseOriginLong;
 1321         panEPSGCodes[2] = EPSGStdParallel1Lat;
 1322         panEPSGCodes[3] = EPSGStdParallel2Lat;
 1323         panEPSGCodes[5] = EPSGFalseOriginEasting;
 1324         panEPSGCodes[6] = EPSGFalseOriginNorthing;
 1325         return TRUE;
 1326 
 1327       case CT_AlbersEqualArea:
 1328         panProjParmId[0] = ProjStdParallel1GeoKey;
 1329         panProjParmId[1] = ProjStdParallel2GeoKey;
 1330         panProjParmId[2] = ProjNatOriginLatGeoKey;
 1331         panProjParmId[3] = ProjNatOriginLongGeoKey;
 1332         panProjParmId[5] = ProjFalseEastingGeoKey;
 1333         panProjParmId[6] = ProjFalseNorthingGeoKey;
 1334 
 1335         panEPSGCodes[0] = EPSGStdParallel1Lat;
 1336         panEPSGCodes[1] = EPSGStdParallel2Lat;
 1337         panEPSGCodes[2] = EPSGFalseOriginLat;
 1338         panEPSGCodes[3] = EPSGFalseOriginLong;
 1339         panEPSGCodes[5] = EPSGFalseOriginEasting;
 1340         panEPSGCodes[6] = EPSGFalseOriginNorthing;
 1341         return TRUE;
 1342 
 1343       case CT_SwissObliqueCylindrical:
 1344         panProjParmId[0] = ProjCenterLatGeoKey;
 1345         panProjParmId[1] = ProjCenterLongGeoKey;
 1346         panProjParmId[5] = ProjFalseEastingGeoKey;
 1347         panProjParmId[6] = ProjFalseNorthingGeoKey;
 1348 
 1349         /* EPSG codes? */
 1350         return TRUE;
 1351 
 1352       case CT_LambertAzimEqualArea:
 1353         panProjParmId[0] = ProjCenterLatGeoKey;
 1354         panProjParmId[1] = ProjCenterLongGeoKey;
 1355         panProjParmId[5] = ProjFalseEastingGeoKey;
 1356         panProjParmId[6] = ProjFalseNorthingGeoKey;
 1357 
 1358         panEPSGCodes[0] = EPSGNatOriginLat;
 1359         panEPSGCodes[1] = EPSGNatOriginLong;
 1360         panEPSGCodes[5] = EPSGFalseEasting;
 1361         panEPSGCodes[6] = EPSGFalseNorthing;
 1362         return TRUE;
 1363 
 1364       case CT_CylindricalEqualArea:
 1365         panProjParmId[0] = ProjStdParallel1GeoKey;
 1366         panProjParmId[1] = ProjNatOriginLongGeoKey;
 1367         panProjParmId[5] = ProjFalseEastingGeoKey;
 1368         panProjParmId[6] = ProjFalseNorthingGeoKey;
 1369 
 1370         panEPSGCodes[0] = EPSGStdParallel1Lat;
 1371         panEPSGCodes[1] = EPSGFalseOriginLong;
 1372         panEPSGCodes[5] = EPSGFalseOriginEasting;
 1373         panEPSGCodes[6] = EPSGFalseOriginNorthing;
 1374         return TRUE;
 1375 
 1376       case CT_Equirectangular:
 1377         panProjParmId[0] = ProjCenterLatGeoKey;
 1378         panProjParmId[1] = ProjCenterLongGeoKey;
 1379         panProjParmId[2] = ProjStdParallel1GeoKey;
 1380         panProjParmId[5] = ProjFalseEastingGeoKey;
 1381         panProjParmId[6] = ProjFalseNorthingGeoKey;
 1382 
 1383         panEPSGCodes[0] = EPSGNatOriginLat;
 1384         panEPSGCodes[1] = EPSGNatOriginLong;
 1385         panEPSGCodes[2] = EPSGStdParallel1Lat;
 1386         panEPSGCodes[5] = EPSGFalseEasting;
 1387         panEPSGCodes[6] = EPSGFalseNorthing;
 1388         return TRUE;
 1389 
 1390       case CT_Ext_Mercator_2SP:
 1391         panProjParmId[0] = ProjNatOriginLatGeoKey;
 1392         panProjParmId[1] = ProjNatOriginLongGeoKey;
 1393         panProjParmId[2] = ProjStdParallel1GeoKey;
 1394         panProjParmId[5] = ProjFalseEastingGeoKey;
 1395         panProjParmId[6] = ProjFalseNorthingGeoKey;
 1396 
 1397         panEPSGCodes[0] = EPSGNatOriginLat;
 1398         panEPSGCodes[1] = EPSGNatOriginLong;
 1399         panEPSGCodes[2] = EPSGStdParallel1Lat;
 1400         panEPSGCodes[5] = EPSGFalseEasting;
 1401         panEPSGCodes[6] = EPSGFalseNorthing;
 1402         return TRUE;
 1403 
 1404       default:
 1405         return FALSE;
 1406     }
 1407 }
 1408 
 1409 /************************************************************************/
 1410 /*                         GTIFGetProjTRFInfo()                         */
 1411 /*                                                                      */
 1412 /*      Transform a PROJECTION_TRF_CODE into a projection method,       */
 1413 /*      and a set of parameters.  The parameters identify will          */
 1414 /*      depend on the returned method, but they will all have been      */
 1415 /*      normalized into degrees and meters.                             */
 1416 /************************************************************************/
 1417 
 1418 int GTIFGetProjTRFInfoEx( void* ctxIn,
 1419                           int nProjTRFCode,
 1420                           char **ppszProjTRFName,
 1421                           short * pnProjMethod,
 1422                           double * padfProjParms )
 1423 
 1424 {
 1425     PJ_CONTEXT* ctx = (PJ_CONTEXT*)ctxIn;
 1426 
 1427     if ((nProjTRFCode >= Proj_UTM_zone_1N && nProjTRFCode <= Proj_UTM_zone_60N) ||
 1428         (nProjTRFCode >= Proj_UTM_zone_1S && nProjTRFCode <= Proj_UTM_zone_60S))
 1429     {
 1430         int bNorth;
 1431         int nZone;
 1432         if (nProjTRFCode <= Proj_UTM_zone_60N)
 1433         {
 1434             bNorth = TRUE;
 1435             nZone = nProjTRFCode - Proj_UTM_zone_1N + 1;
 1436         }
 1437         else
 1438         {
 1439             bNorth = FALSE;
 1440             nZone = nProjTRFCode - Proj_UTM_zone_1S + 1;
 1441         }
 1442 
 1443         if (ppszProjTRFName)
 1444         {
 1445             char szProjTRFName[64];
 1446             sprintf(szProjTRFName, "UTM zone %d%c",
 1447                     nZone, (bNorth) ? 'N' : 'S');
 1448             *ppszProjTRFName = CPLStrdup(szProjTRFName);
 1449         }
 1450 
 1451         if (pnProjMethod)
 1452             *pnProjMethod = 9807;
 1453 
 1454         if (padfProjParms)
 1455         {
 1456             padfProjParms[0] = 0;
 1457             padfProjParms[1] = -183 + 6 * nZone;
 1458             padfProjParms[2] = 0;
 1459             padfProjParms[3] = 0;
 1460             padfProjParms[4] = 0.9996;
 1461             padfProjParms[5] = 500000;
 1462             padfProjParms[6] = (bNorth) ? 0 : 10000000;
 1463         }
 1464 
 1465         return TRUE;
 1466     }
 1467 
 1468     if( nProjTRFCode == KvUserDefined )
 1469         return FALSE;
 1470 
 1471     {
 1472         int     nProjMethod, i, anEPSGCodes[7];
 1473         double  adfProjParms[7];
 1474         char    szCode[12];
 1475         const char* pszMethodCode = NULL;
 1476         int     nCTProjMethod;
 1477         PJ *transf;
 1478 
 1479         sprintf(szCode, "%d", nProjTRFCode);
 1480         transf = proj_create_from_database(
 1481             ctx, "EPSG", szCode, PJ_CATEGORY_COORDINATE_OPERATION, 0, NULL);
 1482         if( !transf )
 1483         {
 1484             return FALSE;
 1485         }
 1486 
 1487         if( proj_get_type(transf) != PJ_TYPE_CONVERSION )
 1488         {
 1489             proj_destroy(transf);
 1490             return FALSE;
 1491         }
 1492 
 1493         /* Get the projection method code */
 1494         proj_coordoperation_get_method_info(ctx, transf,
 1495                                             NULL, /* method name */
 1496                                             NULL, /* method auth name (should be EPSG) */
 1497                                             &pszMethodCode);
 1498         assert( pszMethodCode );
 1499         nProjMethod = atoi(pszMethodCode);
 1500 
 1501 /* -------------------------------------------------------------------- */
 1502 /*      Initialize a definition of what EPSG codes need to be loaded    */
 1503 /*      into what fields in adfProjParms.                               */
 1504 /* -------------------------------------------------------------------- */
 1505         nCTProjMethod = EPSGProjMethodToCTProjMethod( nProjMethod, TRUE );
 1506         SetGTParmIds( nCTProjMethod, nProjMethod, NULL, anEPSGCodes );
 1507 
 1508 /* -------------------------------------------------------------------- */
 1509 /*      Get the parameters for this projection.                         */
 1510 /* -------------------------------------------------------------------- */
 1511 
 1512         for( i = 0; i < 7; i++ )
 1513         {
 1514             double  dfValue = 0.0;
 1515             double  dfUnitConvFactor = 0.0;
 1516             const char *pszUOMCategory = NULL;
 1517             int     nEPSGCode = anEPSGCodes[i];
 1518             int     iEPSG;
 1519             int     nParamCount;
 1520 
 1521             /* Establish default */
 1522             if( nEPSGCode == EPSGAngleRectifiedToSkewedGrid )
 1523                 adfProjParms[i] = 90.0;
 1524             else if( nEPSGCode == EPSGNatOriginScaleFactor
 1525                     || nEPSGCode == EPSGInitialLineScaleFactor
 1526                     || nEPSGCode == EPSGPseudoStdParallelScaleFactor )
 1527                 adfProjParms[i] = 1.0;
 1528             else
 1529                 adfProjParms[i] = 0.0;
 1530 
 1531             /* If there is no parameter, skip */
 1532             if( nEPSGCode == 0 )
 1533                 continue;
 1534 
 1535             nParamCount = proj_coordoperation_get_param_count(ctx, transf);
 1536 
 1537             /* Find the matching parameter */
 1538             for( iEPSG = 0; iEPSG < nParamCount; iEPSG++ )
 1539             {
 1540                 const char* pszParamCode = NULL;
 1541                 proj_coordoperation_get_param(
 1542                     ctx, transf, iEPSG,
 1543                     NULL, /* name */
 1544                     NULL, /* auth name */
 1545                     &pszParamCode,
 1546                     &dfValue,
 1547                     NULL, /* value (string) */
 1548                     &dfUnitConvFactor, /* unit conv factor */
 1549                     NULL, /* unit name */
 1550                     NULL, /* unit auth name */
 1551                     NULL, /* unit code */
 1552                     &pszUOMCategory /* unit category */);
 1553                 assert(pszParamCode);
 1554                 if( atoi(pszParamCode) == nEPSGCode )
 1555                 {
 1556                     break;
 1557                 }
 1558             }
 1559 
 1560             /* not found, accept the default */
 1561             if( iEPSG == nParamCount )
 1562             {
 1563                 /* for CT_ObliqueMercator try alternate parameter codes first */
 1564                 /* because EPSG proj method 9812 uses EPSGFalseXXXXX, but 9815 uses EPSGProjCenterXXXXX */
 1565                 if ( nCTProjMethod == CT_ObliqueMercator && nEPSGCode == EPSGProjCenterEasting )
 1566                     nEPSGCode = EPSGFalseEasting;
 1567                 else if ( nCTProjMethod == CT_ObliqueMercator && nEPSGCode == EPSGProjCenterNorthing )
 1568                     nEPSGCode = EPSGFalseNorthing;
 1569                 /* for CT_PolarStereographic try alternate parameter codes first */
 1570                 /* because EPSG proj method 9829 uses EPSGLatOfStdParallel instead of EPSGNatOriginLat */
 1571                 /* and EPSGOriginLong instead of EPSGNatOriginLong */
 1572                 else if( nCTProjMethod == CT_PolarStereographic && nEPSGCode == EPSGNatOriginLat )
 1573                     nEPSGCode = EPSGLatOfStdParallel;
 1574                 else if( nCTProjMethod == CT_PolarStereographic && nEPSGCode == EPSGNatOriginLong )
 1575                     nEPSGCode = EPSGOriginLong;
 1576                 else
 1577                     continue;
 1578 
 1579                 for( iEPSG = 0; iEPSG < nParamCount; iEPSG++ )
 1580                 {
 1581                     const char* pszParamCode = NULL;
 1582                     proj_coordoperation_get_param(
 1583                         ctx, transf, iEPSG,
 1584                         NULL, /* name */
 1585                         NULL, /* auth name */
 1586                         &pszParamCode,
 1587                         &dfValue,
 1588                         NULL, /* value (string) */
 1589                         &dfUnitConvFactor, /* unit conv factor */
 1590                         NULL, /* unit name */
 1591                         NULL, /* unit auth name */
 1592                         NULL, /* unit code */
 1593                         &pszUOMCategory /* unit category */);
 1594                     assert(pszParamCode);
 1595                     if( atoi(pszParamCode) == nEPSGCode )
 1596                     {
 1597                         break;
 1598                     }
 1599                 }
 1600 
 1601                 if( iEPSG == nParamCount )
 1602                     continue;
 1603             }
 1604 
 1605             assert(pszUOMCategory);
 1606 
 1607             adfProjParms[i] = dfValue * dfUnitConvFactor;
 1608             if( strcmp(pszUOMCategory, "angular") == 0.0 )
 1609             {
 1610                 /* Convert from radians to degrees */
 1611                 adfProjParms[i] *= 180 / M_PI;
 1612             }
 1613         }
 1614 
 1615 /* -------------------------------------------------------------------- */
 1616 /*      Get the name, if requested.                                     */
 1617 /* -------------------------------------------------------------------- */
 1618         if( ppszProjTRFName != NULL )
 1619         {
 1620             const char* pszName = proj_get_name(transf);
 1621             if( !pszName )
 1622             {
 1623                 // shouldn't happen
 1624                 proj_destroy(transf);
 1625                 return FALSE;
 1626             }
 1627             *ppszProjTRFName = CPLStrdup(pszName);
 1628         }
 1629 
 1630 /* -------------------------------------------------------------------- */
 1631 /*      Transfer requested data into passed variables.                  */
 1632 /* -------------------------------------------------------------------- */
 1633         if( pnProjMethod != NULL )
 1634             *pnProjMethod = (short) nProjMethod;
 1635 
 1636         if( padfProjParms != NULL )
 1637         {
 1638             for( i = 0; i < 7; i++ )
 1639                 padfProjParms[i] = adfProjParms[i];
 1640         }
 1641 
 1642         proj_destroy(transf);
 1643 
 1644         return TRUE;
 1645     }
 1646 }
 1647 
 1648 int GTIFGetProjTRFInfo( /* Conversion code */
 1649                         int nProjTRFCode,
 1650                         char **ppszProjTRFName,
 1651                         short * pnProjMethod,
 1652                         double * padfProjParms )
 1653 {
 1654     PJ_CONTEXT* ctx = proj_context_create();
 1655     int ret = GTIFGetProjTRFInfoEx(
 1656         ctx, nProjTRFCode, ppszProjTRFName, pnProjMethod, padfProjParms);
 1657     proj_context_destroy(ctx);
 1658     return ret;
 1659 }
 1660 
 1661 /************************************************************************/
 1662 /*                         GTIFFetchProjParms()                         */
 1663 /*                                                                      */
 1664 /*      Fetch the projection parameters for a particular projection     */
 1665 /*      from a GeoTIFF file, and fill the GTIFDefn structure out        */
 1666 /*      with them.                                                      */
 1667 /************************************************************************/
 1668 
 1669 static void GTIFFetchProjParms( GTIF * psGTIF, GTIFDefn * psDefn )
 1670 
 1671 {
 1672     double dfNatOriginLong = 0.0, dfNatOriginLat = 0.0, dfRectGridAngle = 0.0;
 1673     double dfFalseEasting = 0.0, dfFalseNorthing = 0.0, dfNatOriginScale = 1.0;
 1674     double dfStdParallel1 = 0.0, dfStdParallel2 = 0.0, dfAzimuth = 0.0;
 1675     int iParm;
 1676     int bHaveSP1, bHaveNOS;
 1677 
 1678 /* -------------------------------------------------------------------- */
 1679 /*      Get the false easting, and northing if available.               */
 1680 /* -------------------------------------------------------------------- */
 1681     if( !GTIFKeyGetDOUBLE(psGTIF, ProjFalseEastingGeoKey, &dfFalseEasting, 0, 1)
 1682         && !GTIFKeyGetDOUBLE(psGTIF, ProjCenterEastingGeoKey,
 1683                        &dfFalseEasting, 0, 1)
 1684         && !GTIFKeyGetDOUBLE(psGTIF, ProjFalseOriginEastingGeoKey,
 1685                        &dfFalseEasting, 0, 1) )
 1686         dfFalseEasting = 0.0;
 1687 
 1688     if( !GTIFKeyGetDOUBLE(psGTIF, ProjFalseNorthingGeoKey, &dfFalseNorthing,0,1)
 1689         && !GTIFKeyGetDOUBLE(psGTIF, ProjCenterNorthingGeoKey,
 1690                        &dfFalseNorthing, 0, 1)
 1691         && !GTIFKeyGetDOUBLE(psGTIF, ProjFalseOriginNorthingGeoKey,
 1692                        &dfFalseNorthing, 0, 1) )
 1693         dfFalseNorthing = 0.0;
 1694 
 1695     switch( psDefn->CTProjection )
 1696     {
 1697 /* -------------------------------------------------------------------- */
 1698       case CT_Stereographic:
 1699 /* -------------------------------------------------------------------- */
 1700         if( GTIFKeyGetDOUBLE(psGTIF, ProjNatOriginLongGeoKey,
 1701                        &dfNatOriginLong, 0, 1 ) == 0
 1702             && GTIFKeyGetDOUBLE(psGTIF, ProjFalseOriginLongGeoKey,
 1703                           &dfNatOriginLong, 0, 1 ) == 0
 1704             && GTIFKeyGetDOUBLE(psGTIF, ProjCenterLongGeoKey,
 1705                           &dfNatOriginLong, 0, 1 ) == 0 )
 1706             dfNatOriginLong = 0.0;
 1707 
 1708         if( GTIFKeyGetDOUBLE(psGTIF, ProjNatOriginLatGeoKey,
 1709                        &dfNatOriginLat, 0, 1 ) == 0
 1710             && GTIFKeyGetDOUBLE(psGTIF, ProjFalseOriginLatGeoKey,
 1711                           &dfNatOriginLat, 0, 1 ) == 0
 1712             && GTIFKeyGetDOUBLE(psGTIF, ProjCenterLatGeoKey,
 1713                           &dfNatOriginLat, 0, 1 ) == 0 )
 1714             dfNatOriginLat = 0.0;
 1715 
 1716         if( GTIFKeyGetDOUBLE(psGTIF, ProjScaleAtNatOriginGeoKey,
 1717                        &dfNatOriginScale, 0, 1 ) == 0 )
 1718             dfNatOriginScale = 1.0;
 1719 
 1720         /* notdef: should transform to decimal degrees at this point */
 1721 
 1722         psDefn->ProjParm[0] = dfNatOriginLat;
 1723         psDefn->ProjParmId[0] = ProjCenterLatGeoKey;
 1724         psDefn->ProjParm[1] = dfNatOriginLong;
 1725         psDefn->ProjParmId[1] = ProjCenterLongGeoKey;
 1726         psDefn->ProjParm[4] = dfNatOriginScale;
 1727         psDefn->ProjParmId[4] = ProjScaleAtNatOriginGeoKey;
 1728         psDefn->ProjParm[5] = dfFalseEasting;
 1729         psDefn->ProjParmId[5] = ProjFalseEastingGeoKey;
 1730         psDefn->ProjParm[6] = dfFalseNorthing;
 1731         psDefn->ProjParmId[6] = ProjFalseNorthingGeoKey;
 1732 
 1733         psDefn->nParms = 7;
 1734         break;
 1735 
 1736 /* -------------------------------------------------------------------- */
 1737       case CT_Mercator:
 1738 /* -------------------------------------------------------------------- */
 1739         if( GTIFKeyGetDOUBLE(psGTIF, ProjNatOriginLongGeoKey,
 1740                        &dfNatOriginLong, 0, 1 ) == 0
 1741             && GTIFKeyGetDOUBLE(psGTIF, ProjFalseOriginLongGeoKey,
 1742                           &dfNatOriginLong, 0, 1 ) == 0
 1743             && GTIFKeyGetDOUBLE(psGTIF, ProjCenterLongGeoKey,
 1744                           &dfNatOriginLong, 0, 1 ) == 0 )
 1745             dfNatOriginLong = 0.0;
 1746 
 1747         if( GTIFKeyGetDOUBLE(psGTIF, ProjNatOriginLatGeoKey,
 1748                        &dfNatOriginLat, 0, 1 ) == 0
 1749             && GTIFKeyGetDOUBLE(psGTIF, ProjFalseOriginLatGeoKey,
 1750                           &dfNatOriginLat, 0, 1 ) == 0
 1751             && GTIFKeyGetDOUBLE(psGTIF, ProjCenterLatGeoKey,
 1752                           &dfNatOriginLat, 0, 1 ) == 0 )
 1753             dfNatOriginLat = 0.0;
 1754 
 1755 
 1756         bHaveSP1 = GTIFKeyGetDOUBLE(psGTIF, ProjStdParallel1GeoKey,
 1757                               &dfStdParallel1, 0, 1 );
 1758 
 1759         bHaveNOS = GTIFKeyGetDOUBLE(psGTIF, ProjScaleAtNatOriginGeoKey,
 1760                               &dfNatOriginScale, 0, 1 );
 1761 
 1762         /* Default scale only if dfStdParallel1 isn't defined either */
 1763         if( !bHaveNOS && !bHaveSP1)
 1764         {
 1765             bHaveNOS = TRUE;
 1766             dfNatOriginScale = 1.0;
 1767         }
 1768 
 1769         /* notdef: should transform to decimal degrees at this point */
 1770 
 1771         psDefn->ProjParm[0] = dfNatOriginLat;
 1772         psDefn->ProjParmId[0] = ProjNatOriginLatGeoKey;
 1773         psDefn->ProjParm[1] = dfNatOriginLong;
 1774         psDefn->ProjParmId[1] = ProjNatOriginLongGeoKey;
 1775         if( bHaveSP1 )
 1776         {
 1777             psDefn->ProjParm[2] = dfStdParallel1;
 1778             psDefn->ProjParmId[2] = ProjStdParallel1GeoKey;
 1779         }
 1780         if( bHaveNOS )
 1781         {
 1782             psDefn->ProjParm[4] = dfNatOriginScale;
 1783             psDefn->ProjParmId[4] = ProjScaleAtNatOriginGeoKey;
 1784         }
 1785         psDefn->ProjParm[5] = dfFalseEasting;
 1786         psDefn->ProjParmId[5] = ProjFalseEastingGeoKey;
 1787         psDefn->ProjParm[6] = dfFalseNorthing;
 1788         psDefn->ProjParmId[6] = ProjFalseNorthingGeoKey;
 1789 
 1790         psDefn->nParms = 7;
 1791         break;
 1792 
 1793 /* -------------------------------------------------------------------- */
 1794       case CT_LambertConfConic_1SP:
 1795       case CT_ObliqueStereographic:
 1796       case CT_TransverseMercator:
 1797       case CT_TransvMercator_SouthOriented:
 1798 /* -------------------------------------------------------------------- */
 1799         if( GTIFKeyGetDOUBLE(psGTIF, ProjNatOriginLongGeoKey,
 1800                        &dfNatOriginLong, 0, 1 ) == 0
 1801             && GTIFKeyGetDOUBLE(psGTIF, ProjFalseOriginLongGeoKey,
 1802                           &dfNatOriginLong, 0, 1 ) == 0
 1803             && GTIFKeyGetDOUBLE(psGTIF, ProjCenterLongGeoKey,
 1804                           &dfNatOriginLong, 0, 1 ) == 0 )
 1805             dfNatOriginLong = 0.0;
 1806 
 1807         if( GTIFKeyGetDOUBLE(psGTIF, ProjNatOriginLatGeoKey,
 1808                        &dfNatOriginLat, 0, 1 ) == 0
 1809             && GTIFKeyGetDOUBLE(psGTIF, ProjFalseOriginLatGeoKey,
 1810                           &dfNatOriginLat, 0, 1 ) == 0
 1811             && GTIFKeyGetDOUBLE(psGTIF, ProjCenterLatGeoKey,
 1812                           &dfNatOriginLat, 0, 1 ) == 0 )
 1813             dfNatOriginLat = 0.0;
 1814 
 1815         if( GTIFKeyGetDOUBLE(psGTIF, ProjScaleAtNatOriginGeoKey,
 1816                        &dfNatOriginScale, 0, 1 ) == 0
 1817             /* See https://github.com/OSGeo/gdal/files/1665718/lasinfo.txt */
 1818             && GTIFKeyGetDOUBLE(psGTIF, ProjScaleAtCenterGeoKey,
 1819                        &dfNatOriginScale, 0, 1 ) == 0 )
 1820             dfNatOriginScale = 1.0;
 1821 
 1822         /* notdef: should transform to decimal degrees at this point */
 1823 
 1824         psDefn->ProjParm[0] = dfNatOriginLat;
 1825         psDefn->ProjParmId[0] = ProjNatOriginLatGeoKey;
 1826         psDefn->ProjParm[1] = dfNatOriginLong;
 1827         psDefn->ProjParmId[1] = ProjNatOriginLongGeoKey;
 1828         psDefn->ProjParm[4] = dfNatOriginScale;
 1829         psDefn->ProjParmId[4] = ProjScaleAtNatOriginGeoKey;
 1830         psDefn->ProjParm[5] = dfFalseEasting;
 1831         psDefn->ProjParmId[5] = ProjFalseEastingGeoKey;
 1832         psDefn->ProjParm[6] = dfFalseNorthing;
 1833         psDefn->ProjParmId[6] = ProjFalseNorthingGeoKey;
 1834 
 1835         psDefn->nParms = 7;
 1836         break;
 1837 
 1838 /* -------------------------------------------------------------------- */
 1839       case CT_ObliqueMercator: /* hotine */
 1840       case CT_HotineObliqueMercatorAzimuthCenter:
 1841 /* -------------------------------------------------------------------- */
 1842         if( GTIFKeyGetDOUBLE(psGTIF, ProjNatOriginLongGeoKey,
 1843                        &dfNatOriginLong, 0, 1 ) == 0
 1844             && GTIFKeyGetDOUBLE(psGTIF, ProjFalseOriginLongGeoKey,
 1845                           &dfNatOriginLong, 0, 1 ) == 0
 1846             && GTIFKeyGetDOUBLE(psGTIF, ProjCenterLongGeoKey,
 1847                           &dfNatOriginLong, 0, 1 ) == 0 )
 1848             dfNatOriginLong = 0.0;
 1849 
 1850         if( GTIFKeyGetDOUBLE(psGTIF, ProjNatOriginLatGeoKey,
 1851                        &dfNatOriginLat, 0, 1 ) == 0
 1852             && GTIFKeyGetDOUBLE(psGTIF, ProjFalseOriginLatGeoKey,
 1853                           &dfNatOriginLat, 0, 1 ) == 0
 1854             && GTIFKeyGetDOUBLE(psGTIF, ProjCenterLatGeoKey,
 1855                           &dfNatOriginLat, 0, 1 ) == 0 )
 1856             dfNatOriginLat = 0.0;
 1857 
 1858         if( GTIFKeyGetDOUBLE(psGTIF, ProjAzimuthAngleGeoKey,
 1859                        &dfAzimuth, 0, 1 ) == 0 )
 1860             dfAzimuth = 0.0;
 1861 
 1862         if( GTIFKeyGetDOUBLE(psGTIF, ProjRectifiedGridAngleGeoKey,
 1863                        &dfRectGridAngle, 0, 1 ) == 0 )
 1864             dfRectGridAngle = 90.0;
 1865 
 1866         if( GTIFKeyGetDOUBLE(psGTIF, ProjScaleAtNatOriginGeoKey,
 1867                        &dfNatOriginScale, 0, 1 ) == 0
 1868             && GTIFKeyGetDOUBLE(psGTIF, ProjScaleAtCenterGeoKey,
 1869                           &dfNatOriginScale, 0, 1 ) == 0 )
 1870             dfNatOriginScale = 1.0;
 1871 
 1872         /* notdef: should transform to decimal degrees at this point */
 1873 
 1874         psDefn->ProjParm[0] = dfNatOriginLat;
 1875         psDefn->ProjParmId[0] = ProjCenterLatGeoKey;
 1876         psDefn->ProjParm[1] = dfNatOriginLong;
 1877         psDefn->ProjParmId[1] = ProjCenterLongGeoKey;
 1878         psDefn->ProjParm[2] = dfAzimuth;
 1879         psDefn->ProjParmId[2] = ProjAzimuthAngleGeoKey;
 1880         psDefn->ProjParm[3] = dfRectGridAngle;
 1881         psDefn->ProjParmId[3] = ProjRectifiedGridAngleGeoKey;
 1882         psDefn->ProjParm[4] = dfNatOriginScale;
 1883         psDefn->ProjParmId[4] = ProjScaleAtCenterGeoKey;
 1884         psDefn->ProjParm[5] = dfFalseEasting;
 1885         psDefn->ProjParmId[5] = ProjFalseEastingGeoKey;
 1886         psDefn->ProjParm[6] = dfFalseNorthing;
 1887         psDefn->ProjParmId[6] = ProjFalseNorthingGeoKey;
 1888 
 1889         psDefn->nParms = 7;
 1890         break;
 1891 
 1892 /* -------------------------------------------------------------------- */
 1893       case CT_ObliqueMercator_Laborde:
 1894 /* -------------------------------------------------------------------- */
 1895         if( GTIFKeyGetDOUBLE(psGTIF, ProjNatOriginLongGeoKey,
 1896                        &dfNatOriginLong, 0, 1 ) == 0
 1897             && GTIFKeyGetDOUBLE(psGTIF, ProjFalseOriginLongGeoKey,
 1898                           &dfNatOriginLong, 0, 1 ) == 0
 1899             && GTIFKeyGetDOUBLE(psGTIF, ProjCenterLongGeoKey,
 1900                           &dfNatOriginLong, 0, 1 ) == 0 )
 1901             dfNatOriginLong = 0.0;
 1902 
 1903         if( GTIFKeyGetDOUBLE(psGTIF, ProjNatOriginLatGeoKey,
 1904                        &dfNatOriginLat, 0, 1 ) == 0
 1905             && GTIFKeyGetDOUBLE(psGTIF, ProjFalseOriginLatGeoKey,
 1906                           &dfNatOriginLat, 0, 1 ) == 0
 1907             && GTIFKeyGetDOUBLE(psGTIF, ProjCenterLatGeoKey,
 1908                           &dfNatOriginLat, 0, 1 ) == 0 )
 1909             dfNatOriginLat = 0.0;
 1910 
 1911         if( GTIFKeyGetDOUBLE(psGTIF, ProjAzimuthAngleGeoKey,
 1912                        &dfAzimuth, 0, 1 ) == 0 )
 1913             dfAzimuth = 0.0;
 1914 
 1915         if( GTIFKeyGetDOUBLE(psGTIF, ProjScaleAtNatOriginGeoKey,
 1916                        &dfNatOriginScale, 0, 1 ) == 0
 1917             && GTIFKeyGetDOUBLE(psGTIF, ProjScaleAtCenterGeoKey,
 1918                           &dfNatOriginScale, 0, 1 ) == 0 )
 1919             dfNatOriginScale = 1.0;
 1920 
 1921         /* notdef: should transform to decimal degrees at this point */
 1922 
 1923         psDefn->ProjParm[0] = dfNatOriginLat;
 1924         psDefn->ProjParmId[0] = ProjCenterLatGeoKey;
 1925         psDefn->ProjParm[1] = dfNatOriginLong;
 1926         psDefn->ProjParmId[1] = ProjCenterLongGeoKey;
 1927         psDefn->ProjParm[2] = dfAzimuth;
 1928         psDefn->ProjParmId[2] = ProjAzimuthAngleGeoKey;
 1929         psDefn->ProjParm[4] = dfNatOriginScale;
 1930         psDefn->ProjParmId[4] = ProjScaleAtCenterGeoKey;
 1931         psDefn->ProjParm[5] = dfFalseEasting;
 1932         psDefn->ProjParmId[5] = ProjFalseEastingGeoKey;
 1933         psDefn->ProjParm[6] = dfFalseNorthing;
 1934         psDefn->ProjParmId[6] = ProjFalseNorthingGeoKey;
 1935 
 1936         psDefn->nParms = 7;
 1937         break;
 1938 
 1939 /* -------------------------------------------------------------------- */
 1940       case CT_CassiniSoldner:
 1941       case CT_Polyconic:
 1942 /* -------------------------------------------------------------------- */
 1943         if( GTIFKeyGetDOUBLE(psGTIF, ProjNatOriginLongGeoKey,
 1944                        &dfNatOriginLong, 0, 1 ) == 0
 1945             && GTIFKeyGetDOUBLE(psGTIF, ProjFalseOriginLongGeoKey,
 1946                           &dfNatOriginLong, 0, 1 ) == 0
 1947             && GTIFKeyGetDOUBLE(psGTIF, ProjCenterLongGeoKey,
 1948                           &dfNatOriginLong, 0, 1 ) == 0 )
 1949             dfNatOriginLong = 0.0;
 1950 
 1951         if( GTIFKeyGetDOUBLE(psGTIF, ProjNatOriginLatGeoKey,
 1952                        &dfNatOriginLat, 0, 1 ) == 0
 1953             && GTIFKeyGetDOUBLE(psGTIF, ProjFalseOriginLatGeoKey,
 1954                           &dfNatOriginLat, 0, 1 ) == 0
 1955             && GTIFKeyGetDOUBLE(psGTIF, ProjCenterLatGeoKey,
 1956                           &dfNatOriginLat, 0, 1 ) == 0 )
 1957             dfNatOriginLat = 0.0;
 1958 
 1959         if( GTIFKeyGetDOUBLE(psGTIF, ProjScaleAtNatOriginGeoKey,
 1960                        &dfNatOriginScale, 0, 1 ) == 0
 1961             && GTIFKeyGetDOUBLE(psGTIF, ProjScaleAtCenterGeoKey,
 1962                           &dfNatOriginScale, 0, 1 ) == 0 )
 1963             dfNatOriginScale = 1.0;
 1964 
 1965         /* notdef: should transform to decimal degrees at this point */
 1966 
 1967         psDefn->ProjParm[0] = dfNatOriginLat;
 1968         psDefn->ProjParmId[0] = ProjNatOriginLatGeoKey;
 1969         psDefn->ProjParm[1] = dfNatOriginLong;
 1970         psDefn->ProjParmId[1] = ProjNatOriginLongGeoKey;
 1971         psDefn->ProjParm[4] = dfNatOriginScale;
 1972         psDefn->ProjParmId[4] = ProjScaleAtNatOriginGeoKey;
 1973         psDefn->ProjParm[5] = dfFalseEasting;
 1974         psDefn->ProjParmId[5] = ProjFalseEastingGeoKey;
 1975         psDefn->ProjParm[6] = dfFalseNorthing;
 1976         psDefn->ProjParmId[6] = ProjFalseNorthingGeoKey;
 1977 
 1978         psDefn->nParms = 7;
 1979         break;
 1980 
 1981 /* -------------------------------------------------------------------- */
 1982       case CT_AzimuthalEquidistant:
 1983       case CT_MillerCylindrical:
 1984       case CT_Gnomonic:
 1985       case CT_LambertAzimEqualArea:
 1986       case CT_Orthographic:
 1987       case CT_NewZealandMapGrid:
 1988 /* -------------------------------------------------------------------- */
 1989         if( GTIFKeyGetDOUBLE(psGTIF, ProjNatOriginLongGeoKey,
 1990                        &dfNatOriginLong, 0, 1 ) == 0
 1991             && GTIFKeyGetDOUBLE(psGTIF, ProjFalseOriginLongGeoKey,
 1992                           &dfNatOriginLong, 0, 1 ) == 0
 1993             && GTIFKeyGetDOUBLE(psGTIF, ProjCenterLongGeoKey,
 1994                           &dfNatOriginLong, 0, 1 ) == 0 )
 1995             dfNatOriginLong = 0.0;
 1996 
 1997         if( GTIFKeyGetDOUBLE(psGTIF, ProjNatOriginLatGeoKey,
 1998                        &dfNatOriginLat, 0, 1 ) == 0
 1999             && GTIFKeyGetDOUBLE(psGTIF, ProjFalseOriginLatGeoKey,
 2000                           &dfNatOriginLat, 0, 1 ) == 0
 2001             && GTIFKeyGetDOUBLE(psGTIF, ProjCenterLatGeoKey,
 2002                           &dfNatOriginLat, 0, 1 ) == 0 )
 2003             dfNatOriginLat = 0.0;
 2004 
 2005         /* notdef: should transform to decimal degrees at this point */
 2006 
 2007         psDefn->ProjParm[0] = dfNatOriginLat;
 2008         psDefn->ProjParmId[0] = ProjCenterLatGeoKey;
 2009         psDefn->ProjParm[1] = dfNatOriginLong;
 2010         psDefn->ProjParmId[1] = ProjCenterLongGeoKey;
 2011         psDefn->ProjParm[5] = dfFalseEasting;
 2012         psDefn->ProjParmId[5] = ProjFalseEastingGeoKey;
 2013         psDefn->ProjParm[6] = dfFalseNorthing;
 2014         psDefn->ProjParmId[6] = ProjFalseNorthingGeoKey;
 2015 
 2016         psDefn->nParms = 7;
 2017         break;
 2018 
 2019 /* -------------------------------------------------------------------- */
 2020       case CT_Equirectangular:
 2021 /* -------------------------------------------------------------------- */
 2022         if( GTIFKeyGetDOUBLE(psGTIF, ProjNatOriginLongGeoKey,
 2023                        &dfNatOriginLong, 0, 1 ) == 0
 2024             && GTIFKeyGetDOUBLE(psGTIF, ProjFalseOriginLongGeoKey,
 2025                           &dfNatOriginLong, 0, 1 ) == 0
 2026             && GTIFKeyGetDOUBLE(psGTIF, ProjCenterLongGeoKey,
 2027                           &dfNatOriginLong, 0, 1 ) == 0 )
 2028             dfNatOriginLong = 0.0;
 2029 
 2030         if( GTIFKeyGetDOUBLE(psGTIF, ProjNatOriginLatGeoKey,
 2031                        &dfNatOriginLat, 0, 1 ) == 0
 2032             && GTIFKeyGetDOUBLE(psGTIF, ProjFalseOriginLatGeoKey,
 2033                           &dfNatOriginLat, 0, 1 ) == 0
 2034             && GTIFKeyGetDOUBLE(psGTIF, ProjCenterLatGeoKey,
 2035                           &dfNatOriginLat, 0, 1 ) == 0 )
 2036             dfNatOriginLat = 0.0;
 2037 
 2038         if( GTIFKeyGetDOUBLE(psGTIF, ProjStdParallel1GeoKey,
 2039                        &dfStdParallel1, 0, 1 ) == 0 )
 2040             dfStdParallel1 = 0.0;
 2041 
 2042         /* notdef: should transform to decimal degrees at this point */
 2043 
 2044         psDefn->ProjParm[0] = dfNatOriginLat;
 2045         psDefn->ProjParmId[0] = ProjCenterLatGeoKey;
 2046         psDefn->ProjParm[1] = dfNatOriginLong;
 2047         psDefn->ProjParmId[1] = ProjCenterLongGeoKey;
 2048         psDefn->ProjParm[2] = dfStdParallel1;
 2049         psDefn->ProjParmId[2] = ProjStdParallel1GeoKey;
 2050         psDefn->ProjParm[5] = dfFalseEasting;
 2051         psDefn->ProjParmId[5] = ProjFalseEastingGeoKey;
 2052         psDefn->ProjParm[6] = dfFalseNorthing;
 2053         psDefn->ProjParmId[6] = ProjFalseNorthingGeoKey;
 2054 
 2055         psDefn->nParms = 7;
 2056         break;
 2057 
 2058 /* -------------------------------------------------------------------- */
 2059       case CT_Robinson:
 2060       case CT_Sinusoidal:
 2061       case CT_VanDerGrinten:
 2062 /* -------------------------------------------------------------------- */
 2063         if( GTIFKeyGetDOUBLE(psGTIF, ProjNatOriginLongGeoKey,
 2064                        &dfNatOriginLong, 0, 1 ) == 0
 2065             && GTIFKeyGetDOUBLE(psGTIF, ProjFalseOriginLongGeoKey,
 2066                           &dfNatOriginLong, 0, 1 ) == 0
 2067             && GTIFKeyGetDOUBLE(psGTIF, ProjCenterLongGeoKey,
 2068                           &dfNatOriginLong, 0, 1 ) == 0 )
 2069             dfNatOriginLong = 0.0;
 2070 
 2071         /* notdef: should transform to decimal degrees at this point */
 2072 
 2073         psDefn->ProjParm[1] = dfNatOriginLong;
 2074         psDefn->ProjParmId[1] = ProjCenterLongGeoKey;
 2075         psDefn->ProjParm[5] = dfFalseEasting;
 2076         psDefn->ProjParmId[5] = ProjFalseEastingGeoKey;
 2077         psDefn->ProjParm[6] = dfFalseNorthing;
 2078         psDefn->ProjParmId[6] = ProjFalseNorthingGeoKey;
 2079 
 2080         psDefn->nParms = 7;
 2081         break;
 2082 
 2083 /* -------------------------------------------------------------------- */
 2084       case CT_PolarStereographic:
 2085 /* -------------------------------------------------------------------- */
 2086         if( GTIFKeyGetDOUBLE(psGTIF, ProjStraightVertPoleLongGeoKey,
 2087                        &dfNatOriginLong, 0, 1 ) == 0
 2088             && GTIFKeyGetDOUBLE(psGTIF, ProjNatOriginLongGeoKey,
 2089                           &dfNatOriginLong, 0, 1 ) == 0
 2090             && GTIFKeyGetDOUBLE(psGTIF, ProjFalseOriginLongGeoKey,
 2091                           &dfNatOriginLong, 0, 1 ) == 0
 2092             && GTIFKeyGetDOUBLE(psGTIF, ProjCenterLongGeoKey,
 2093                           &dfNatOriginLong, 0, 1 ) == 0 )
 2094             dfNatOriginLong = 0.0;
 2095 
 2096         if( GTIFKeyGetDOUBLE(psGTIF, ProjNatOriginLatGeoKey,
 2097                        &dfNatOriginLat, 0, 1 ) == 0
 2098             && GTIFKeyGetDOUBLE(psGTIF, ProjFalseOriginLatGeoKey,
 2099                           &dfNatOriginLat, 0, 1 ) == 0
 2100             && GTIFKeyGetDOUBLE(psGTIF, ProjCenterLatGeoKey,
 2101                           &dfNatOriginLat, 0, 1 ) == 0 )
 2102             dfNatOriginLat = 0.0;
 2103 
 2104         if( GTIFKeyGetDOUBLE(psGTIF, ProjScaleAtNatOriginGeoKey,
 2105                        &dfNatOriginScale, 0, 1 ) == 0
 2106             && GTIFKeyGetDOUBLE(psGTIF, ProjScaleAtCenterGeoKey,
 2107                           &dfNatOriginScale, 0, 1 ) == 0 )
 2108             dfNatOriginScale = 1.0;
 2109 
 2110         /* notdef: should transform to decimal degrees at this point */
 2111 
 2112         psDefn->ProjParm[0] = dfNatOriginLat;
 2113         psDefn->ProjParmId[0] = ProjNatOriginLatGeoKey;;
 2114         psDefn->ProjParm[1] = dfNatOriginLong;
 2115         psDefn->ProjParmId[1] = ProjStraightVertPoleLongGeoKey;
 2116         psDefn->ProjParm[4] = dfNatOriginScale;
 2117         psDefn->ProjParmId[4] = ProjScaleAtNatOriginGeoKey;
 2118         psDefn->ProjParm[5] = dfFalseEasting;
 2119         psDefn->ProjParmId[5] = ProjFalseEastingGeoKey;
 2120         psDefn->ProjParm[6] = dfFalseNorthing;
 2121         psDefn->ProjParmId[6] = ProjFalseNorthingGeoKey;
 2122 
 2123         psDefn->nParms = 7;
 2124         break;
 2125 
 2126 /* -------------------------------------------------------------------- */
 2127       case CT_LambertConfConic_2SP:
 2128 /* -------------------------------------------------------------------- */
 2129         if( GTIFKeyGetDOUBLE(psGTIF, ProjStdParallel1GeoKey,
 2130                        &dfStdParallel1, 0, 1 ) == 0 )
 2131             dfStdParallel1 = 0.0;
 2132 
 2133         if( GTIFKeyGetDOUBLE(psGTIF, ProjStdParallel2GeoKey,
 2134                        &dfStdParallel2, 0, 1 ) == 0 )
 2135             dfStdParallel1 = 0.0;
 2136 
 2137         if( GTIFKeyGetDOUBLE(psGTIF, ProjNatOriginLongGeoKey,
 2138                        &dfNatOriginLong, 0, 1 ) == 0
 2139             && GTIFKeyGetDOUBLE(psGTIF, ProjFalseOriginLongGeoKey,
 2140                           &dfNatOriginLong, 0, 1 ) == 0
 2141             && GTIFKeyGetDOUBLE(psGTIF, ProjCenterLongGeoKey,
 2142                           &dfNatOriginLong, 0, 1 ) == 0 )
 2143             dfNatOriginLong = 0.0;
 2144 
 2145         if( GTIFKeyGetDOUBLE(psGTIF, ProjNatOriginLatGeoKey,
 2146                        &dfNatOriginLat, 0, 1 ) == 0
 2147             && GTIFKeyGetDOUBLE(psGTIF, ProjFalseOriginLatGeoKey,
 2148                           &dfNatOriginLat, 0, 1 ) == 0
 2149             && GTIFKeyGetDOUBLE(psGTIF, ProjCenterLatGeoKey,
 2150                           &dfNatOriginLat, 0, 1 ) == 0 )
 2151             dfNatOriginLat = 0.0;
 2152 
 2153         /* notdef: should transform to decimal degrees at this point */
 2154 
 2155         psDefn->ProjParm[0] = dfNatOriginLat;
 2156         psDefn->ProjParmId[0] = ProjFalseOriginLatGeoKey;
 2157         psDefn->ProjParm[1] = dfNatOriginLong;
 2158         psDefn->ProjParmId[1] = ProjFalseOriginLongGeoKey;
 2159         psDefn->ProjParm[2] = dfStdParallel1;
 2160         psDefn->ProjParmId[2] = ProjStdParallel1GeoKey;
 2161         psDefn->ProjParm[3] = dfStdParallel2;
 2162         psDefn->ProjParmId[3] = ProjStdParallel2GeoKey;
 2163         psDefn->ProjParm[5] = dfFalseEasting;
 2164         psDefn->ProjParmId[5] = ProjFalseEastingGeoKey;
 2165         psDefn->ProjParm[6] = dfFalseNorthing;
 2166         psDefn->ProjParmId[6] = ProjFalseNorthingGeoKey;
 2167 
 2168         psDefn->nParms = 7;
 2169         break;
 2170 
 2171 /* -------------------------------------------------------------------- */
 2172       case CT_AlbersEqualArea:
 2173       case CT_EquidistantConic:
 2174 /* -------------------------------------------------------------------- */
 2175         if( GTIFKeyGetDOUBLE(psGTIF, ProjStdParallel1GeoKey,
 2176                        &dfStdParallel1, 0, 1 ) == 0 )
 2177             dfStdParallel1 = 0.0;
 2178 
 2179         if( GTIFKeyGetDOUBLE(psGTIF, ProjStdParallel2GeoKey,
 2180                        &dfStdParallel2, 0, 1 ) == 0 )
 2181             dfStdParallel2 = 0.0;
 2182 
 2183         if( GTIFKeyGetDOUBLE(psGTIF, ProjNatOriginLongGeoKey,
 2184                        &dfNatOriginLong, 0, 1 ) == 0
 2185             && GTIFKeyGetDOUBLE(psGTIF, ProjFalseOriginLongGeoKey,
 2186                           &dfNatOriginLong, 0, 1 ) == 0
 2187             && GTIFKeyGetDOUBLE(psGTIF, ProjCenterLongGeoKey,
 2188                           &dfNatOriginLong, 0, 1 ) == 0 )
 2189             dfNatOriginLong = 0.0;
 2190 
 2191         if( GTIFKeyGetDOUBLE(psGTIF, ProjNatOriginLatGeoKey,
 2192                        &dfNatOriginLat, 0, 1 ) == 0
 2193             && GTIFKeyGetDOUBLE(psGTIF, ProjFalseOriginLatGeoKey,
 2194                           &dfNatOriginLat, 0, 1 ) == 0
 2195             && GTIFKeyGetDOUBLE(psGTIF, ProjCenterLatGeoKey,
 2196                           &dfNatOriginLat, 0, 1 ) == 0 )
 2197             dfNatOriginLat = 0.0;
 2198 
 2199         /* notdef: should transform to decimal degrees at this point */
 2200 
 2201         psDefn->ProjParm[0] = dfStdParallel1;
 2202         psDefn->ProjParmId[0] = ProjStdParallel1GeoKey;
 2203         psDefn->ProjParm[1] = dfStdParallel2;
 2204         psDefn->ProjParmId[1] = ProjStdParallel2GeoKey;
 2205         psDefn->ProjParm[2] = dfNatOriginLat;
 2206         psDefn->ProjParmId[2] = ProjNatOriginLatGeoKey;
 2207         psDefn->ProjParm[3] = dfNatOriginLong;
 2208         psDefn->ProjParmId[3] = ProjNatOriginLongGeoKey;
 2209         psDefn->ProjParm[5] = dfFalseEasting;
 2210         psDefn->ProjParmId[5] = ProjFalseEastingGeoKey;
 2211         psDefn->ProjParm[6] = dfFalseNorthing;
 2212         psDefn->ProjParmId[6] = ProjFalseNorthingGeoKey;
 2213 
 2214         psDefn->nParms = 7;
 2215         break;
 2216 
 2217 /* -------------------------------------------------------------------- */
 2218       case CT_CylindricalEqualArea:
 2219 /* -------------------------------------------------------------------- */
 2220         if( GTIFKeyGetDOUBLE(psGTIF, ProjStdParallel1GeoKey,
 2221                        &dfStdParallel1, 0, 1 ) == 0 )
 2222             dfStdParallel1 = 0.0;
 2223 
 2224         if( GTIFKeyGetDOUBLE(psGTIF, ProjNatOriginLongGeoKey,
 2225                        &dfNatOriginLong, 0, 1 ) == 0
 2226             && GTIFKeyGetDOUBLE(psGTIF, ProjFalseOriginLongGeoKey,
 2227                           &dfNatOriginLong, 0, 1 ) == 0
 2228             && GTIFKeyGetDOUBLE(psGTIF, ProjCenterLongGeoKey,
 2229                           &dfNatOriginLong, 0, 1 ) == 0 )
 2230             dfNatOriginLong = 0.0;
 2231 
 2232         /* notdef: should transform to decimal degrees at this point */
 2233 
 2234         psDefn->ProjParm[0] = dfStdParallel1;
 2235         psDefn->ProjParmId[0] = ProjStdParallel1GeoKey;
 2236         psDefn->ProjParm[1] = dfNatOriginLong;
 2237         psDefn->ProjParmId[1] = ProjNatOriginLongGeoKey;
 2238         psDefn->ProjParm[5] = dfFalseEasting;
 2239         psDefn->ProjParmId[5] = ProjFalseEastingGeoKey;
 2240         psDefn->ProjParm[6] = dfFalseNorthing;
 2241         psDefn->ProjParmId[6] = ProjFalseNorthingGeoKey;
 2242 
 2243         psDefn->nParms = 7;
 2244         break;
 2245     }
 2246 
 2247 /* -------------------------------------------------------------------- */
 2248 /*      Normalize any linear parameters into meters.  In GeoTIFF        */
 2249 /*      the linear projection parameter tags are normally in the        */
 2250 /*      units of the coordinate system described.                       */
 2251 /* -------------------------------------------------------------------- */
 2252     for( iParm = 0; iParm < psDefn->nParms; iParm++ )
 2253     {
 2254         switch( psDefn->ProjParmId[iParm] )
 2255         {
 2256           case ProjFalseEastingGeoKey:
 2257           case ProjFalseNorthingGeoKey:
 2258           case ProjFalseOriginEastingGeoKey:
 2259           case ProjFalseOriginNorthingGeoKey:
 2260           case ProjCenterEastingGeoKey:
 2261           case ProjCenterNorthingGeoKey:
 2262             if( psDefn->UOMLengthInMeters != 0
 2263                 && psDefn->UOMLengthInMeters != 1.0 )
 2264             {
 2265                 psDefn->ProjParm[iParm] *= psDefn->UOMLengthInMeters;
 2266             }
 2267             break;
 2268 
 2269           default:
 2270             break;
 2271         }
 2272     }
 2273 }
 2274 
 2275 /************************************************************************/
 2276 /*                            GTIFGetDefn()                             */
 2277 /************************************************************************/
 2278 
 2279 /**
 2280 @param psGTIF GeoTIFF information handle as returned by GTIFNew.
 2281 @param psDefn Pointer to an existing GTIFDefn structure allocated by GTIFAllocDefn().
 2282 
 2283 @return TRUE if the function has been successful, otherwise FALSE.
 2284 
 2285 This function reads the coordinate system definition from a GeoTIFF file,
 2286 and <i>normalizes</i> it into a set of component information using
 2287 definitions from the EPSG database as provided by the PROJ library.
 2288 This function is intended to simplify correct support for
 2289 reading files with defined PCS (Projected Coordinate System) codes that
 2290 wouldn't otherwise be directly known by application software by reducing
 2291 it to the underlying projection method, parameters, datum, ellipsoid,
 2292 prime meridian and units.<p>
 2293 
 2294 The application should pass a pointer to an existing uninitialized
 2295 GTIFDefn structure, and GTIFGetDefn() will fill it in.  The function
 2296 currently always returns TRUE but in the future will return FALSE if
 2297 the database is not found.  In any event, all geokeys actually found in the
 2298 file will be copied into the GTIFDefn.  However, if the database isn't
 2299 found, codes implied by other codes will not be set properly.<p>
 2300 
 2301 GTIFGetDefn() will not generally work if the EPSG derived database cannot
 2302 be found.<p>
 2303 
 2304 The normalization methodology operates by fetching tags from the GeoTIFF
 2305 file, and then setting all other tags implied by them in the structure.  The
 2306 implied relationships are worked out by reading definitions from the
 2307 various EPSG derived database tables.<p>
 2308 
 2309 For instance, if a PCS (ProjectedCSTypeGeoKey) is found in the GeoTIFF file
 2310 this code is used to lookup a record in the database.
 2311 For example given the PCS 26746 we can find the name
 2312 (NAD27 / California zone VI), the GCS 4257 (NAD27), and the ProjectionCode
 2313 10406 (California CS27 zone VI).  The GCS, and ProjectionCode can in turn
 2314 be looked up in other tables until all the details of units, ellipsoid,
 2315 prime meridian, datum, projection (LambertConfConic_2SP) and projection
 2316 parameters are established.  A full listgeo dump of a file
 2317 for this result might look like the following, all based on a single PCS
 2318 value:<p>
 2319 
 2320 <pre>
 2321 % listgeo -norm ~/data/geotiff/pci_eg/spaf27.tif
 2322 Geotiff_Information:
 2323    Version: 1
 2324    Key_Revision: 1.0
 2325    Tagged_Information:
 2326       ModelTiepointTag (2,3):
 2327          0                0                0
 2328          1577139.71       634349.176       0
 2329       ModelPixelScaleTag (1,3):
 2330          195.509321       198.32184        0
 2331       End_Of_Tags.
 2332    Keyed_Information:
 2333       GTModelTypeGeoKey (Short,1): ModelTypeProjected
 2334       GTRasterTypeGeoKey (Short,1): RasterPixelIsArea
 2335       ProjectedCSTypeGeoKey (Short,1): PCS_NAD27_California_VI
 2336       End_Of_Keys.
 2337    End_Of_Geotiff.
 2338 
 2339 PCS = 26746 (NAD27 / California zone VI)
 2340 Projection = 10406 (California CS27 zone VI)
 2341 Projection Method: CT_LambertConfConic_2SP
 2342    ProjStdParallel1GeoKey: 33.883333
 2343    ProjStdParallel2GeoKey: 32.766667
 2344    ProjFalseOriginLatGeoKey: 32.166667
 2345    ProjFalseOriginLongGeoKey: -116.233333
 2346    ProjFalseEastingGeoKey: 609601.219202
 2347    ProjFalseNorthingGeoKey: 0.000000
 2348 GCS: 4267/NAD27
 2349 Datum: 6267/North American Datum 1927
 2350 Ellipsoid: 7008/Clarke 1866 (6378206.40,6356583.80)
 2351 Prime Meridian: 8901/Greenwich (0.000000)
 2352 Projection Linear Units: 9003/US survey foot (0.304801m)
 2353 </pre>
 2354 
 2355 Note that GTIFGetDefn() does not inspect or return the tiepoints and scale.
 2356 This must be handled separately as it normally would.  It is intended to
 2357 simplify capture and normalization of the coordinate system definition.
 2358 Note that GTIFGetDefn() also does the following things:
 2359 
 2360 <ol>
 2361 <li> Convert all angular values to decimal degrees.
 2362 <li> Convert all linear values to meters.
 2363 <li> Return the linear units and conversion to meters for the tiepoints and
 2364 scale (though the tiepoints and scale remain in their native units).
 2365 <li> When reading projection parameters a variety of differences between
 2366 different GeoTIFF generators are handled, and a normalized set of parameters
 2367 for each projection are always returned.
 2368 </ol>
 2369 
 2370 Code fields in the GTIFDefn are filled with KvUserDefined if there is not
 2371 value to assign.  The parameter lists for each of the underlying projection
 2372 transform methods can be found at the
 2373 <a href="http://www.remotesensing.org/geotiff/proj_list">Projections</a>
 2374 page.  Note that nParms will be set based on the maximum parameter used.
 2375 Some of the parameters may not be used in which case the
 2376 GTIFDefn::ProjParmId[] will
 2377 be zero.  This is done to retain correspondence to the EPSG parameter
 2378 numbering scheme.<p>
 2379 
 2380 The
 2381 <a href="http://www.remotesensing.org/cgi-bin/cvsweb.cgi/~checkout~/osrs/geotiff/libgeotiff/geotiff_proj4.c">geotiff_proj4.c</a> module distributed with libgeotiff can
 2382 be used as an example of code that converts a GTIFDefn into another projection
 2383 system.<p>
 2384 
 2385 @see GTIFKeySet()
 2386 
 2387 */
 2388 
 2389 int GTIFGetDefn( GTIF * psGTIF, GTIFDefn * psDefn )
 2390 
 2391 {
 2392     int     i;
 2393     short   nGeogUOMLinear;
 2394     double  dfInvFlattening;
 2395 
 2396     if( !GTIFGetPROJContext(psGTIF, TRUE, NULL) )
 2397     {
 2398         return FALSE;
 2399     }
 2400 
 2401 /* -------------------------------------------------------------------- */
 2402 /*      Initially we default all the information we can.                */
 2403 /* -------------------------------------------------------------------- */
 2404     psDefn->DefnSet = 1;
 2405     psDefn->Model = KvUserDefined;
 2406     psDefn->PCS = KvUserDefined;
 2407     psDefn->GCS = KvUserDefined;
 2408     psDefn->UOMLength = KvUserDefined;
 2409     psDefn->UOMLengthInMeters = 1.0;
 2410     psDefn->UOMAngle = KvUserDefined;
 2411     psDefn->UOMAngleInDegrees = 1.0;
 2412     psDefn->Datum = KvUserDefined;
 2413     psDefn->Ellipsoid = KvUserDefined;
 2414     psDefn->SemiMajor = 0.0;
 2415     psDefn->SemiMinor = 0.0;
 2416     psDefn->PM = KvUserDefined;
 2417     psDefn->PMLongToGreenwich = 0.0;
 2418 #if !defined(GEO_NORMALIZE_DISABLE_TOWGS84)
 2419     psDefn->TOWGS84Count = 0;
 2420     memset( psDefn->TOWGS84, 0, sizeof(psDefn->TOWGS84) );
 2421 #endif
 2422 
 2423     psDefn->ProjCode = KvUserDefined;
 2424     psDefn->Projection = KvUserDefined;
 2425     psDefn->CTProjection = KvUserDefined;
 2426 
 2427     psDefn->nParms = 0;
 2428     for( i = 0; i < MAX_GTIF_PROJPARMS; i++ )
 2429     {
 2430         psDefn->ProjParm[i] = 0.0;
 2431         psDefn->ProjParmId[i] = 0;
 2432     }
 2433 
 2434     psDefn->MapSys = KvUserDefined;
 2435     psDefn->Zone = 0;
 2436 
 2437 /* -------------------------------------------------------------------- */
 2438 /*      Do we have any geokeys?                                         */
 2439 /* -------------------------------------------------------------------- */
 2440     {
 2441         int     nKeyCount = 0;
 2442         int     anVersion[3];
 2443         GTIFDirectoryInfo( psGTIF, anVersion, &nKeyCount );
 2444 
 2445         if( nKeyCount == 0 )
 2446         {
 2447             psDefn->DefnSet = 0;
 2448             return FALSE;
 2449         }
 2450     }
 2451 
 2452 /* -------------------------------------------------------------------- */
 2453 /*  Try to get the overall model type.              */
 2454 /* -------------------------------------------------------------------- */
 2455     GTIFKeyGetSSHORT(psGTIF,GTModelTypeGeoKey,&(psDefn->Model));
 2456 
 2457 /* -------------------------------------------------------------------- */
 2458 /*  Extract the Geog units.                     */
 2459 /* -------------------------------------------------------------------- */
 2460     nGeogUOMLinear = 9001; /* Linear_Meter */
 2461     if( GTIFKeyGetSSHORT(psGTIF, GeogLinearUnitsGeoKey, &nGeogUOMLinear) == 1 )
 2462     {
 2463         psDefn->UOMLength = nGeogUOMLinear;
 2464     }
 2465 
 2466 /* -------------------------------------------------------------------- */
 2467 /*      Try to get a PCS.                                               */
 2468 /* -------------------------------------------------------------------- */
 2469     if( GTIFKeyGetSSHORT(psGTIF,ProjectedCSTypeGeoKey, &(psDefn->PCS)) == 1
 2470         && psDefn->PCS != KvUserDefined )
 2471     {
 2472         /*
 2473          * Translate this into useful information.
 2474          */
 2475         GTIFGetPCSInfoEx( psGTIF->pj_context,
 2476                           psDefn->PCS, NULL, &(psDefn->ProjCode),
 2477                           &(psDefn->UOMLength), &(psDefn->GCS) );
 2478     }
 2479 
 2480 /* -------------------------------------------------------------------- */
 2481 /*       If we have the PCS code, but didn't find it in the database    */
 2482 /*      (likely because we can't find them) we will try some ``jiffy    */
 2483 /*      rules'' for UTM and state plane.                                */
 2484 /* -------------------------------------------------------------------- */
 2485     if( psDefn->PCS != KvUserDefined && psDefn->ProjCode == KvUserDefined )
 2486     {
 2487         int nMapSys, nZone;
 2488         int nGCS = psDefn->GCS;
 2489 
 2490         nMapSys = GTIFPCSToMapSys( psDefn->PCS, &nGCS, &nZone );
 2491         if( nMapSys != KvUserDefined )
 2492         {
 2493             psDefn->ProjCode = (short) GTIFMapSysToProj( nMapSys, nZone );
 2494             psDefn->GCS = (short) nGCS;
 2495         }
 2496     }
 2497 
 2498 /* -------------------------------------------------------------------- */
 2499 /*      If the Proj_ code is specified directly, use that.              */
 2500 /* -------------------------------------------------------------------- */
 2501     if( psDefn->ProjCode == KvUserDefined )
 2502         GTIFKeyGetSSHORT(psGTIF, ProjectionGeoKey, &(psDefn->ProjCode));
 2503 
 2504     if( psDefn->ProjCode != KvUserDefined )
 2505     {
 2506         /*
 2507          * We have an underlying projection transformation value.  Look
 2508          * this up.  For a PCS of ``WGS 84 / UTM 11'' the transformation
 2509          * would be Transverse Mercator, with a particular set of options.
 2510          * The nProjTRFCode itself would correspond to the name
 2511          * ``UTM zone 11N'', and doesn't include datum info.
 2512          */
 2513         GTIFGetProjTRFInfoEx( psGTIF->pj_context,
 2514                               psDefn->ProjCode, NULL, &(psDefn->Projection),
 2515                               psDefn->ProjParm );
 2516 
 2517         /*
 2518          * Set the GeoTIFF identity of the parameters.
 2519          */
 2520         psDefn->CTProjection = (short)
 2521             EPSGProjMethodToCTProjMethod( psDefn->Projection, FALSE );
 2522 
 2523         SetGTParmIds( EPSGProjMethodToCTProjMethod(psDefn->Projection, TRUE),
 2524                       psDefn->Projection,
 2525                       psDefn->ProjParmId, NULL);
 2526         psDefn->nParms = 7;
 2527     }
 2528 
 2529 /* -------------------------------------------------------------------- */
 2530 /*      Try to get a GCS.  If found, it will override any implied by    */
 2531 /*      the PCS.                                                        */
 2532 /* -------------------------------------------------------------------- */
 2533     GTIFKeyGetSSHORT(psGTIF, GeographicTypeGeoKey, &(psDefn->GCS));
 2534     if( psDefn->GCS < 1 || psDefn->GCS >= KvUserDefined )
 2535         psDefn->GCS = KvUserDefined;
 2536 
 2537 /* -------------------------------------------------------------------- */
 2538 /*      Derive the datum, and prime meridian from the GCS.              */
 2539 /* -------------------------------------------------------------------- */
 2540     if( psDefn->GCS != KvUserDefined )
 2541     {
 2542         GTIFGetGCSInfoEx( psGTIF->pj_context,
 2543                           psDefn->GCS, NULL, &(psDefn->Datum), &(psDefn->PM),
 2544                           &(psDefn->UOMAngle) );
 2545     }
 2546 
 2547 /* -------------------------------------------------------------------- */
 2548 /*      Handle the GCS angular units.  GeogAngularUnitsGeoKey           */
 2549 /*      overrides the GCS or PCS setting.                               */
 2550 /* -------------------------------------------------------------------- */
 2551     GTIFKeyGetSSHORT(psGTIF, GeogAngularUnitsGeoKey, &(psDefn->UOMAngle));
 2552     if( psDefn->UOMAngle != KvUserDefined )
 2553     {
 2554         GTIFGetUOMAngleInfoEx( psGTIF->pj_context,
 2555                                psDefn->UOMAngle, NULL,
 2556                                &(psDefn->UOMAngleInDegrees) );
 2557     }
 2558 
 2559 /* -------------------------------------------------------------------- */
 2560 /*      Check for a datum setting, and then use the datum to derive     */
 2561 /*      an ellipsoid.                                                   */
 2562 /* -------------------------------------------------------------------- */
 2563     GTIFKeyGetSSHORT(psGTIF, GeogGeodeticDatumGeoKey, &(psDefn->Datum));
 2564 
 2565     if( psDefn->Datum != KvUserDefined )
 2566     {
 2567         GTIFGetDatumInfoEx( psGTIF->pj_context,
 2568                             psDefn->Datum, NULL, &(psDefn->Ellipsoid) );
 2569     }
 2570 
 2571 /* -------------------------------------------------------------------- */
 2572 /*      Check for an explicit ellipsoid.  Use the ellipsoid to          */
 2573 /*      derive the ellipsoid characteristics, if possible.              */
 2574 /* -------------------------------------------------------------------- */
 2575     GTIFKeyGetSSHORT(psGTIF, GeogEllipsoidGeoKey, &(psDefn->Ellipsoid));
 2576 
 2577     if( psDefn->Ellipsoid != KvUserDefined )
 2578     {
 2579         GTIFGetEllipsoidInfoEx( psGTIF->pj_context,
 2580                                 psDefn->Ellipsoid, NULL,
 2581                                 &(psDefn->SemiMajor), &(psDefn->SemiMinor) );
 2582     }
 2583 
 2584 /* -------------------------------------------------------------------- */
 2585 /*      Check for overridden ellipsoid parameters.  It would be nice    */
 2586 /*      to warn if they conflict with provided information, but for     */
 2587 /*      now we just override.                                           */
 2588 /* -------------------------------------------------------------------- */
 2589     CPL_IGNORE_RET_VAL_INT(GTIFKeyGetDOUBLE(psGTIF, GeogSemiMajorAxisGeoKey, &(psDefn->SemiMajor), 0, 1 ));
 2590     CPL_IGNORE_RET_VAL_INT(GTIFKeyGetDOUBLE(psGTIF, GeogSemiMinorAxisGeoKey, &(psDefn->SemiMinor), 0, 1 ));
 2591 
 2592     if( GTIFKeyGetDOUBLE(psGTIF, GeogInvFlatteningGeoKey, &dfInvFlattening,
 2593                    0, 1 ) == 1 )
 2594     {
 2595         if( dfInvFlattening != 0.0 )
 2596             psDefn->SemiMinor =
 2597                 psDefn->SemiMajor * (1 - 1.0/dfInvFlattening);
 2598         else
 2599             psDefn->SemiMinor = psDefn->SemiMajor;
 2600     }
 2601 
 2602 /* -------------------------------------------------------------------- */
 2603 /*      Get the prime meridian info.                                    */
 2604 /* -------------------------------------------------------------------- */
 2605     GTIFKeyGetSSHORT(psGTIF, GeogPrimeMeridianGeoKey, &(psDefn->PM));
 2606 
 2607     if( psDefn->PM != KvUserDefined )
 2608     {
 2609         GTIFGetPMInfoEx( psGTIF->pj_context,
 2610                          psDefn->PM, NULL, &(psDefn->PMLongToGreenwich) );
 2611     }
 2612     else
 2613     {
 2614         CPL_IGNORE_RET_VAL_INT(GTIFKeyGetDOUBLE(psGTIF, GeogPrimeMeridianLongGeoKey,
 2615                    &(psDefn->PMLongToGreenwich), 0, 1 ));
 2616 
 2617         psDefn->PMLongToGreenwich =
 2618             GTIFAngleToDD( psDefn->PMLongToGreenwich,
 2619                            psDefn->UOMAngle );
 2620     }
 2621 
 2622 /* -------------------------------------------------------------------- */
 2623 /*      Get the TOWGS84 parameters.                                     */
 2624 /* -------------------------------------------------------------------- */
 2625 #if !defined(GEO_NORMALIZE_DISABLE_TOWGS84)
 2626     psDefn->TOWGS84Count =
 2627         (short)GTIFKeyGetDOUBLE(psGTIF, GeogTOWGS84GeoKey, psDefn->TOWGS84, 0, 7 );
 2628 #endif
 2629 
 2630 /* -------------------------------------------------------------------- */
 2631 /*      Have the projection units of measure been overridden?  We       */
 2632 /*      should likely be doing something about angular units too,       */
 2633 /*      but these are very rarely not decimal degrees for actual        */
 2634 /*      file coordinates.                                               */
 2635 /* -------------------------------------------------------------------- */
 2636     GTIFKeyGetSSHORT(psGTIF,ProjLinearUnitsGeoKey,&(psDefn->UOMLength));
 2637 
 2638     if( psDefn->UOMLength != KvUserDefined )
 2639     {
 2640         GTIFGetUOMLengthInfoEx( psGTIF->pj_context,
 2641                                 psDefn->UOMLength, NULL,
 2642                                 &(psDefn->UOMLengthInMeters) );
 2643     }
 2644     else
 2645     {
 2646         CPL_IGNORE_RET_VAL_INT(GTIFKeyGetDOUBLE(psGTIF,ProjLinearUnitSizeGeoKey,&(psDefn->UOMLengthInMeters),0,1));
 2647     }
 2648 
 2649 /* -------------------------------------------------------------------- */
 2650 /*      Handle a variety of user defined transform types.               */
 2651 /* -------------------------------------------------------------------- */
 2652     if( GTIFKeyGetSSHORT(psGTIF,ProjCoordTransGeoKey,
 2653                    &(psDefn->CTProjection)) == 1)
 2654     {
 2655         GTIFFetchProjParms( psGTIF, psDefn );
 2656     }
 2657 
 2658 /* -------------------------------------------------------------------- */
 2659 /*      Try to set the zoned map system information.                    */
 2660 /* -------------------------------------------------------------------- */
 2661     psDefn->MapSys = GTIFProjToMapSys( psDefn->ProjCode, &(psDefn->Zone) );
 2662 
 2663 /* -------------------------------------------------------------------- */
 2664 /*      If this is UTM, and we were unable to extract the projection    */
 2665 /*      parameters from the database just set them directly now,        */
 2666 /*      since it's pretty easy, and a common case.                      */
 2667 /* -------------------------------------------------------------------- */
 2668     if( (psDefn->MapSys == MapSys_UTM_North
 2669          || psDefn->MapSys == MapSys_UTM_South)
 2670         && psDefn->CTProjection == KvUserDefined )
 2671     {
 2672         psDefn->CTProjection = CT_TransverseMercator;
 2673         psDefn->nParms = 7;
 2674         psDefn->ProjParmId[0] = ProjNatOriginLatGeoKey;
 2675         psDefn->ProjParm[0] = 0.0;
 2676 
 2677         psDefn->ProjParmId[1] = ProjNatOriginLongGeoKey;
 2678         psDefn->ProjParm[1] = psDefn->Zone*6 - 183.0;
 2679 
 2680         psDefn->ProjParmId[4] = ProjScaleAtNatOriginGeoKey;
 2681         psDefn->ProjParm[4] = 0.9996;
 2682 
 2683         psDefn->ProjParmId[5] = ProjFalseEastingGeoKey;
 2684         psDefn->ProjParm[5] = 500000.0;
 2685 
 2686         psDefn->ProjParmId[6] = ProjFalseNorthingGeoKey;
 2687 
 2688         if( psDefn->MapSys == MapSys_UTM_North )
 2689             psDefn->ProjParm[6] = 0.0;
 2690         else
 2691             psDefn->ProjParm[6] = 10000000.0;
 2692     }
 2693 
 2694     return TRUE;
 2695 }
 2696 
 2697 /************************************************************************/
 2698 /*                            GTIFDecToDMS()                            */
 2699 /*                                                                      */
 2700 /*      Convenient function to translate decimal degrees to DMS         */
 2701 /*      format for reporting to a user.                                 */
 2702 /************************************************************************/
 2703 
 2704 const char *GTIFDecToDMS( double dfAngle, const char * pszAxis,
 2705                           int nPrecision )
 2706 
 2707 {
 2708     int     nDegrees, nMinutes;
 2709     double  dfSeconds;
 2710     char    szFormat[30];
 2711     static char szBuffer[50];
 2712     const char  *pszHemisphere = NULL;
 2713     double  dfRound;
 2714     int     i;
 2715 
 2716     if( !(dfAngle >= -360 && dfAngle <= 360) )
 2717         return "";
 2718 
 2719     dfRound = 0.5/60;
 2720     for( i = 0; i < nPrecision; i++ )
 2721         dfRound = dfRound * 0.1;
 2722 
 2723     nDegrees = (int) ABS(dfAngle);
 2724     nMinutes = (int) ((ABS(dfAngle) - nDegrees) * 60 + dfRound);
 2725     if( nMinutes == 60 )
 2726     {
 2727         nDegrees ++;
 2728         nMinutes = 0;
 2729     }
 2730     dfSeconds = ABS((ABS(dfAngle) * 3600 - nDegrees*3600 - nMinutes*60));
 2731 
 2732     if( EQUAL(pszAxis,"Long") && dfAngle < 0.0 )
 2733         pszHemisphere = "W";
 2734     else if( EQUAL(pszAxis,"Long") )
 2735         pszHemisphere = "E";
 2736     else if( dfAngle < 0.0 )
 2737         pszHemisphere = "S";
 2738     else
 2739         pszHemisphere = "N";
 2740 
 2741     sprintf( szFormat, "%%3dd%%2d\'%%%d.%df\"%s",
 2742              nPrecision+3, nPrecision, pszHemisphere );
 2743     sprintf( szBuffer, szFormat, nDegrees, nMinutes, dfSeconds );
 2744 
 2745     return szBuffer;
 2746 }
 2747 
 2748 /************************************************************************/
 2749 /*                           GTIFPrintDefn()                            */
 2750 /*                                                                      */
 2751 /*      Report the contents of a GTIFDefn structure ... mostly for      */
 2752 /*      debugging.                                                      */
 2753 /************************************************************************/
 2754 
 2755 void GTIFPrintDefnEx( GTIF *psGTIF, GTIFDefn * psDefn, FILE * fp )
 2756 
 2757 {
 2758     GTIFGetPROJContext(psGTIF, TRUE, NULL);
 2759 
 2760 /* -------------------------------------------------------------------- */
 2761 /*      Do we have anything to report?                                  */
 2762 /* -------------------------------------------------------------------- */
 2763     if( !psDefn->DefnSet )
 2764     {
 2765         fprintf( fp, "No GeoKeys found.\n" );
 2766         return;
 2767     }
 2768 
 2769 /* -------------------------------------------------------------------- */
 2770 /*      Get the PCS name if possible.                                   */
 2771 /* -------------------------------------------------------------------- */
 2772     if( psDefn->PCS != KvUserDefined )
 2773     {
 2774         char    *pszPCSName = NULL;
 2775 
 2776         if( psGTIF->pj_context )
 2777         {
 2778             GTIFGetPCSInfoEx( psGTIF->pj_context,
 2779                               psDefn->PCS, &pszPCSName, NULL, NULL, NULL );
 2780         }
 2781         if( pszPCSName == NULL )
 2782             pszPCSName = CPLStrdup("name unknown");
 2783 
 2784         fprintf( fp, "PCS = %d (%s)\n", psDefn->PCS, pszPCSName );
 2785         CPLFree( pszPCSName );
 2786     }
 2787 
 2788 /* -------------------------------------------------------------------- */
 2789 /*  Dump the projection code if possible.               */
 2790 /* -------------------------------------------------------------------- */
 2791     if( psDefn->ProjCode != KvUserDefined )
 2792     {
 2793         char    *pszTRFName = NULL;
 2794 
 2795         if( psGTIF->pj_context )
 2796         {
 2797             GTIFGetProjTRFInfoEx( psGTIF->pj_context,
 2798                                   psDefn->ProjCode, &pszTRFName, NULL, NULL );
 2799         }
 2800         if( pszTRFName == NULL )
 2801             pszTRFName = CPLStrdup("");
 2802 
 2803         fprintf( fp, "Projection = %d (%s)\n",
 2804                  psDefn->ProjCode, pszTRFName );
 2805 
 2806         CPLFree( pszTRFName );
 2807     }
 2808 
 2809 /* -------------------------------------------------------------------- */
 2810 /*      Try to dump the projection method name, and parameters if possible.*/
 2811 /* -------------------------------------------------------------------- */
 2812     if( psDefn->CTProjection != KvUserDefined )
 2813     {
 2814         const char *pszProjectionMethodName =
 2815             GTIFValueNameEx(psGTIF,
 2816                             ProjCoordTransGeoKey,
 2817                             psDefn->CTProjection);
 2818         int     i;
 2819 
 2820         if( pszProjectionMethodName == NULL )
 2821             pszProjectionMethodName = "(unknown)";
 2822 
 2823         fprintf( fp, "Projection Method: %s\n", pszProjectionMethodName );
 2824 
 2825         for( i = 0; i < psDefn->nParms; i++ )
 2826         {
 2827             char* pszName;
 2828             if( psDefn->ProjParmId[i] == 0 )
 2829                 continue;
 2830 
 2831             pszName = GTIFKeyName((geokey_t) psDefn->ProjParmId[i]);
 2832             if( pszName == NULL )
 2833                 pszName = "(unknown)";
 2834 
 2835             if( i < 4 )
 2836             {
 2837                 char    *pszAxisName;
 2838 
 2839                 if( strstr(pszName,"Long") != NULL )
 2840                     pszAxisName = "Long";
 2841                 else if( strstr(pszName,"Lat") != NULL )
 2842                     pszAxisName = "Lat";
 2843                 else
 2844                     pszAxisName = "?";
 2845 
 2846                 fprintf( fp, "   %s: %f (%s)\n",
 2847                          pszName, psDefn->ProjParm[i],
 2848                          GTIFDecToDMS( psDefn->ProjParm[i], pszAxisName, 2 ) );
 2849             }
 2850             else if( i == 4 )
 2851                 fprintf( fp, "   %s: %f\n", pszName, psDefn->ProjParm[i] );
 2852             else
 2853                 fprintf( fp, "   %s: %f m\n", pszName, psDefn->ProjParm[i] );
 2854         }
 2855     }
 2856 
 2857 /* -------------------------------------------------------------------- */
 2858 /*      Report the GCS name, and number.                                */
 2859 /* -------------------------------------------------------------------- */
 2860     if( psDefn->GCS != KvUserDefined )
 2861     {
 2862         char    *pszName = NULL;
 2863 
 2864         if( psGTIF->pj_context )
 2865         {
 2866             GTIFGetGCSInfoEx( psGTIF->pj_context,
 2867                               psDefn->GCS, &pszName, NULL, NULL, NULL );
 2868         }
 2869         if( pszName == NULL )
 2870             pszName = CPLStrdup("(unknown)");
 2871 
 2872         fprintf( fp, "GCS: %d/%s\n", psDefn->GCS, pszName );
 2873         CPLFree( pszName );
 2874     }
 2875 
 2876 /* -------------------------------------------------------------------- */
 2877 /*      Report the datum name.                                          */
 2878 /* -------------------------------------------------------------------- */
 2879     if( psDefn->Datum != KvUserDefined )
 2880     {
 2881         char    *pszName = NULL;
 2882 
 2883         if( psGTIF->pj_context )
 2884         {
 2885             GTIFGetDatumInfoEx( psGTIF->pj_context,
 2886                                 psDefn->Datum, &pszName, NULL );
 2887         }
 2888         if( pszName == NULL )
 2889             pszName = CPLStrdup("(unknown)");
 2890 
 2891         fprintf( fp, "Datum: %d/%s\n", psDefn->Datum, pszName );
 2892         CPLFree( pszName );
 2893     }
 2894 
 2895 /* -------------------------------------------------------------------- */
 2896 /*      Report the ellipsoid.                                           */
 2897 /* -------------------------------------------------------------------- */
 2898     if( psDefn->Ellipsoid != KvUserDefined )
 2899     {
 2900         char    *pszName = NULL;
 2901 
 2902         if( psGTIF->pj_context )
 2903         {
 2904             GTIFGetEllipsoidInfoEx( psGTIF->pj_context,
 2905                                     psDefn->Ellipsoid, &pszName, NULL, NULL );
 2906         }
 2907         if( pszName == NULL )
 2908             pszName = CPLStrdup("(unknown)");
 2909 
 2910         fprintf( fp, "Ellipsoid: %d/%s (%.2f,%.2f)\n",
 2911                  psDefn->Ellipsoid, pszName,
 2912                  psDefn->SemiMajor, psDefn->SemiMinor );
 2913         CPLFree( pszName );
 2914     }
 2915 
 2916 /* -------------------------------------------------------------------- */
 2917 /*      Report the prime meridian.                                      */
 2918 /* -------------------------------------------------------------------- */
 2919     if( psDefn->PM != KvUserDefined )
 2920     {
 2921         char    *pszName = NULL;
 2922 
 2923         if( psGTIF->pj_context )
 2924         {
 2925             GTIFGetPMInfoEx( psGTIF->pj_context,
 2926                              psDefn->PM, &pszName, NULL );
 2927         }
 2928 
 2929         if( pszName == NULL )
 2930             pszName = CPLStrdup("(unknown)");
 2931 
 2932         fprintf( fp, "Prime Meridian: %d/%s (%f/%s)\n",
 2933                  psDefn->PM, pszName,
 2934                  psDefn->PMLongToGreenwich,
 2935                  GTIFDecToDMS( psDefn->PMLongToGreenwich, "Long", 2 ) );
 2936         CPLFree( pszName );
 2937     }
 2938 
 2939 /* -------------------------------------------------------------------- */
 2940 /*      Report TOWGS84 parameters.                                      */
 2941 /* -------------------------------------------------------------------- */
 2942 #if !defined(GEO_NORMALIZE_DISABLE_TOWGS84)
 2943     if( psDefn->TOWGS84Count > 0 )
 2944     {
 2945         int i;
 2946 
 2947         fprintf( fp, "TOWGS84: " );
 2948 
 2949         for( i = 0; i < psDefn->TOWGS84Count; i++ )
 2950         {
 2951             if( i > 0 )
 2952                 fprintf( fp, "," );
 2953             fprintf( fp, "%g", psDefn->TOWGS84[i] );
 2954         }
 2955 
 2956         fprintf( fp, "\n" );
 2957     }
 2958 #endif
 2959 
 2960 /* -------------------------------------------------------------------- */
 2961 /*      Report the projection units of measure (currently just          */
 2962 /*      linear).                                                        */
 2963 /* -------------------------------------------------------------------- */
 2964     if( psDefn->UOMLength != KvUserDefined )
 2965     {
 2966         char    *pszName = NULL;
 2967 
 2968         if( psGTIF->pj_context )
 2969         {
 2970             GTIFGetUOMLengthInfoEx(
 2971                 psGTIF->pj_context, psDefn->UOMLength, &pszName, NULL );
 2972         }
 2973         if( pszName == NULL )
 2974             pszName = CPLStrdup( "(unknown)" );
 2975 
 2976         fprintf( fp, "Projection Linear Units: %d/%s (%fm)\n",
 2977                  psDefn->UOMLength, pszName, psDefn->UOMLengthInMeters );
 2978         CPLFree( pszName );
 2979     }
 2980     else
 2981     {
 2982         fprintf( fp, "Projection Linear Units: User-Defined (%fm)\n",
 2983                  psDefn->UOMLengthInMeters );
 2984     }
 2985 }
 2986 
 2987 void GTIFPrintDefn( GTIFDefn * psDefn, FILE * fp )
 2988 {
 2989     GTIF *psGTIF = GTIFNew(NULL);
 2990     GTIFPrintDefnEx(psGTIF, psDefn, fp);
 2991     GTIFFree(psGTIF);
 2992 }
 2993 
 2994 /************************************************************************/
 2995 /*                           GTIFFreeMemory()                           */
 2996 /*                                                                      */
 2997 /*      Externally visible function to free memory allocated within     */
 2998 /*      geo_normalize.c.                                                */
 2999 /************************************************************************/
 3000 
 3001 void GTIFFreeMemory( char * pMemory )
 3002 
 3003 {
 3004     if( pMemory != NULL )
 3005         VSIFree( pMemory );
 3006 }
 3007 
 3008 /************************************************************************/
 3009 /*                           GTIFAllocDefn()                            */
 3010 /*                                                                      */
 3011 /*      This allocates a GTIF structure in such a way that the          */
 3012 /*      calling application doesn't need to know the size and           */
 3013 /*      initializes it appropriately.                                   */
 3014 /************************************************************************/
 3015 
 3016 GTIFDefn *GTIFAllocDefn()
 3017 {
 3018     return (GTIFDefn *) CPLCalloc(sizeof(GTIFDefn),1);
 3019 }
 3020 
 3021 /************************************************************************/
 3022 /*                            GTIFFreeDefn()                            */
 3023 /*                                                                      */
 3024 /*      Free a GTIF structure allocated by GTIFAllocDefn().             */
 3025 /************************************************************************/
 3026 
 3027 void GTIFFreeDefn( GTIFDefn *defn )
 3028 {
 3029     VSIFree( defn );
 3030 }
 3031 
 3032 /************************************************************************/
 3033 /*                       GTIFAttachPROJContext()                        */
 3034 /*                                                                      */
 3035 /*      Attach an existing PROJ context to the GTIF handle, but         */
 3036 /*      ownership of the context remains to the caller.                 */
 3037 /************************************************************************/
 3038 
 3039 void GTIFAttachPROJContext( GTIF *psGTIF, void* pjContext )
 3040 {
 3041     if( psGTIF->own_pj_context )
 3042     {
 3043         proj_context_destroy(psGTIF->pj_context);
 3044     }
 3045     psGTIF->own_pj_context = FALSE;
 3046     psGTIF->pj_context = (PJ_CONTEXT*) pjContext;
 3047 }
 3048 
 3049 /************************************************************************/
 3050 /*                         GTIFGetPROJContext()                         */
 3051 /*                                                                      */
 3052 /*      Return the PROJ context attached to the GTIF handle.            */
 3053 /*      If it has not yet been instantiated and instantiateIfNeeded=TRUE*/
 3054 /*      then, it will be instantiated (and owned by GTIF handle).       */
 3055 /************************************************************************/
 3056 
 3057 void *GTIFGetPROJContext( GTIF *psGTIF, int instantiateIfNeeded,
 3058                           int* out_gtif_own_pj_context )
 3059 {
 3060     if( psGTIF->pj_context || !instantiateIfNeeded )
 3061     {
 3062         if( out_gtif_own_pj_context )
 3063         {
 3064             *out_gtif_own_pj_context = psGTIF->own_pj_context;
 3065         }
 3066         return psGTIF->pj_context;
 3067     }
 3068     psGTIF->pj_context = proj_context_create();
 3069     psGTIF->own_pj_context = psGTIF->pj_context != NULL;
 3070     if( out_gtif_own_pj_context )
 3071     {
 3072         *out_gtif_own_pj_context = psGTIF->own_pj_context;
 3073     }
 3074     return psGTIF->pj_context;
 3075 }
 3076 
 3077 
 3078 void GTIFDeaccessCSV( void )
 3079 {
 3080     /* No operation */
 3081 }
 3082 
 3083 #ifndef GDAL_COMPILATION
 3084 void SetCSVFilenameHook( const char *(*CSVFileOverride)(const char *) )
 3085 {
 3086     (void)CSVFileOverride;
 3087     /* No operation */
 3088 }
 3089 #endif