0

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

Has anyone been able to define a geocentric projection on a custom ellipsoid? If so, how?

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.

Update

I was able to define a pipeline to make the conversion work as follows. Thanks @snowman2.

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.

3
  • 1
    Did you try Transformer.from_pipeline? – snowman2 Jul 3 '20 at 2:18
  • @snowman2 Updated the post. The pipeline works, but the arguments are not making sense when compared to the doc string. Specifically, the yy arg should default to degrees, but it appears to be in radians. Setting radians=False does not appear to have any effect. – Jay Laura Jul 3 '20 at 2:53
  • 1
    No unit conversions happen when using the pipeline in the current version. When version 3 comes out, the unit conversions will happen automatically. – snowman2 Jul 3 '20 at 12:22
1

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, privacy policy and cookie policy

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