[gdal-dev] Standard Parallel 1 - being changed but not displayed?
Jonathan Moules
jonathan-lists at lightpear.com
Wed Nov 18 04:53:04 PST 2015
Hi Even,Ok, that clears things up a little. I guess the real question I have at this point is - which is correct in this circumstance?
The original GeoTIFF?
ArcGIS?
GDAL's conversion?
None of the above? I get the feeling it's all a bit subjective.
>From your comments, I'm surprised this doesn't happen more often with various packages\files. I thought projections were nice and reasonably well defined things.
Thanks again,
Jonathan
---- On Tue, 17 Nov 2015 10:33:37 -0800 Even Rouault<even.rouault at spatialys.com> wrote ----
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
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/gdal-dev/attachments/20151118/6a88de1f/attachment.html>
More information about the gdal-dev
mailing list