[Geotiff] polar-stereo, proj4, and geotiff encoding...

[brownrigg at sunflower.com]

The so-called "normalized form"  specifying polar-stereographic project in geotiff seems to be incorrect.  Indeed, the "Random Outstanding Issues" page indicates concerns over the meaning of proj4 parameters and their mapping onto geotif tags. (http://www.remotesensing.org/geotiff/proj_list/random_issues.html)

I'm referring to in particular the meaning and interplay of +lat_0, +lat_ts, and +k_0.  From reading through John Synder's text and Gerald Evenden's libproj docs, +lat_0, along with +lon_0, establishes the origin of the projected space. +k_0 specifies a scaling factor that generally applies for stereographic projections, but in the special case of polar-stereographic projections, specifying a +lat_ts accomplishes the same thing. Indeed, the PROJ4 library ignores +k_0 if +lat_ts is present. The geometrical interpretation of lat_ts is directly analogous to standard parallels; it is the latitude at which the projection-plane intersects the Earth. I've seen specs where the scale factor is given, and others where the latitude-of-true scale is given;  either one does the job, but calculating one from the other is anything but straightforward.

Now consider the geotiff mapping. If lat_0 is determined to be a polar value, +lat_ts gets (erroneously) mapped onto +lat_0, as seen in ogr_srs_proj4.cpp in GDAL:

(lines 484-493, v. 1.6.0)

    else if( EQUAL(pszProj,"stere") 
             && ABS(OSR_GDV( papszNV, "lat_0", 0.0 ) + 90) < 0.001 )
    {
        SetPS( OSR_GDV( papszNV, "lat_ts", -90.0 ), 
               OSR_GDV( papszNV, "lon_0", 0.0 ), 
               OSR_GDV( papszNV, "k", 1.0 ), 
               OSR_GDV( papszNV, "x_0", 0.0 ), 
               OSR_GDV( papszNV, "y_0", 0.0 ) );
    }

and is reflected in GDAL's .../frmts/gtiff/libgeotiff/geo_normalize.c

(line 1227)
      case CT_LambertConfConic_1SP:
      case CT_Mercator:
      case CT_ObliqueStereographic:
      case CT_PolarStereographic:
      case CT_TransverseMercator:
      case CT_TransvMercator_SouthOriented:
        panProjParmId[0] = ProjNatOriginLatGeoKey;
        panProjParmId[1] = ProjNatOriginLongGeoKey;
        panProjParmId[4] = ProjScaleAtNatOriginGeoKey;
        panProjParmId[5] = ProjFalseEastingGeoKey;
        panProjParmId[6] = ProjFalseNorthingGeoKey;

Thus for our polar projection, the origin is no longer at the pole, any scaling intended via +lat_ts is now lost, and any scaling given via +k_0 will be based off the unintended, non-polar origin. According to Synder & Evenden, the following should be equivalent, but this geotiff mapping effectively guarantees they aren't:

     +proj=ups +south   /* Universal Polar stereographic,  defined in PROJ4 */
     +proj=stere +lat_0=-90 +k_0=.994   /* see Synder, pg 157 */
     +proj=stere +lat_0=-90 +lat_ts=-81.11453  /* Synder, pg 157 */

It seems to me that the geotiff mapping requires carrying along both +lat_0 and +lat_ts (where +lat_ts defaults to +lat_0 if not specified), especially if coordinate transformations are ultimately being handed off to PROJ4. 

FWIW...
Rick Brownrigg