[PROJ] Cassini-Soldner ellipsoid form bug ?

[Even Rouault]

Michael,

digging a bit, I see that the formulas that PROJ uses come from Snyder, 
which are also the same used in EPSG guidance note 7-2 
(https://www.iogp.org/wp-content/uploads/2019/09/373-07-02.pdf), use a 
truncated series expansions that must be valid in a (large, but not 
global) restricted domain around the projection center (here 0,0), and 
diverge beyond.  The formulas used by GeographicLib in 
https://github.com/geographiclib/geographiclib/blob/main/src/CassiniSoldner.cpp#L34 
are totally different looking, using geodesics computation, and I 
presume implement Charles' remark in https://arxiv.org/pdf/1109.4448 
"Two geodesic projections, the azimuthal equidistant projection and the 
Cassini-Soldner projection, are simple to write and their domain of 
applicability is not artificially restricted, as would be the case, for 
example, if the series expansion for the Cassini-Soldner projection were 
used (Snyder, 1987, §13)". And looking into Snyder's Map Projections, 
"For the ellipsoidal form, a set of series approximations is given for 
use in a zone extending 3° to 4° of longitude from the central 
meridian". So given that's the ones used by PROJ, no wonder that at 149° 
you get garbage. I suppose someone could port GeographicLib's formulas 
to PROJ (interestingly the inverse method using the geodesics is a 
2-liner and doesn't require iterative solving like done in PROJ), 
assuming PROJ's geodesic.c contains all the needed base functions

Even

Le 23/01/2026 à 03:19, Michael Sumner via PROJ a écrit :
> I'm having trouble with Cassini-Soldner:
>
> # Good point - agreeable coords
> echo "-6 -71" | proj +proj=cass +R=6378137        # spherical 
>  -217097.18 -7914445.72
> echo "-6 -71" | proj +proj=cass +ellps=WGS84      # ellipsoidal 
> -217749.76 -7891343.23
>
> # Bad point - ellipsoidal diverges
> echo "149 45" | proj +proj=cass +R=6378137        # spherical: 
> 2377511.89 14538559.49
> echo "149 45" | proj +proj=cass +ellps=WGS84      # ellipsoidal: 
> -2738664.72 28023005.14 (WRONG)
>
> That second coordinate seems wildly wrong, I've compared with 
> GeographicLib which shows
>
> lon,lat,x,y
> -6.000000,-71.000000,-217749.757184,-7891343.215151
> 149.000000,45.000000,2381489.299010,14528505.93346
>
> I've laid out a longer description exploring this (and includes the 
> C++ wrapper for GeographicLib):
>
> https://github.com/mdsumner/projcassini
>
> I hope that's enough to show the problem, I found this very confusing 
> but the output of GeographicLib makes me more confident.
>
> Cheers, Mike
>
> -- 
> Michael Sumner
> Research Software Engineer
> Australian Antarctic Division
> Hobart, Australia
> 0438489030
> e-mail: mdsumner at gmail.com
>
> _______________________________________________
> PROJ mailing list
> PROJ at lists.osgeo.org
> https://lists.osgeo.org/mailman/listinfo/proj

-- 
http://www.spatialys.com
My software is free, but my time generally not.
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/proj/attachments/20260123/9a36a418/attachment.htm>