[Proj] Convert a planetographic latitude to a planetocentric latitude

[Even Rouault]

> Good catch. I wonder why this was changed to the current implementation.
> The differing ellipsoid parameters for CRS1 and CRS2 can be regarded as
> datum information (albeit rudimentary datums).

I tracked down this to

commit d6636b805e0e8136f427ca79c57b21d421539169
Author: Frank Warmerdam <warmerdam at pobox.com>
Date:   Mon Nov 26 00:21:59 2007 +0000

    Modified PJ structure to hold a_orig, es_orig, ellipsoid definition before
    adjustment for spherical projections.
    Modified pj_datum_transform() to use the original ellipsoid parameters,
    not the ones adjusted for spherical projections.
    Modified pj_datum_transform() to not attempt any datum shift via
    geocentric coordinates if the source *or* destination are raw ellipsoids
    (ie. PJD_UNKNOWN).  All per PROJ bug #1602, GDAL bug #2025.
    
    
    git-svn-id: http://svn.osgeo.org/metacrs/proj/trunk@1406 
4e78687f-474d-0410-85f9-8d5e500ac6b2


The GDAL ticket is still there:
https://trac.osgeo.org/gdal/ticket/2025


There are 2 bugfixes mixed in that commit, but the specific bit for the issue 
discussed here is at the beginning of pj_datum_transform()

/* -------------------------------------------------------------------- */
/*      We cannot do any meaningful datum transformation if either      */
/*      the source or destination are of an unknown datum type          */
/*      (ie. only a +ellps declaration, no +datum).  This is new        */
/*      behavior for PROJ 4.6.0.                                        */
/* -------------------------------------------------------------------- */
    if( src->datum_type == PJD_UNKNOWN
        || dst->datum_type == PJD_UNKNOWN )
        return 0;


I guess the intent was that it would be obvious to users that something is 
wrong when they see that their coordinates do not change, whereas if silent
geocentric transformation was done, they would have the feeling that the
"right thing" is done, whereas they should probably have provided a horizontal 
grid or +towgs84 parameters to have a good result

Both behaviour can be discussed, but at this point of history of the cs2cs / 
pj_transform API, it might perhaps be prudent to freeze the behaviour of that 
legacy API to the current one (especially if we remove it at some future 
point)
I can recall a number of complaints/surprises when the 4.6.0 behaviour was 
introduced and people got different results than with previous versions. (That 
said, I guess they must have now added their +towgs84=0,0,0 parameter to get 
the pre-4.6 results, so if we change again, that shouldn't affect them.)

We can perhaps reconsider that for the proj_transform() / cct case.

Even

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