[Proj] Is it possible to store a WGS84 / Pseudo Mercator image as a GeoTIFF?

[Even Rouault]

Le dimanche 08 février 2015 11:26:35, Mikael Rittri a écrit :
> And if it is possible, since when or what version of
> GeoTIFF? And do most GeoTIFF readers support it?
> 
> (I apologize if this is the wrong mailing list for the questions.)
> 
> I am writing a report on Pseudo Mercator, and I realized
> that I don't know how well it is supported in GeoTIFF.

Mikael,

Would rather be a question the GeoTIFF list indeed.

Regarding libgeotiff, 
http://trac.osgeo.org/geotiff/changeset/1685 brought support for decoding 
EPSG:3857 in the "normalization" code (GTIFGetDefn()) 5 years ago.

But even for software using libgeotiff, it pretty much depends if they can 
interpret what libgeotiff returns and how the encoder encoded EPSG:3857.

For example, GDAL, when passed an explicit EPSG code basically only writes the 
EPSG code in the ProjectedCSTypeGeoKey key:

For example, if I do "gdal_translate foo.tif foo_3857.tif -a_srs EPSG:3857"

And then "listgeo foo_3857.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,25): "WGS 84 / Pseudo-Mercator"
      GeogCitationGeoKey (Ascii,7): "WGS 84"
      GeogAngularUnitsGeoKey (Short,1): Angular_Degree
      ProjectedCSTypeGeoKey (Short,1): Unknown-3857
      ProjLinearUnitsGeoKey (Short,1): Linear_Meter
      End_Of_Keys.
   End_Of_Geotiff.

PCS = 3857 (WGS 84 / Pseudo-Mercator)
Projection = 3856 (Popular Visualisation Pseudo-Mercator)
Projection Method: CT_Mercator
   ProjNatOriginLatGeoKey: 0.000000 (  0d 0' 0.00"N)
   ProjNatOriginLongGeoKey: 0.000000 (  0d 0' 0.00"E)
   ProjScaleAtNatOriginGeoKey: 1.000000
   ProjFalseEastingGeoKey: 0.000000 m
   ProjFalseNorthingGeoKey: 0.000000 m
GCS: 4326/WGS 84
Datum: 6326/World Geodetic System 1984
Ellipsoid: 7030/WGS 84 (6378137.00,6356752.31)
Prime Meridian: 8901/Greenwich (0.000000/  0d 0' 0.00"E)
Projection Linear Units: 9001/metre (1.000000m)

All the decoding of projection method and parameters is done by 
geo_normalize.c from the EPSG dictionnaries.

And even with that, it is not sufficient as we must fix the ellipsoid decoding 
from WGS84 to the sphere for proj.4 correct operation. So the GeoTIFF reader 
in GDAL has for example this extra hack after GTIFGetDefn():
{{{
            if (psDefn->Projection == 1024 || psDefn->Projection == 9841) // 
override hack for google mercator. 
            {
                oSRS.SetExtension( "PROJCS", "PROJ4",  
                                   "+proj=merc +a=6378137 +b=6378137 
+lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext  
+no_defs" ); 
            }
}}}


On the contrary, if I generate a GeoTIFF with the proj.4 string with:
gdal_translate foo.tif foo_pseudomercator.tif -a_srs "+proj=merc +a=6378137 
+b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m 
+nadgrids=@null +wktext  +no_defs"

listgeo foo_pseudomercator.tif gives:

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,8): "unnamed"
      GeographicTypeGeoKey (Short,1): User-Defined
      GeogCitationGeoKey (Ascii,83): "GCS Name = unnamed ellipse|Datum = 
unknown|Ellipsoid = unnamed|Primem = Greenwich|"
      GeogGeodeticDatumGeoKey (Short,1): User-Defined
      GeogAngularUnitsGeoKey (Short,1): Angular_Degree
      GeogEllipsoidGeoKey (Short,1): User-Defined
      GeogSemiMajorAxisGeoKey (Double,1): 6378137          
      GeogSemiMinorAxisGeoKey (Double,1): 6378137          
      GeogPrimeMeridianLongGeoKey (Double,1): 0                
      ProjectedCSTypeGeoKey (Short,1): User-Defined
      ProjectionGeoKey (Short,1): User-Defined
      ProjCoordTransGeoKey (Short,1): CT_Mercator
      ProjLinearUnitsGeoKey (Short,1): Linear_Meter
      ProjStdParallel1GeoKey (Double,1): 0                
      ProjNatOriginLongGeoKey (Double,1): 0                
      ProjNatOriginLatGeoKey (Double,1): 0                
      ProjFalseEastingGeoKey (Double,1): 0                
      ProjFalseNorthingGeoKey (Double,1): 0                
      End_Of_Keys.
   End_Of_Geotiff.

Projection Method: CT_Mercator
   ProjNatOriginLatGeoKey: 0.000000 (  0d 0' 0.00"N)
   ProjNatOriginLongGeoKey: 0.000000 (  0d 0' 0.00"E)
   ProjStdParallel1GeoKey: 0.000000 (  0d 0' 0.00"N)
   ProjFalseEastingGeoKey: 0.000000 m
   ProjFalseNorthingGeoKey: 0.000000 m
Projection Linear Units: 9001/metre (1.000000m)

On reading, GDAL will not really identify it as being PseudoMercator but the 
georeferencing seems to be OK.

Even



> 
> Best regards,
> 
> Mikael Rittri
> Carmenta
> Sweden
> http://www.carmenta.com
> _______________________________________________
> Proj mailing list
> Proj at lists.maptools.org
> http://lists.maptools.org/mailman/listinfo/proj

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