[gdal-dev] gdaltransform and cs2cs don't agree ...

[Even Rouault]

Selon Jean-Claude Repetto <jrepetto at free.fr>:

Jean-Claude,

I've investigated a bit. Here are my findings.

gdaltransform asks PROJ.4 to expand the SRS definition "+init=IGNF:MILLER" into
a full PROJ.4 string.

The lookup in the IGNF file results in :
+proj=mill +towgs84=0.0000,0.0000,0.0000,0.0000,0.0000,0.0000,0.000000
+a=6378137.0000 +rf=298.2572221010000 +lon_0=0.000000000 +x_0=0.000 +y_0=0.000
+units=m +no_defs

Then OSR interprets the string to build its WKT representation, which leads to :
PROJCS["unnamed",
    GEOGCS["unnamed ellipse",
        DATUM["unknown",
            SPHEROID["unnamed",6378137,298.257222101],
            TOWGS84[0,0,0,0,0,0,0]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433]],
    PROJECTION["Miller_Cylindrical"],
    PARAMETER["latitude_of_center",0],
    PARAMETER["longitude_of_center",0],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0],
    UNIT["Meter",1]]

Which is OK up to that point.

Now when building the transformation between the 2 SRS, OSR builds a new PROJ.4
string from the WKT (not using the original IGNF:MILLER expansion) and here what
we get :
+proj=mill +lat_0=0 +lon_0=0 +x_0=0 +y_0=0 +R_A +ellps=GRS80
+towgs84=0,0,0,0,0,0,0 +units=m +no_defs

The difference (+a=6378137.0000 +rf=298.2572221010000) vs +ellps=GRS80 is OK
since those are 2 equivalent definitions.

The issue is the added +R_A which is not found in the original IGNF file.

Now if I test(note the removed +R_A and the added +wktext) :

gdaltransform -s_srs "+proj=mill +lat_0=0 +lon_0=0 +x_0=0 +y_0=0 +ellps=GRS80
+towgs84=0,0,0,0,0,0,0 +units=m +no_defs +wktext" -t_srs EPSG:4326

I get the same results as in the cs2cs command line.

I have no idea why we add +R_A when writting the PROJ.4 string of a Miller
projection. SVN history points to a commit that dates back to 10 years ago :
http://trac.osgeo.org/gdal/changeset/1862/trunk/ogr/ogr_srs_proj4.cpp

Perhaps Frank will remember... ?

In any case, you can open a Trac ticket to record the issue.

Best regards,

Even


> Hello,
> I have found a case where gdaltransform and cs2cs don't give the same
> results :
>
> $ gdaltransform -s_srs +init=IGNF:MILLER -t_srs EPSG:4326
> -20037504 10018752
> 179.798600354466 72.943850852954 12400.9953186931
>
> $ cs2cs  -f "%.11f" +init=IGNF:MILLER +to +init=epsg:4326
> -20037504 10018752
> -179.99996098806	72.78284247593 0.00000000000
>
>
> Can somebody confirm the bug ? If it is not a bug, is there an explanation ?
>
> I am using PROJ-4.7.0 and GDAL-1.8.0.
>
> Regards,
> Jean-Claude
> _______________________________________________
> gdal-dev mailing list
> gdal-dev at lists.osgeo.org
> http://lists.osgeo.org/mailman/listinfo/gdal-dev
>