gs48 / alsk / gs50: `P->a` is reassigned in the projection setup but `P->ra` keeps the ellipsoid's reciprocal, so forward and inverse disagree

Open
  • Is something wrong?
#4,933 2 comments 0 reactions 0 assignees View on GitHub

Maintainers usually reply within 1 day

@manduinca is already working on this.

Since Nov 5, 2026.

  • #4941 by @manduinca — open

Assessment

Difficulty
3/5
Estimated time
1-2 days
Newbie friendliness
68/100
Issue type
Bug
Clarity
Clearly specified
Activity status
Active
Tech stack
cpp
Domain
backend-api-design, testing-qa

Research direction

Start with the parameter assignments in src/projections/mod_ster.cpp for gs48, alsk and gs50, then compare how nzmg.cpp resets derived ellipsoid parameters. Check the existing projection cases in test/gie/builtins.gie and add round-trip coverage using the default ellipsoid. Done means forward and inverse results agree for all three projections with the default ellipsoid and with the supplied ellipsoid or radius options.

Written by the indexing model from the issue text.

Description

bug
AI usage
  • Fable 5.1 supported the creation of this ticket.
  • This ticket has been (mostly) written by myself (human), and is aimed at being read by humans.
Example of problem
+proj=alsk                     fwd (140 W, 60 N) =   665040.800  -383489.801   inv -> -139.99988618  59.99995133
+proj=alsk +ellps=clrk66       fwd (140 W, 60 N) =   665040.800  -383489.801   inv -> -140.00000000  60.00000000
+proj=alsk +R=6370997          fwd (140 W, 60 N) =   662599.055  -382708.191   inv -> -140.00000000  60.00000000
+proj=gs50                     fwd (100 W, 40 N) =  1672352.402  -347454.748   inv ->  -99.99980068  39.99992783
+proj=gs50 +ellps=clrk66       fwd (100 W, 40 N) =  1672352.402  -347454.748   inv -> -100.00000000  40.00000000
+proj=gs48                     fwd (120 W, 45 N) = -1863484.184   894632.102   inv -> -119.97122854  44.99640670
+proj=gs48 +ellps=clrk66       fwd (120 W, 45 N) = -1863548.718   894402.079   inv -> -119.97095025  44.99637412
+proj=gs48 +R=6370997          fwd (120 W, 45 N) = -1857681.173   915305.567   inv -> -120.00000000  45.00000000

The misses are the radius ratios: for alsk and gs50 on the default GRS80 the setup switches to Clarke 1866
(a = 6378206.4, 6378137 / 6378206.4 = 0.999989), about 7 m and 17 m at those points; for gs48 the setup forces the
sphere a = 6370997 (6378137 / 6370997 = 1.00112), about 2.3 km. Passing +ellps=clrk66 to alsk/gs50, or
+R=6370997 to gs48, pre-sets the same a, and the round trips close to 1e-8 degrees. (echo -140 60 | proj +proj=alsk | invproj +proj=alsk -f %.8f shows that at the command line.)

Problem description

pj_init_ctx runs pj_calc_ellipsoid_params() (init.cpp:584), which sets P->ra = 1. / P->a (ell_set.cpp:618),
before the projection's own setup. PJ_PROJECTION(gs48) then assigns P->a = 6370997. (mod_ster.cpp:186), and alsk
and gs50 assign P->a = 6378206.4 or 6370997. (lines 214, 218, 250, 254) without touching P->ra. The forward
multiplies by P->a (fwd.cpp:142-143) and the inverse by P->ra (inv.cpp:92-93), so a round trip through any of the
three comes back scaled by a_old / a_new unless the user happened to pass the very radius the projection forces.
nzmg (nzmg.cpp:117, P->ra = 1. / (P->a = 6378388.0);) and healpix (healpix.cpp:667-668) do reset P->ra when
they change P->a; mod_ster.cpp seems to have forgotten this step.

Expected Output

Forward and inverse are each other's inverse for any +ellps/+R the user passes, as they are for nzmg.

Proposed fix

Reset the reciprocal wherever P->a is reassigned in mod_ster.cpp, as nzmg does:

    P->a = 6370997.;
    P->ra = 1. / P->a;

or recompute the derived parameters with pj_calc_ellipsoid_params(P, P->a, P->es) after the reassignment, which also covers P->b, P->rb, P->one_es and P->rone_es. (Note: three more sites for alsk/gs50's 6378206.4 and 6370997.)

A round-trip test for the three projections on the default ellipsoid in test/gie/builtins.gie
(they currently have forward expectations only, or expectations taken on the forced radius) would keep this from happening again.

Environment Information
  • PROJ version (proj): 9.7.0
  • Operation System Information: MacOS Sequoia 15.7.7
Installation method
  • conda, apt-get, from source, etc...
Dominant language
C++
Stars
2k
Forks
918
Avg merge
1d 11h
Merged PRs (30d)
20

Getting set up

First steps

  1. Read the whole issue, then the project's contributing guide.
  2. Comment on the issue to say you are picking it up — it saves two people doing the same work.
  3. Fork the repository and make your change on a branch.
  4. Open a pull request that references the issue number.

More from OSGeo/PROJ

All issues in OSGeo/PROJ

Similar issues

More C++ issues

Get new issues in your inbox

A short digest of beginner-friendly GitHub issues.