I am trying to use Pyproj to perform the conversion between ocentric and ographic latitudes. I can do this using the cct
command line tool:
echo 0 45 | cct -z0 +proj=geoc +a=3396190 +b=3376200
0.0000000000 44.6617680466 0.0000 inf
based on the docs here.
When attempting to do this with Pyproj, I would expect to be able to define the projects and then apply the transformation. Something akin to:
from pyproj import Proj, transform
ocentric = Proj(proj='geoc', a=3396190, b=3376200)
ographic = Proj(proj='latlon', a=3396190, b=3376200)
ographic(0,45)
I can get the ographic projection to instantiate without an issue, but the ocentric projection fails with an invalid CRS error: CRSError: Input is not a CRS: +proj=geoc +a=3396190 +b=3376200 +type=crs
If you have been able to define a geocentric projection on a custom ellipsoid, then how did you do it?
Note: I also tried:
from pyproj import Proj, transform
ocentric = Proj(proj='latlon', geoc=1, a=3396190, b=3376200)
ographic = Proj(proj='latlon', a=3396190, b=3376200)
transform(ocentric, ographic, 0, 45)
and I get an error that the input is not a transformation. I believe that geoc
is being silently ignored in this example.
With help from @snowman2, I was able to define a pipeline to make the conversion work as follows.
trans = pyproj.transformer.Transformer.from_pipeline('+proj=pipeline +a=3396190 +b=3376200 +step +proj=geoc')
lon, lat = trans.transform(0, np.radians(45), errcheck=True)
np.degrees(lat) = #44.6617680466192
Note the weird np.radians call. When I run this transformation using 45 (as in degrees, which should be the default because radians=False is the transform default) I see: ProjError: transform error: latitude or longitude exceeded limits
.