PROJ Cassini-Soldner Ellipsoidal Forward Projection Bug?

Summary

PROJ’s ellipsoidal Cassini-Soldner (+proj=cass) forward projection produces incorrect results for points far from the central meridian. The spherical formulation works correctly, as does the ellipsoidal implementation in GeographicLib.

Environment

• PROJ build 2026-01-23 6a529ad433d87742ba7163e94fa2fdf043432db1
• GeographicLib (C++ library, used as ground truth)
• Tested January 2026

Test Cases

Four test points (lon, lat), two that work correctly and two that fail:

Point	Longitude	Latitude	Description
Good 1	-6	-71	Near central meridian
Good 2	-17	20	Near central meridian
Bad 1	146	-44	~150° from CM
Bad 2	149	45	~150° from CM

Results

Good point (-6, -71) — all implementations agree:

Source	x	y
GeographicLib (ellipsoidal)	-217749.76	-7891343.22
PROJ spherical	-217097.18	-7914445.72
PROJ ellipsoidal	-217749.76	-7891343.23

Bad point (149, 45) — PROJ ellipsoidal diverges catastrophically:

Source	x	y
GeographicLib (ellipsoidal)	2381489.30	14528505.93
PROJ spherical	2377511.89	14538559.49
PROJ ellipsoidal	-2738664.72	28023005.14

The PROJ ellipsoidal result has the wrong sign on x and nearly double the y value.

Note that PROJ’s spherical result is closer to the correct ellipsoidal answer than PROJ’s own ellipsoidal implementation.

Test coordinates (lon, lat)

test_pts <- matrix(c(-6, -17, 146, 149, -71, 20, -44, 45), ncol = 2)
colnames(test_pts) <- c("lon", "lat")

cass_ellps <- "+proj=cass +ellps=WGS84"
cass_spher <- "+proj=cass +R=6378137"
library(gdalraster)
#> GDAL 3.13.0dev-0e4f0cadf6 (released 2026-01-20), GEOS 3.12.1, PROJ 9.7.0
library(geographiclib) # gh: hypertidy/geographiclib
pts <- do.call(cbind, maps::map(plot = F)[1:2])

## these points are wrong
plot(ellps_xy <- transform_xy(test_pts, srs_from = "EPSG:4326", srs_to = cass_ellps), 
     xlim = c(-5000000,  5000000), ylim = c( -56460000,  41400000), cex  = 4, pch = 4, 
     xlab = "", ylab = "")

## red + is correct
points(karney_xy <- cassini_fwd(test_pts, lon0 = 0, lat0= 0)[1:2], col = "red", cex  = 2, pch = 3)
points(spher_xy <- transform_xy(test_pts, srs_from = "EPSG:4326", srs_to = cass_spher), pch = 19)
legend("bottomright", 
       legend = c("PROJ ellipsoidal (wrong)", "GeographicLib (correct)", "PROJ spherical"),
       pch = c(4, 3, 19),
       col = c("black", "red", "black"),
       pt.cex = c(2, 1.2, 1))

Reproduction

PROJ CLI

# 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)

GeographicLib C++ (ground truth)

#include <GeographicLib/CassiniSoldner.hpp>
#include <iostream>
#include <iomanip>

int main() {
    using namespace GeographicLib;
    
    const Geodesic geod(Constants::WGS84_a(), Constants::WGS84_f());
    CassiniSoldner cass(0.0, 0.0, geod);
    
    double test[][2] = {
        {-6, -71}, {-17, 20},
        {146, -44}, {149, 45}
    };
    
    std::cout << std::fixed << std::setprecision(6);
    std::cout << "lon,lat,x,y\n";
    
    for (auto& pt : test) {
        double x, y;
        cass.Forward(pt[1], pt[0], x, y);
        std::cout << pt[0] << "," << pt[1] << "," << x << "," << y << "\n";
    }
    
    return 0;
}

Compile: g++ -o cassini_test cassini_test.cpp -lGeographic

writeLines(system("./cassini_test", intern = TRUE))
#> lon,lat,x,y
#> -6.000000,-71.000000,-217749.757184,-7891343.215151
#> -17.000000,20.000000,-1775850.311294,2305579.711463
#> 146.000000,-44.000000,2644630.861016,-14532990.972920
#> 149.000000,45.000000,2381489.299010,14528505.933463

Related Issues

• #4385 / #4386 — Fixed a non-convergence bug in the inverse Cassini solver (January 2025). That fix addressed the iterative inverse when easting or northing equals the false origin. This is a separate bug in the forward ellipsoidal formulas.

Visual Evidence

Plotting world coastlines transformed via: - PROJ ellipsoidal +proj=cass +ellps=WGS84 (green) - GeographicLib cassini_fwd() (red) - PROJ spherical +proj=cass +R=6378137 (black)

pts <- do.call(cbind, maps::map(plot = F)[1:2])
## these points are wrong
plot(transform_xy(pts, srs_from = "EPSG:4326", srs_to = cass_ellps), 
cex  = 1, pch = 19, col = "green", asp = 1)
#> Warning in .transform_xy(pts_in, srs_from, srs_to): 1972 point(s) had missing
#> values, NA returned in that case

## red + is correct
points(cassini_fwd(pts, lon0 = 0, lat0= 0)[1:2], col = "red", cex  = 2, pch = 19)
points(transform_xy(pts, srs_from = "EPSG:4326", srs_to = cass_spher), pch = 19)
#> Warning in .transform_xy(pts_in, srs_from, srs_to): 1972 point(s) had missing
#> values, NA returned in that case
legend("bottomright", 
       legend = c("PROJ ellipsoidal (wrong)", "GeographicLib (correct)", "PROJ spherical"),
       pch = c(19, 19, 19),
       col = c("green", "red", "black"),
       pt.cex = c(1, 2, 1))

The spherical and GeographicLib results align well. The PROJ ellipsoidal results diverge dramatically at high latitudes and far from the central meridian.