[Proj] Convert a planetographic latitude to a planetocentric latitude

[Even Rouault]

On samedi 18 août 2018 21:30:22 CEST Jean-Christophe Malapert wrote:
> Hello,
> 
> I am quite new to the use of  Proj4. I spend a few days to investigate a
> problem but I have no clue how to solve it. So that's why I am writing to
> you.  I hope you will be able to answer me.
> 
> I have 2 different CRS for the same planet : one is based on the sphere
> (CRS1, radius=3396190) and another one based on the ellipse (CRS2,
> a=3396190 b=3376200). When I give the (longitude, latitude) in the CRS1 as
> inputs ( planetocentric latitude = planetographic latitude), I want to get
> the planetographic latitude expressed in the CRS2. So, I wrote the
> following command :
> cs2cs +proj=latlong +a=3396190 +b=3396190 +to +proj=latlong +a=3396190
> +b=3376200
> 
> > 0 80  # longitude latitude
> 
> I get the following result "0dE 80dN 0.000" . The latitude should change
> but it is not the case !!! Should I add some options ?  Is it possible to
> do this kind of transformation with proj4 ?

You are running into the case 3) of
https://proj4.org/faq.html#how-do-i-debug-problems-with-nad27-nad83-datum-shifting

"The cs2cs command and the underlying pj_transform() API know how to do a grid
shift as part of a more complex coordinate transformation; however, it is
imperative that both the source and destination coordinate system be defined
with appropriate datum information. [...] If either the input or output are not"
"identified as having a datum, the datum shifting (and *ellipsoid change*) step "
"is just quietly skipped!

The grid shift mentions do not apply to your case, but the mention of no ellipsoid change
does apply since you've only described ellipsoids with +a and +b
So basically PROJ will just reuse the geodetic coordinates unchanged

The fix is simple, although a bit counterintuitive when operating on non-Earth bodies,
add a +towgs84=0,0,0 clause to each CRS that will cause a transform from geodetic to
geocentric coordinates.

$ echo "0 80 0" | cs2cs -f "%.8f"  +proj=latlong +a=3396190 +b=3396190 +towgs84=0,0,0 +to +proj=latlong +a=3396190 +b=3376200 +towgs84=0,0,0
0.00000000	80.11439544 19392.34167562

With PROJ 5 or later, you can also use the cct utility with:

$ echo "0 80 0" | cct +proj=pipeline +step +proj=cart +a=3396190 +b=3396190 +step +inv +proj=cart +a=3396190 +b=3376200
  0.0000000000   80.1143954395    19392.3417           inf

Here you create a pipeline that explicitly transforms from geodetic coordinates on CS1 to
geocentric, and then from those geocentric to CS2. 

Even

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