[gdal-dev] [PATCH] Support Mercator_2SP in GeoTIFF

[Trent Piepho]

Mercator_[12]SP are both the same projection in GeoTIFF.  Lacking a
definition in the official specification, the intention appears to be that
for Mercator_2SP the latitude of true scale (lat_ts) should be specified in
ProjStdParallel1GeoKey and for Mercator_2SP the scale at origin (k) should
be specified in ProjScaleAtNatOriginGeoKey.

Current behavior when creating a GTiff with 2SP is to drop lat_ts and create
a default k value of 1.0, which will product an incorrect projection.  When
reading, a default k of 1.0 is used and lat_ts is ignored, also changing the
projection.

This patches changes the behavior to supply ProjScaleAtNatOriginGeoKey for
Mercator_1SP or ProjStdParallel1GeoKey for Mercator_2SP when writing.

For reading either or both values are supplied when present by libgeotiff's
normalization code.  If neither is present, a default k of 1.0 is created.

GDAL's GTiff driver is changed to create a 2SP projection is lat_ts is
present, 1SP otherwise.  It emits a warning of both last_ts and scale are
present (and ignores scale).
---
 gdal/frmts/gtiff/gt_wkt_srs.cpp             |   28 ++++++++++---
 gdal/frmts/gtiff/libgeotiff/geo_normalize.c |   59 ++++++++++++++++++++++++++-
 2 files changed, 80 insertions(+), 7 deletions(-)

diff --git a/gdal/frmts/gtiff/gt_wkt_srs.cpp b/gdal/frmts/gtiff/gt_wkt_srs.cpp
index 5934df5..0afb101 100644
--- a/gdal/frmts/gtiff/gt_wkt_srs.cpp
+++ b/gdal/frmts/gtiff/gt_wkt_srs.cpp
@@ -737,9 +737,21 @@ char *GTIFGetOGISDefn( GTIF *hGTIF, GTIFDefn * psDefn )
             break;
 
           case CT_Mercator:
-            oSRS.SetMercator( adfParm[0], adfParm[1],
-                              adfParm[4],
-                              adfParm[5], adfParm[6] );
+            /* If a lat_ts was specified use 2SP, otherwise use 1SP */
+            if (psDefn->ProjParmId[2] == ProjStdParallel1GeoKey)
+            {
+                if (psDefn->ProjParmId[4] == ProjScaleAtNatOriginGeoKey)
+                    CPLError( CE_Warning, CPLE_AppDefined,
+                              "Mercator projection should not define both StdParallel1 and ScaleAtNatOrigin.\n"
+                              "Using StdParallel1 and ignoring ScaleAtNatOrigin.\n" );
+                oSRS.SetMercator2SP( adfParm[2],
+                                     adfParm[0], adfParm[1],
+                                     adfParm[5], adfParm[6]);
+            }
+            else
+                oSRS.SetMercator( adfParm[0], adfParm[1],
+                                  adfParm[4],
+                                  adfParm[5], adfParm[6] );
                               
             if (psDefn->Projection == 1024 || psDefn->Projection == 9841) // override hack for google mercator. 
             {
@@ -1510,14 +1522,18 @@ int GTIFSetFromOGISDefn( GTIF * psGTIF, const char *pszOGCWKT )
         GTIFKeySet(psGTIF, ProjNatOriginLongGeoKey, TYPE_DOUBLE, 1,
                    poSRS->GetNormProjParm( SRS_PP_CENTRAL_MERIDIAN, 0.0 ) );
         
-        GTIFKeySet(psGTIF, ProjScaleAtNatOriginGeoKey, TYPE_DOUBLE, 1,
-                   poSRS->GetNormProjParm( SRS_PP_SCALE_FACTOR, 1.0 ) );
-        
         GTIFKeySet(psGTIF, ProjFalseEastingGeoKey, TYPE_DOUBLE, 1,
                    poSRS->GetProjParm( SRS_PP_FALSE_EASTING, 0.0 ) );
         
         GTIFKeySet(psGTIF, ProjFalseNorthingGeoKey, TYPE_DOUBLE, 1,
                    poSRS->GetProjParm( SRS_PP_FALSE_NORTHING, 0.0 ) );
+
+        if( EQUAL(pszProjection,SRS_PT_MERCATOR_2SP) )
+            GTIFKeySet(psGTIF, ProjStdParallel1GeoKey, TYPE_DOUBLE, 1,
+                       poSRS->GetNormProjParm( SRS_PP_STANDARD_PARALLEL_1, 0.0 ) );
+        else
+            GTIFKeySet(psGTIF, ProjScaleAtNatOriginGeoKey, TYPE_DOUBLE, 1,
+                       poSRS->GetNormProjParm( SRS_PP_SCALE_FACTOR, 1.0 ) );
     }
     
     else if( EQUAL(pszProjection,SRS_PT_OBLIQUE_STEREOGRAPHIC) )
diff --git a/gdal/frmts/gtiff/libgeotiff/geo_normalize.c b/gdal/frmts/gtiff/libgeotiff/geo_normalize.c
index 4761ea0..75705ee 100644
--- a/gdal/frmts/gtiff/libgeotiff/geo_normalize.c
+++ b/gdal/frmts/gtiff/libgeotiff/geo_normalize.c
@@ -1495,8 +1495,65 @@ static void GTIFFetchProjParms( GTIF * psGTIF, GTIFDefn * psDefn )
         break;
 
 /* -------------------------------------------------------------------- */
-      case CT_LambertConfConic_1SP:
       case CT_Mercator:
+/* -------------------------------------------------------------------- */
+        if( GTIFKeyGet(psGTIF, ProjNatOriginLongGeoKey,
+                       &dfNatOriginLong, 0, 1 ) == 0
+            && GTIFKeyGet(psGTIF, ProjFalseOriginLongGeoKey,
+                          &dfNatOriginLong, 0, 1 ) == 0
+            && GTIFKeyGet(psGTIF, ProjCenterLongGeoKey,
+                          &dfNatOriginLong, 0, 1 ) == 0 )
+            dfNatOriginLong = 0.0;
+
+        if( GTIFKeyGet(psGTIF, ProjNatOriginLatGeoKey,
+                       &dfNatOriginLat, 0, 1 ) == 0
+            && GTIFKeyGet(psGTIF, ProjFalseOriginLatGeoKey,
+                          &dfNatOriginLat, 0, 1 ) == 0
+            && GTIFKeyGet(psGTIF, ProjCenterLatGeoKey,
+                          &dfNatOriginLat, 0, 1 ) == 0 )
+            dfNatOriginLat = 0.0;
+
+
+        if( GTIFKeyGet(psGTIF, ProjStdParallel1GeoKey,
+                       &dfStdParallel1, 0, 1 ) == 0 )
+            dfStdParallel1 = NAN;
+
+        if( GTIFKeyGet(psGTIF, ProjScaleAtNatOriginGeoKey,
+                       &dfNatOriginScale, 0, 1 ) == 0 )
+        {
+            /* Default only if dfStdParallel1 isn't defined */
+            if( isnan(dfStdParallel1) )
+                dfNatOriginScale = 1.0;
+            else
+                dfNatOriginScale = NAN;
+        }
+
+        /* notdef: should transform to decimal degrees at this point */
+
+        psDefn->ProjParm[0] = dfNatOriginLat;
+        psDefn->ProjParmId[0] = ProjNatOriginLatGeoKey;
+        psDefn->ProjParm[1] = dfNatOriginLong;
+        psDefn->ProjParmId[1] = ProjNatOriginLongGeoKey;
+        if( !isnan(dfStdParallel1) )
+        {
+            psDefn->ProjParm[2] = dfStdParallel1;
+            psDefn->ProjParmId[2] = ProjStdParallel1GeoKey;
+        }
+        if( !isnan(dfNatOriginScale) )
+        {
+            psDefn->ProjParm[4] = dfNatOriginScale;
+            psDefn->ProjParmId[4] = ProjScaleAtNatOriginGeoKey;
+        }
+        psDefn->ProjParm[5] = dfFalseEasting;
+        psDefn->ProjParmId[5] = ProjFalseEastingGeoKey;
+        psDefn->ProjParm[6] = dfFalseNorthing;
+        psDefn->ProjParmId[6] = ProjFalseNorthingGeoKey;
+
+        psDefn->nParms = 7;
+        break;
+
+/* -------------------------------------------------------------------- */
+      case CT_LambertConfConic_1SP:
       case CT_ObliqueStereographic:
       case CT_TransverseMercator:
       case CT_TransvMercator_SouthOriented:
-- 
1.7.10.4