[Proj] Convert a planetographic latitude to a planetocentric latitude
even.rouault at spatialys.com
Sat Aug 18 15:18:08 EST 2018
On samedi 18 août 2018 21:30:22 CEST Jean-Christophe Malapert wrote:
> 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
> > 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
"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
$ 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.
Spatialys - Geospatial professional services
More information about the Proj