Erroneuous use of "latitude of origin" in Polar Stereographic projection

[desruisseaux]

This issue does not apply to GeoTIFF 1.1 specification, but it applies to the list of projection parameters that should be part of a GeoTIFF specification according issue #5. It may also apply to GDAL implementation (to be confirmed).

The http://geotiff.maptools.org/proj_list/ web site describes the mapping of GeoTIFF projection parameters to EPSG codes or PROJ arguments, as implemented by GDAL. We should consider putting this material, possibly with some corrections, in the GeoTIFF specification for avoiding ambiguities that reduce interoperability. However that material may need some corrections. An issue with Mercator projection was raised in above link. Another issue is described below.

The page about Polar Stereographic (GeoTIFF code 15: CT_PolarStereographic) maps that projection to EPSG::9810, which is Polar Stereographic (variant A). The EPSG registry lists the following parameters for that projection:

• Latitude of natural origin
• Longitude of natural origin
• Scale factor at natural origin
• False easting
• False northing

However, the EPSG guidance notes Coordinate Conversions and Transformations including Formulas — May 2022 said the following about the Latitude of natural origin parameter:

In the basic variant (variant A) the latitude of origin is either the north pole or the south pole, at which is defined a scale factor at the natural origin (…snip…) In variant A the parameter Latitude of natural origin is used only to identify which hemisphere case is required. The only valid entries are ±90° or equivalent in alternative angle units.

But the use of that parameter in geotiff.maptools.org, where they give an example setting that parameter to 71°N, contradicts the EPSG statement. I think that geotiff.maptools.org actually wants to set the Latitude of standard parallel, not the latitude of natural origin. That parameter is defined in EPSG Polar Stereographic (variant B) projection, not variant A. My impression seems confirmed by the example given by geotiff.maptools.org where they set the 71°N value to PROJ +lat_ts argument rather than +lat_0.

Proposed corrections

In geotiff.maptools.org

In geotiff.maptools.org web site, or in the equivalent material that would be put in OGC GeoTIFF specification:

• In "EPSG Code", replace 9810 by 9829.
• Replace the list of parameters by this one from EPSG registry:
  • Latitude of standard parallel
  • Longitude of origin
  • False easting
  • False northing
• Map Latitude of standard parallel parameter to ProjStdParallel1GeoKey (GeoTIFF code 3078) instead of ProjNatOriginLatGeoKey (GeoTIFF code 3081).
• In the example mapping to PROJ arguments, replace +lat_ts=Latitude at natural origin by +lat_ts=Latitude of standard parallel.

In GDAL

If geotiff.maptools.org reflects GDAL implementation (it seems to be the case for GDAL 3.4.3 according a test we did with gdal_translate), then the correction would be:

• If the ProjScaleAtNatOriginGeoKey has a value different than 1, keep current behaviour (see next section).
• Otherwise when writing the parameters of projection code 15:
  • Omit the ProjScaleAtNatOriginGeoKey parameter (code 3092).
  • Replace the parameter code 3081 by 3078 (associated to the same value).

In GeoTIFF specification

Proposed heuristic for reducing the risk of compatibility break:

• If the scale factor is different than 1, then GeoTIFF projection 15 shall be interpreted as Polar Stereographic (variant A) and the expected latitude parameter is ProjNatOriginLatGeoKey (code 3081). This is the current behaviour except that the latitude should be restricted to -90° or +90°.
• Otherwise the GeoTIFF projection 15 shall be interpreted as Polar Stereographic (variant B) and the expected latitude parameter is ProjStdParallel1GeoKey (code 3078). The ProjScaleAtNatOriginGeoKey (code 3092) parameter shall not be present.

[rouault]

I've reviewed things a bit. Both variant A and variant B are handled with the mapping of http://geotiff.maptools.org/proj_list/polar_stereographic.html. NatOriginLat is (ab)used to store the Latitude of standard parallel in the case of variant B. On reading, the criterion to decide on the variant is based on the ScaleAtNatOrigin. If it is 1, then this is variant B and NatOriginLat is used at the Latitude of standard parallel, otherwise it is variant A.

While the above proposal matches more the EPSG registry, it would cause compatibility issues with a > 10 year industry practice.

A less risky resolution would be to edit http://geotiff.maptools.org/proj_list/polar_stereographic.html to document the current mapping w.r.t the parameters of the EPSG registry

[rouault]

libgeotiff repository contains https://github.com/OSGeo/libgeotiff/blob/master/libgeotiff/test/data/polar_stereographic_variant_a.tif and https://github.com/OSGeo/libgeotiff/blob/master/libgeotiff/test/data/polar_stereographic_variant_b.tif samples files which are recognized as:

$ listgeo ~/libgeotiff/libgeotiff/test/data/polar_stereographic_variant_a.tif
Geotiff_Information:
   Version: 1
   Key_Revision: 1.0
   Tagged_Information:
      ModelTiepointTag (2,3):
         0                 0                 0                
         440720            3751320           0                
      ModelPixelScaleTag (1,3):
         60                60                0                
      End_Of_Tags.
   Keyed_Information:
      GTModelTypeGeoKey (Short,1): ModelTypeProjected
      GTRasterTypeGeoKey (Short,1): RasterPixelIsArea
      GTCitationGeoKey (Ascii,21): "RSRGD2000 / RSPS2000"
      GeographicTypeGeoKey (Short,1): Code-4764 (RSRGD2000)
      GeogCitationGeoKey (Ascii,10): "RSRGD2000"
      GeogAngularUnitsGeoKey (Short,1): Angular_Degree
      GeogSemiMajorAxisGeoKey (Double,1): 6378137          
      GeogInvFlatteningGeoKey (Double,1): 298.257222101    
      ProjectedCSTypeGeoKey (Short,1): User-Defined
      ProjectionGeoKey (Short,1): User-Defined
      ProjCoordTransGeoKey (Short,1): CT_PolarStereographic
      ProjLinearUnitsGeoKey (Short,1): Linear_Meter
      ProjNatOriginLatGeoKey (Double,1): -90              
      ProjFalseEastingGeoKey (Double,1): 5000000          
      ProjFalseNorthingGeoKey (Double,1): 1000000          
      ProjScaleAtNatOriginGeoKey (Double,1): 0.994            
      ProjStraightVertPoleLongGeoKey (Double,1): 180              
      End_Of_Keys.
   End_Of_Geotiff.

$ gdalinfo ~/libgeotiff/libgeotiff/test/data/polar_stereographic_variant_a.tif
Driver: GTiff/GeoTIFF
Files: /home/even/libgeotiff/libgeotiff/test/data/polar_stereographic_variant_a.tif
Size is 1, 1
Coordinate System is:
PROJCRS["RSRGD2000 / RSPS2000",
    BASEGEOGCRS["RSRGD2000",
        DATUM["Ross Sea Region Geodetic Datum 2000",
            ELLIPSOID["GRS 1980",6378137,298.257222101004,
                LENGTHUNIT["metre",1]]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        ID["EPSG",4764]],
    CONVERSION["Polar Stereographic (variant A)",
        METHOD["Polar Stereographic (variant A)",
            ID["EPSG",9810]],
        PARAMETER["Latitude of natural origin",-90,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8801]],
        PARAMETER["Longitude of natural origin",180,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["Scale factor at natural origin",0.994,
            SCALEUNIT["unity",1],
            ID["EPSG",8805]],
        PARAMETER["False easting",5000000,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",1000000,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["(E)",north,
            MERIDIAN[90,
                ANGLEUNIT["degree",0.0174532925199433,
                    ID["EPSG",9122]]],
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["(N)",north,
            MERIDIAN[0,
                ANGLEUNIT["degree",0.0174532925199433,
                    ID["EPSG",9122]]],
            ORDER[2],
            LENGTHUNIT["metre",1]]]

$ listgeo ~/libgeotiff/libgeotiff/test/data/polar_stereographic_variant_b.tif
Geotiff_Information:
   Version: 1
   Key_Revision: 1.0
   Tagged_Information:
      ModelTiepointTag (2,3):
         0                 0                 0                
         440720            3751320           0                
      ModelPixelScaleTag (1,3):
         60                60                0                
      End_Of_Tags.
   Keyed_Information:
      GTModelTypeGeoKey (Short,1): ModelTypeProjected
      GTRasterTypeGeoKey (Short,1): RasterPixelIsArea
      GTCitationGeoKey (Ascii,50): "WGS 84 / Australian Antarctic Polar Stereographic"
      GeographicTypeGeoKey (Short,1): GCS_WGS_84
      GeogCitationGeoKey (Ascii,7): "WGS 84"
      GeogAngularUnitsGeoKey (Short,1): Angular_Degree
      GeogSemiMajorAxisGeoKey (Double,1): 6378137          
      GeogInvFlatteningGeoKey (Double,1): 298.257223563    
      ProjectedCSTypeGeoKey (Short,1): User-Defined
      ProjectionGeoKey (Short,1): User-Defined
      ProjCoordTransGeoKey (Short,1): CT_PolarStereographic
      ProjLinearUnitsGeoKey (Short,1): Linear_Meter
      ProjNatOriginLatGeoKey (Double,1): -71              
      ProjFalseEastingGeoKey (Double,1): 6000000          
      ProjFalseNorthingGeoKey (Double,1): 6000000          
      ProjScaleAtNatOriginGeoKey (Double,1): 1                
      ProjStraightVertPoleLongGeoKey (Double,1): 70               
      End_Of_Keys.
   End_Of_Geotiff.

$ gdalinfo ~/libgeotiff/libgeotiff/test/data/polar_stereographic_variant_b.tif
Driver: GTiff/GeoTIFF
Files: /home/even/libgeotiff/libgeotiff/test/data/polar_stereographic_variant_b.tif
Size is 1, 1
Coordinate System is:
PROJCRS["WGS 84 / Australian Antarctic Polar Stereographic",
    BASEGEOGCRS["WGS 84",
        DATUM["World Geodetic System 1984",
            ELLIPSOID["WGS 84",6378137,298.257223563,
                LENGTHUNIT["metre",1]]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433]],
        ID["EPSG",4326]],
    CONVERSION["Polar Stereographic (variant B)",
        METHOD["Polar Stereographic (variant B)",
            ID["EPSG",9829]],
        PARAMETER["Latitude of standard parallel",-71,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8832]],
        PARAMETER["Longitude of origin",70,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8833]],
        PARAMETER["False easting",6000000,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",6000000,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["(E)",north,
            MERIDIAN[90,
                ANGLEUNIT["degree",0.0174532925199433,
                    ID["EPSG",9122]]],
            ORDER[1],
            LENGTHUNIT["metre",1]],
        AXIS["(N)",north,
            MERIDIAN[0,
                ANGLEUNIT["degree",0.0174532925199433,
                    ID["EPSG",9122]]],
            ORDER[2],
            LENGTHUNIT["metre",1]]]