[gdal-dev] Standard Parallel 1 - being changed but not displayed?

Even Rouault even.rouault at spatialys.com
Tue Nov 17 10:33:57 PST 2015 

Le mardi 17 novembre 2015 19:13:52, Jonathan Moules a écrit :
> Hi Even,
> Thanks for the reply.
> 
> 
> I guess GDAL chose Mercator_1SP because that's the one that's explicitly
> defined in the projection text.

No, it is not explicitly defined, and GDAL ignores PCSCitationGeoKey to 
determine the projection method. It will use ProjCoordTransGeoKey = 
CT_Mercator for that, and in GDAL 1.11, it can only map it to Mercator_1SP

> 
> 
> I've had a quick look at the two files using listgeo, and highlighted the
> different lines with **** prefix below (it seems GDAL removes 6 tags, and
> adds one).
> 
> 
> The before:
>    Keyed_Information:
>       GTModelTypeGeoKey (Short,1): ModelTypeProjected
>       GTRasterTypeGeoKey (Short,1): RasterPixelIsArea
>       GTCitationGeoKey (Ascii,9): "Mercator"
>       GeographicTypeGeoKey (Short,1): GCS_WGS_84
>       GeogCitationGeoKey (Ascii,7): "WGS 84"
>   ****GeogGeodeticDatumGeoKey (Short,1): Datum_WGS84
>   ****GeogLinearUnitsGeoKey (Short,1): Linear_Meter
>       GeogAngularUnitsGeoKey (Short,1): Angular_Degree
>   ****GeogEllipsoidGeoKey (Short,1): Ellipse_WGS_84
>       GeogSemiMajorAxisGeoKey (Double,1): 6378137
>   ****GeogSemiMinorAxisGeoKey (Double,1): 6356752.31424518
>       GeogInvFlatteningGeoKey (Double,1): 298.257223563
>       ProjectedCSTypeGeoKey (Short,1): User-Defined
>   ****PCSCitationGeoKey (Ascii,9): "Mercator"
>       ProjectionGeoKey (Short,1): User-Defined
>       ProjCoordTransGeoKey (Short,1): CT_Mercator
>       ProjLinearUnitsGeoKey (Short,1): Linear_Meter
>   ****ProjStdParallel1GeoKey (Double,1): 60
>       ProjNatOriginLongGeoKey (Double,1): 0
>       ProjNatOriginLatGeoKey (Double,1): 0
>       ProjFalseEastingGeoKey (Double,1): 0
>       ProjFalseNorthingGeoKey (Double,1): 0
>       End_Of_Keys.
> 
> ------------
> The after:
>    Keyed_Information:
>       GTModelTypeGeoKey (Short,1): ModelTypeProjected
>       GTRasterTypeGeoKey (Short,1): RasterPixelIsArea
>       GTCitationGeoKey (Ascii,9): "Mercator"
>       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_Mercator
>       ProjLinearUnitsGeoKey (Short,1): Linear_Meter
>       ProjNatOriginLongGeoKey (Double,1): 0
>       ProjNatOriginLatGeoKey (Double,1): 0
>       ProjFalseEastingGeoKey (Double,1): 0
>       ProjFalseNorthingGeoKey (Double,1): 0
>   ****ProjScaleAtNatOriginGeoKey (Double,1): 1
>       End_Of_Keys.
>    End_Of_Geotiff.
> 
> 

Actually when trying that with GDAL 2.0, you would now get :

PROJCS["Mercator",
    GEOGCS["WGS 84",
        DATUM["WGS_1984",
            SPHEROID["WGS 84",6378137,298.257223563,
                AUTHORITY["EPSG","7030"]],
            AUTHORITY["EPSG","6326"]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433],
        AUTHORITY["EPSG","4326"]],
    PROJECTION["Mercator_2SP"],
    PARAMETER["standard_parallel_1",60],
    PARAMETER["central_meridian",0],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]]]

due to the changes done in https://trac.osgeo.org/gdal/ticket/5791 

The new logic to determine Mercator_1SP vs _2SP is best explainted by the code 
itself ;-)

{{{
            /* 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] );
}}}

That new behaviour is more consistant with the geotiff keys (and should be the 
one you get with ArcGIS initialy), although it perhaps reveals a problem with 
the encoding of the geotiff file itself if you say that the Mercator_1SP 
interpretation with scale=1 leads to correct geopositioning.

> 
> I'm not clear on why the things that changed were changed, but I can see
> that GDAL removed the ProjStdParallel1GeoKey and value. Should not that
> have been displayed by GDALSRSInfo for the Before file, even if it was
> wrongly set for this given projection?

GDAL is all about abstraction. It interprets at its best the information from 
the geotiff keys to a well-known projection method, and thus can ignore 
silently some keys that it doesn't need.

> 
> 
> Thanks,
> Jonathan
> 
> 
> ---- On Tue, 17 Nov 2015 09:18:49 -0800 Even
> Rouault<even.rouault at spatialys.com> wrote ----
> 
> Le mardi 17 novembre 2015 17:55:15, Jonathan Moules a écrit :
> > Hi List,
> > I have a Geotiff which includes this projection:
> >
> > PROJ.4 : '+proj=merc +lon_0=0 +k=1 +x_0=0 +y_0=0 +datum=WGS84 +units=m
> > +no_defs '
> >
> >
> > OGC WKT :
> > PROJCS["Mercator",
> > GEOGCS["WGS 84",
> > DATUM["WGS_1984",
> > SPHEROID["WGS 84",6378137,298.257223563,
> > AUTHORITY["EPSG","7030"]],
> > AUTHORITY["EPSG","6326"]],
> > PRIMEM["Greenwich",0],
> > UNIT["degree",0.0174532925199433],
> > AUTHORITY["EPSG","4326"]],
> > PROJECTION["Mercator_1SP"],
> > PARAMETER["central_meridian",0],
> > PARAMETER["scale_factor",1],
> > PARAMETER["false_easting",0],
> > PARAMETER["false_northing",0],
> > UNIT["metre",1,
> > AUTHORITY["EPSG","9001"]]]
> >
> >
> >
> >
> > If I load this raster into ArcGIS, it displays in the wrong place (a
> few > thousand kilometres North).
> >
> >
> > I then run it through gdal_translate (GDAL 1.11.1), with no flags:
> > * gdal_trainslate input.tif output.tif
> >
> > For output.tif, GDALSRSInfo shows that the projection is identical,
> but now > the file loads correctly in ArcGIS. The same file works fine
> in QGIS both > before and after the "translation".
> >
> >
> > Looking at the projection info in ArcGIS, it displays one difference:
> > Before (not working):
> > Standard_parallel_1 = 60
> >
> >
> > After (working):
> > Standard_parallel_1 = 0
> >
> >
> > But I don't see anything about those in either of the GDALSRSInfo
> outputs. >
> >
> > So my questions:
> > - What is gdal_translate doing to the file to "fix" it?
> > - If it is something to do with Standard Parallel 1 - why isn't this
> > component of the projection exposed by GDAL?
> 
> Yes, in Mercator_1SP, there's no Standard Parallel 1, this is for
> Mercator_2SP.
> 
> See
> http://www.remotesensing.org/geotiff/proj_list/mercator_1sp.html
> http://www.remotesensing.org/geotiff/proj_list/mercator_2sp.html
> 
> I guess your original geotiff file has some unusual formulation which is
> detected as Mercator_1SP by GDAL, and probably Mercator_2SP by ArCGIS.
> 
> You could try with the listgeo utility that comes with libgeotiff to
> display the geotiff keys.
> 
> >
> > Thoughts welcome.
> > Thanks,
> > Jonathan

-- 
Spatialys - Geospatial professional services
http://www.spatialys.com

More information about the gdal-dev mailing list