1

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.

0

1 Answer 1

3

Based on the input from @snowman2 I was able to find a solution:

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

The inverse is then:

itrans = pyproj.transformer.Transformer.from_pipeline('+proj=pipeline +a=3396190 +b=3376200 +step +proj=geoc +inv')
nlon, nlat = itrans.transform(lon, lat, errcheck=True, radians=True)
print(nlon, np.degrees(nlat))

The unit conversion is expected because of how the pipeline works. The proj pipeline documenation was really good for helping with the syntax of the Proj string.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.