Mercator 1SP PROJ Translation Wrong

[warmerdam]

A file with this listgeo report:

Geotiff_Information:
   Version: 1
   Key_Revision: 1.0
   Tagged_Information:
      ModelTiepointTag (2,3):
         874              434              0                
         -1237.52006      4190457.16       0                
      ModelPixelScaleTag (1,3):
         -1237.52006      1237.52006       0                
      End_Of_Tags.
   Keyed_Information:
      GTModelTypeGeoKey (Short,1): ModelTypeProjected
      GTRasterTypeGeoKey (Short,1): RasterPixelIsArea
      GeogGeodeticDatumGeoKey (Short,1): User-Defined
      GeogPrimeMeridianGeoKey (Short,1): PM_Greenwich
      GeogAngularUnitsGeoKey (Short,1): Angular_Degree
      GeogEllipsoidGeoKey (Short,1): User-Defined
      GeogSemiMajorAxisGeoKey (Double,1): 6378137          
      GeogSemiMinorAxisGeoKey (Double,1): 6356752.31       
      ProjectedCSTypeGeoKey (Short,1): User-Defined
      PCSCitationGeoKey (Ascii,14): "METOC-Halifax"
      ProjectionGeoKey (Short,1): User-Defined
      ProjCoordTransGeoKey (Short,1): CT_Mercator
      ProjLinearUnitsGeoKey (Short,1): Linear_Meter
      ProjNatOriginLongGeoKey (Double,1): -59.5            
      ProjNatOriginLatGeoKey (Double,1): 39               
      ProjFalseEastingGeoKey (Double,1): 0                
      ProjFalseNorthingGeoKey (Double,1): 0                
      End_Of_Keys.
   End_Of_Geotiff.

Projection Method: CT_Mercator
   ProjNatOriginLatGeoKey: 39.000000 ( 39d 0' 0.00"N)
   ProjNatOriginLongGeoKey: -59.500000 ( 59d30' 0.00"W)
   ProjScaleAtNatOriginGeoKey: 1.000000
   ProjFalseEastingGeoKey: 0.000000 m
   ProjFalseNorthingGeoKey: 0.000000 m
Prime Meridian: 8901/Greenwich (0.000000/  0d 0' 0.00"E)
Projection Linear Units: 9001/metre (1.000000m)

Corner Coordinates:
Upper Left    ( 1080355.009, 4727540.859)  ( 47d 1'42.86"W, 48d 0' 1.17"N)
Lower Left    ( 1080355.009, 3652135.931)  ( 47d 1'42.86"W, 38d59'19.87"N)
Upper Right   (-1084067.569, 4727540.859)  ( 72d 0'51.43"W, 48d 0' 1.17"N)
Lower Right   (-1084067.569, 3652135.931)  ( 72d 0'51.43"W, 38d59'19.87"N)
Center        (   -1856.280, 4189838.395)  ( 59d31'17.14"W, 43d39'52.33"N)

Produces this gdalinfo report:

Driver: GTiff/GeoTIFF
Size is 1749, 869
Coordinate System is:
PROJCS["unnamed",
    GEOGCS["unnamed",
        DATUM["unknown",
            SPHEROID["unnamed",6378137,298.257223560493]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433]],
    PROJECTION["Mercator_1SP"],
    PARAMETER["latitude_of_origin",39],
    PARAMETER["central_meridian",-59.5],
    PARAMETER["scale_factor",1],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]]]
Origin = (1080355.008799013914540,4727540.859266426414251)
Pixel Size = (-1237.520055898069131,-1237.520055898069359)
...
mapProjection = "METOC-Halifax"
mapUpperLeft = [48.0003,-72]
mapUpperRight = [48.0003,-47.0286]
mapLowerLeft = [39,-72]
mapLowerRight = [39,-47.0286]
...
Corner Coordinates:
Upper Left  ( 1080355.009, 4727540.859) ( 49d47'42.02"W, 39d13'45.43"N)
Lower Left  ( 1080355.009, 3652135.931) ( 49d47'42.02"W, 31d19'11.08"N)
Upper Right (-1084067.569, 4727540.859) ( 69d14'18.04"W, 39d13'45.43"N)
Lower Right (-1084067.569, 3652135.931) ( 69d14'18.04"W, 31d19'11.08"N)
Center      (   -1856.280, 4189838.395) ( 59d31'0.03"W, 35d22'19.83"N)

Note the different (wrong) corners from gdalinfo. It appears that listgeo is (correctly) treating the Geotiff ProjNatOriginLatGeoKey as a PROJ.4 "latitude of true scale" (lat_ts) in it's proj.4 translation, but GDAL is not.

The ​http://geotiff.maptools.org/proj_list/mercator_1sp.html page seems to indicate that the latitude_of_natural_origin has no proj.4 analog though it does mention lat_ts. Changing ogr_srs_proj4.cpp to set the value in lat_ts seems to produce the expected results.

[warmerdam]

After comparing the Mercator 1SP and Mercator 2SP cases, I am now hopelessly confused.

[rsignell]

Folks,

I think I see what the problem is here for Merc 1SP. The user specifies "lat_ts", which is easier than figuring out what "k_0" should be at the equator to make k equal to one at the latitude of true scale. Right now GDAL is saving lat_0=0, which is correct, but k_0=1, which is incorrect. So GDAL just needs to compute "k_0" when "lat_ts" is input. According to Snyder 1987, the formula for "k" for a spheroid is:

k=sqrt(1-(e*sin(lat))^2)/cos(lat)

so if k=1 at lat_ts, then

k_0=cos(lat)/sqrt(1-(e*sin(lat))^2)

Reference: Snyder 1987, page 44, equation [7-8] (page 53 of the PDF version - I've placed a copy at ​http://stellwagen.er.usgs.gov/rps/share/PP_1395_ocr.pdf)

Test: If I use the above formula on my dataset with +lat_ts=41.65 m with the Clark 1866 ellipsoid, with e=0.0822719, I get k_0=0.748337828698051. If I specify this to gdal_translate instead of +lat_ts=41.65, I get the correct result (to within 1 m) when I then use gdalwarp to warp the Mercator geotiff to a geographic geotiff.

CORRECT RESULT (1 STEP using +lat_ts):
gdalwarp -s_srs "+proj=merc +lat_ts=41.65 +lon_0=-70.3166667 +ellps=clrk66" -t_srs EPSG:4326 mb_backpc30m.tif geo1.tif

INCORRECT RESULT (2 STEP using +lat_ts):
gdal_translate -a_srs "+proj=merc +lat_ts=41.65 +lon_0=-70.3166667 +ellps=clrk66" mb_backpc30m.tif merc2.tif
gdalwarp -t_srs EPSG:4326 merc2.tif geo2.tif

CORRECT RESULT (2 STEP using +k_0)
$ gdal_translate -a_srs "+proj=merc +k_0=0.748337828698051 +lon_0=-70.3166667 +ellps=clrk66" mb_backpc30m.tif merc3.tif
gdalwarp -t_srs EPSG:4326 merc3.tif geo3.tif
$

Here's the proof:

$ gdalinfo geo1.tif | grep Upper
Upper Left  ( -70.9001641,  42.7991819) ( 70d54'0.59"W, 42d47'57.05"N)
Upper Right ( -70.0334114,  42.7991819) ( 70d 2'0.28"W, 42d47'57.05"N)

$ gdalinfo geo2.tif | grep Upper
Upper Left  ( -70.7533199,  33.4420528) ( 70d45'11.95"W, 33d26'31.39"N)
Upper Right ( -70.1045443,  33.4420528) ( 70d 6'16.36"W, 33d26'31.39"N)

$ gdalinfo geo3.tif | grep Upper
Upper Left  ( -70.9001640,  42.7991739) ( 70d54'0.59"W, 42d47'57.03"N)
Upper Right ( -70.0334114,  42.7991739) ( 70d 2'0.28"W, 42d47'57.03"N)

The latitudes from Methods 1 and 3 are within 0.02", or 0.6 meters.

-Rich

[warmerdam]

Rich / Allan,

On further review, i'm not becoming convinced that this file is actually Mercator 2SP and that the ProjNatOriginLatGeoKey value of 39 should be a ProjStdParallel1GeoKey instead. Furthermore, I've come to the conclusion that much of the geotiff code doesn't understand the difference between the 1SP and 2SP parameters.

I have a bunch of changes in the queue to address the flaws in libgeotiff and gdal around mercator 1 sp and mercator 2sp, but I thought I should give you a chance to comment before I do so.

In particular it seems that there is no way for me to know that for this file the natural origin lat was intended to be a standard parallel. So my changes won't help this file work properly. The user will just have to override the coordinate system.

I'm also going to follow up with an email to the geotiff list.

[warmerdam]

GeoTIFF list email:

Folks,

I was recently provided with a GeoTIFF file with this set of keys:

Geotiff_Information:
   Version: 1
   Key_Revision: 1.0
   Tagged_Information:
      ModelTiepointTag (2,3):
         874              434              0                
         -1237.52006      4190457.16       0                
      ModelPixelScaleTag (1,3):
         -1237.52006      1237.52006       0                
      End_Of_Tags.
   Keyed_Information:
      GTModelTypeGeoKey (Short,1): ModelTypeProjected
      GTRasterTypeGeoKey (Short,1): RasterPixelIsArea
      GeogGeodeticDatumGeoKey (Short,1): User-Defined
      GeogPrimeMeridianGeoKey (Short,1): PM_Greenwich
      GeogAngularUnitsGeoKey (Short,1): Angular_Degree
      GeogEllipsoidGeoKey (Short,1): User-Defined
      GeogSemiMajorAxisGeoKey (Double,1): 6378137          
      GeogSemiMinorAxisGeoKey (Double,1): 6356752.31       
      ProjectedCSTypeGeoKey (Short,1): User-Defined
      PCSCitationGeoKey (Ascii,14): "METOC-Halifax"
      ProjectionGeoKey (Short,1): User-Defined
      ProjCoordTransGeoKey (Short,1): CT_Mercator
      ProjLinearUnitsGeoKey (Short,1): Linear_Meter
      ProjNatOriginLongGeoKey (Double,1): -59.5            
      ProjNatOriginLatGeoKey (Double,1): 39               
      ProjFalseEastingGeoKey (Double,1): 0                
      ProjFalseNorthingGeoKey (Double,1): 0                
      End_Of_Keys.
   End_Of_Geotiff.

It seems it is intended to be a Mercator 2SP projection with 39 as 
the standard parallel (or latitude of true scale).  It seems to me that
encoding the standard parallel as NatOriginLat is an error. 

The mercator pages I have been maintaining seem quite mixed up on this
topic:

  http://geotiff.maptools.org/proj_list/mercator_1sp.html
  http://geotiff.maptools.org/proj_list/mercator_2sp.html

I'm seeking to clarify that Mercator 1SP would have the GeoTIFF 
parameters:

  CT_Mercator 
  NatOriginLong (aka Central Meridian)
  ScaleAtNatOrigin
  FalseEasting
  FalseNorthing

and that Mercator 2SP would have the geotiff parameters:

  CT_Mercator
  NatOriginLong 
  StdParallel1  (ie. latitude of true scale)
  FalseEasting
  FalseNorthing

Essentially the forms are equivelent as the 2SP version can be turned
into the 1SP version by computing the equivelent scale at the equator
that corresponds with having the described standard parallel being the
latitude of true scale. 
  
Matters have been muddied (in my mind at least) by the NatOriginLat 
parameter.  Is NatOriginLat intended to be the same as StdParalle1? 
Is it just another way of establishing what latitude is "zero northing"
(aka an alternate formulation of FalseNorthing)? Does it have some 
other meaning?

There does not seem to be any apparent mapping for NatOriginLat in
PROJ.4 or GCTP, so for my purposes if it is a distinct parameter I 
will just have to discard it since I have no way of using it. 

Input welcome.  

Retarget to 1.5. This is too scary to address in 1.4.

[rsignell]

Having reread Snyder's info on "MERCATOR PROJECTION WITH AN0THER STANDARD PARALLEL" (page 47, PDF page 56 using ​http://stellwagen.er.usgs.gov/rps/share/PP_1395_ocr.pdf) here's what I think:

The name Mercator 1SP should only be used when lat_ts=0 (the standard parallel is on the equator). This is the only situation where the Mercator projection has only one standard parallel. If the latitude of true scale is NOT on the equator, then there will be two standard parallels, at +lat_ts and at -lat_ts, which is Mercator 2SP. So for both Mercator 1SP and Mercator 2SP, lat_0=0, and

k_0=cos(lat_ts)/sqrt(1-(e*sin(lat_ts))^2)

At the equator, of course, lat_ts=0, and k_0=cos(0)=1.

I don't see why there is a distinction between Mercator 1SP and 2SP. Hey, come to think of it, there *isn't* a distinction in proj4.

[allanx13]

Frank,

A way to 'fix' the problem with the original image that I sent you has been found. gdalinfo reports that the x pixel size is -1237.520055898069131. Somebody here at DRDC found that by switching that value to the positive, and then running gdalwarp ("gdalwarp –t_srs WGS84 infile.tif outfile.tif"), the newly generated image displays correctly in OpenEV and other GDAL-based applications.

Apparently, PCI Geomatica and some other software ignore the negative sign in the x pixel size field and display the image properly; however, as I understand it, this goes against the GeoTIFF standard which supports negative pixel sizes. So, GDAL handles the negative value case properly, but, in doing so, ends up displaying the image improperly. Again, this *may* be a defect in the image itself, but we can't confirm that with the source of the imagery.

If anything here is unclear or if you have any questions, just let me know!

Allan

[allanx13]

In addition to my previous comment, we've now found that the solution works with GDAL 1.3.2, but not with 1.4.0. The 1.4.0 version of gdalwarp produces an image with incorrect corner coordinates.

[Even Rouault]

Removing obsolete milestone

[Even Rouault]

This issue appears to be solved. gdalinfo now reports the same thing as listgeo