Align code / tests with the most recent PROJ

.gitignore

@@ -36,3 +36,5 @@ _SUCCESS
 tags
 .vscode
 .metals
+
+.claude

CHANGELOG.md

@@ -7,8 +7,43 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 ## [Unreleased]
 
+### Changed
+- `+f` is the flattening and `+rf` the reciprocal flattening, as in PROJ. The two were swapped, so definitions that passed a `1/f` value under `+f` (e.g. `+f=298.257222101`) built a nonsensical ellipsoid; they must now use `+rf`. PROJ rejects such a `+f` value outright
+- Transverse Mercator (`+proj=tmerc`) uses the exact Poder/Engsager algorithm for ellipsoids, which has been PROJ's default since 6.0, and accepts `+approx` to select the Evenden/Snyder series. Inside a normal zone the two agree to well under a millimetre; 40 degrees from the central meridian they differ by ~1.5 km, and beyond ~80 degrees the series returns garbage where the exact algorithm reports the point as outside the projection domain
+- Equidistant Cylindrical (`+proj=eqc`) implements the ellipsoidal method (EPSG:1028) as well as the spherical one (EPSG:1029), matching PROJ 9.8.0. Northings on an ellipsoid are now meridional arc lengths rather than `a * phi`
+
 ### Fixed
 - Oblique Mercator: compute the natural-origin (uc) offset from the central-line azimuth (`+alpha`) instead of the rectified bearing (`+gamma`), so the projection centre maps to the false easting/northing when `+gamma` differs from `+alpha` (e.g. an explicit `+gamma=0` with a non-zero azimuth)
+- Oblique Mercator: `+alpha` alone no longer behaves as `+gamma=0`, and `+gamma` alone no longer falls back to an azimuth of -45 degrees. Both parameters are now tracked as given, as PROJ does, and the two-point form (`+lat_1`/`+lon_1`/`+lat_2`/`+lon_2`) is reachable
+- Oblique Mercator: use `atan2` in the forward transform, as PROJ does; the previous `atan(y/x)` plus a fixed `+PI` correction picked the wrong branch for points more than 90 degrees from the central line, and the near-meridian fallback had a spurious `B` factor
+- `+R` defines a sphere. It used to set only the semi-major axis and leave the ellipsoid's eccentricity in place, which shifted every projection with an ellipsoidal branch (`merc`, `tmerc`, `laea`, `aea`, `cea`, `leac`, `lcc`, `stere`, `sterea`, `somerc`, `aeqd`, `omerc`, `geos`, `cass`) by tens to hundreds of kilometres
+- Mercator: `+lat_ts` was ignored; it now sets the scale factor as in PROJ (`cos(lat_ts)` on a sphere, `msfn(lat_ts)` on an ellipsoid)
+- Equal Area Cylindrical: a supplied `+k_0` was overwritten by `cos(lat_ts)` even when `+lat_ts` was absent
+- The Snyder conics (`euler`, `murd1`, `murd2`, `murd3`, `pconic`, `vitk1`) ignored `+lat_1`/`+lat_2` and used hardcoded 30/60 degree parallels; their shared inverse also used the unshifted northing in `atan2` and dropped the sign flip for a negative cone constant
+- Equidistant Conic (`+proj=eqdc`) was not PROJ's algorithm at all: it used a Lambert-Conformal-style formulation with a hardcoded eccentricity of 0.822719, a unit radius, and hardcoded standard parallels, and its inverse was never wired into the framework. Rewritten from PROJ's `eqdc.cpp`
+- Bonne: `+lat_1` was ignored and the standard parallel was pinned at 90 degrees (the Werner limit)
+- Polyconic: the ellipsoidal branch was unreachable (`spherical` was forced true), its forward transform used an uninitialised value in place of the longitude, and its inverse used `1/es` where PROJ uses `1 - es`
+- Cassini: sign error in the ellipsoidal forward easting series (`C1 - ...` instead of `C1 + ...`)
+- Rectangular Polyconic: `+lat_ts` and `+lat_0` were ignored
+- Transverse Mercator: the spherical branch applied the scale factor twice to the easting
+- Winkel Tripel: `+lat_1` was ignored, and the Winkel constant was being passed as the latitude of origin
+- Krovak: PROJ's fixed defaults (Bessel 1841, `lat_0` 49°30′N, `lon_0` 42°30′ of Ferro − 17°40′, `k_0` 0.9999) are applied when the definition omits them, and `+czech` selects the native westing/southing orientation
+- Bipolar Conic: dropped a spurious `lon_0` default of −90 degrees that PROJ does not have
+- Near-sided Perspective (`+proj=nsper`): `+h` threw `NoSuchElementException`, the aspect (oblique/equatorial/polar) was hardcoded to equatorial, and the far-side-of-the-globe domain check was commented out
+- Space Oblique Mercator for Landsat (`+proj=lsat`): `+lsat` and `+path` were rejected as unsupported and the satellite/path were hardcoded to 1/120
+- Hammer: `+W` and `+M` were rejected as unsupported, and the initialisation overwrote whatever was set with the defaults
+- Lagrange: `+W` was rejected as unsupported, and the default was 1.4 instead of PROJ's 2
+- Urmaev Flat-Polar Sinusoidal (`+proj=urmfps`): `+n` was rejected as unsupported
+- Lambert Equal Area Conic (`+proj=leac`): `+south` threw `NoSuchElementException`
+- `putp2`, `nell`, `mbtfpq`, `mbt_fps`: the Newton iteration ran on the output northing instead of the latitude, so it never converged and the unconverged latitude was used
+- `wag1`/`urmfps`, `wag2`, `mbtfpp`: the transformed latitude was computed and then discarded, and the raw latitude used in its place
+- Loximuthal: a missing `abs()` in the near-parallel test sent every point south of `+lat_1` down the wrong branch
+- Robinson: latitudes landing exactly on a 5-degree table node selected the node below through `floor()` rounding
+- `+units=ch`, `+units=fath`, `+units=link` and `+units=us-ch` were silently treated as metres because those units, though defined, were missing from the lookup table; the Indian units (`ind-yd`, `ind-ft`, `ind-ch`) were added
+- `+vunits` and `+vto_meter` are accepted instead of raising `UnsupportedParameterException`. The EPSG resource ships 158 compound (horizontal + vertical) definitions carrying `+vunits`, so 144 of them - `OSGB36 / British National Grid + ODN height`, `Amersfoort / RD New + NAP height`, the `ETRS89 / NTM zone N`, `ETRS89 / UTM zone N` and `NAD83(CSRS) / UTM zone N` families, `SWEREF99 xx + RH2000`, `Tokyo + JSLD`, and others - could not be constructed at all. The vertical unit is recorded but **not** applied: proj4j carries the vertical ordinate in metres, as it did before. The horizontal transform is unaffected either way, which is all PROJ does with it for the 2D part
+
+### Added
+- `ProjAlignmentTest`, pinning 67 cases against values generated from raw PROJ pipelines and cross-checked on PROJ 9.5.1 and 9.6.0
 
 ## [1.4.3] - 2026-06-02
 

core/src/main/java/org/locationtech/proj4j/parser/DatumParameters.java

@@ -116,20 +116,38 @@ public void setES(double es) {
         this.es = es;
     }
 
+    /**
+     * Sets the reciprocal (inverse) flattening, {@code +rf}.
+     * PROJ: {@code f = 1/rf; es = 2f - f*f}.
+     */
     public void setRF(double rf) {
         ellipsoid = null;  // force user-defined ellipsoid
-        es = rf * (2. - rf);
+        double f = 1. / rf;
+        es = f * (2. - f);
     }
 
     public void setR_A() {
         ellipsoid = null;  // force user-defined ellipsoid
         a *= 1. - es * (SIXTH + es * (RA4 + es * RA6));
     }
 
+    /**
+     * Sets the flattening, {@code +f}.
+     * PROJ: {@code es = 2f - f*f}.
+     */
     public void setF(double f) {
         ellipsoid = null;  // force user-defined ellipsoid
-        double rf = 1.0 / f;
-        es = rf * (2. - rf);
+        es = f * (2. - f);
+    }
+
+    /**
+     * Sets the radius of the sphere, {@code +R}.
+     * PROJ: "specifying R overrules everything" - the ellipsoid becomes a sphere.
+     */
+    public void setR(double r) {
+        ellipsoid = null;  // force user-defined ellipsoid
+        a = r;
+        es = 0;
     }
 
     public double getA() {

core/src/main/java/org/locationtech/proj4j/parser/Proj4Keyword.java

@@ -38,6 +38,8 @@ public class Proj4Keyword {
     public static final String lat_0 = "lat_0";
     public static final String lat_1 = "lat_1";
     public static final String lat_2 = "lat_2";
+    public static final String lon_1 = "lon_1";
+    public static final String lon_2 = "lon_2";
     public static final String lon_0 = "lon_0";
     public static final String lonc = "lonc";
     public static final String pm = "pm";
@@ -57,6 +59,8 @@ public class Proj4Keyword {
 
     public static final String south = "south";
     public static final String to_meter = "to_meter";
+    public static final String vunits = "vunits";
+    public static final String vto_meter = "vto_meter";
     public static final String towgs84 = "towgs84";
     public static final String units = "units";
     public static final String x_0 = "x_0";
@@ -68,6 +72,13 @@ public class Proj4Keyword {
     public static final String no_defs = "no_defs";
     public static final String wktext = "wktext";
     public static final String no_uoff = "no_uoff";
+    public static final String czech = "czech";
+    public static final String approx = "approx";
+    public static final String W = "W";
+    public static final String M = "M";
+    public static final String n = "n";
+    public static final String lsat = "lsat";
+    public static final String path = "path";
 
 
     private static Set<String> supportedParams = null;
@@ -95,6 +106,8 @@ public static synchronized Set supportedParameters() {
             supportedParams.add(lat_0);
             supportedParams.add(lat_1);
             supportedParams.add(lat_2);
+            supportedParams.add(lon_1);
+            supportedParams.add(lon_2);
             supportedParams.add(lon_0);
             supportedParams.add(lonc);
 
@@ -105,13 +118,22 @@ public static synchronized Set supportedParameters() {
             supportedParams.add(south);
             supportedParams.add(towgs84);
             supportedParams.add(to_meter);
+            supportedParams.add(vunits);      // vertical unit of a compound CRS; parsed, not applied
+            supportedParams.add(vto_meter);   // vertical unit of a compound CRS; parsed, not applied
             supportedParams.add(units);
             supportedParams.add(nadgrids);
             supportedParams.add(pm);
             supportedParams.add(axis);
 
             supportedParams.add(gamma);       // Just for Oblique Mercator projection
             supportedParams.add(no_uoff);     // Just for Oblique Mercator projection
+            supportedParams.add(czech);       // Just for the Krovak projection
+            supportedParams.add(approx);      // Transverse Mercator algorithm selection
+            supportedParams.add(W);           // Hammer, Lagrange
+            supportedParams.add(M);           // Hammer
+            supportedParams.add(n);           // Urmaev Flat-Polar Sinusoidal
+            supportedParams.add(lsat);        // Space Oblique Mercator for Landsat
+            supportedParams.add(path);        // Space Oblique Mercator for Landsat
             supportedParams.add(zone);        // Just for Transverse Mercator projection
 
             supportedParams.add(title);       // no-op

core/src/main/java/org/locationtech/proj4j/parser/Proj4Parser.java

@@ -24,6 +24,12 @@
 import org.locationtech.proj4j.datum.Ellipsoid;
 import org.locationtech.proj4j.datum.Grid;
 import org.locationtech.proj4j.proj.ExtendedTransverseMercatorProjection;
+import org.locationtech.proj4j.proj.HammerProjection;
+import org.locationtech.proj4j.proj.KrovakProjection;
+import org.locationtech.proj4j.proj.LandsatProjection;
+import org.locationtech.proj4j.proj.ObliqueMercatorProjection;
+import org.locationtech.proj4j.proj.LagrangeProjection;
+import org.locationtech.proj4j.proj.UrmaevFlatPolarSinusoidalProjection;
 import org.locationtech.proj4j.proj.Projection;
 import org.locationtech.proj4j.proj.TransverseMercatorProjection;
 import org.locationtech.proj4j.units.Angle;
@@ -142,6 +148,22 @@ private Projection parseProjection(Map params, Ellipsoid ellipsoid) {
         s = (String) params.get(Proj4Keyword.to_meter);
         if (s != null)
             projection.setFromMetres(1.0 / Double.parseDouble(s));
+
+        /*
+         * Vertical unit of a compound (horizontal + vertical) CRS. Recorded so those definitions
+         * parse - the EPSG resource ships 158 of them - but not applied: proj4j keeps the vertical
+         * ordinate in metres. The horizontal transform PROJ performs is the same either way.
+         */
+        s = (String) params.get(Proj4Keyword.vunits);
+        if (s != null) {
+            Unit vunit = Units.findUnits(s);
+            if (vunit != null)
+                projection.setVerticalUnit(vunit);
+        }
+
+        s = (String) params.get(Proj4Keyword.vto_meter);
+        if (s != null)
+            projection.setVerticalFromMetres(1.0 / Double.parseDouble(s));
 
         s = (String) params.get(Proj4Keyword.h);
         if (s != null) {
@@ -168,11 +190,51 @@ private Projection parseProjection(Map params, Ellipsoid ellipsoid) {
 
         // this must be done last, since behaviour depends on other params being set (eg +south)
         if (projection instanceof TransverseMercatorProjection) {
+            if (params.containsKey(Proj4Keyword.approx))
+                ((TransverseMercatorProjection) projection).setApprox(true);
             s = (String) params.get(Proj4Keyword.zone);
             if (s != null)
                 ((TransverseMercatorProjection) projection).setUTMZone(Integer
                         .parseInt(s));
         }
+        if (projection instanceof HammerProjection) {
+            s = (String) params.get(Proj4Keyword.W);
+            if (s != null)
+                ((HammerProjection) projection).setW(Double.parseDouble(s));
+            s = (String) params.get(Proj4Keyword.M);
+            if (s != null)
+                ((HammerProjection) projection).setM(Double.parseDouble(s));
+        }
+        if (projection instanceof LagrangeProjection) {
+            s = (String) params.get(Proj4Keyword.W);
+            if (s != null)
+                ((LagrangeProjection) projection).setW(Double.parseDouble(s));
+        }
+        if (projection instanceof UrmaevFlatPolarSinusoidalProjection) {
+            s = (String) params.get(Proj4Keyword.n);
+            if (s != null)
+                ((UrmaevFlatPolarSinusoidalProjection) projection).setN(Double.parseDouble(s));
+        }
+        if (projection instanceof ObliqueMercatorProjection) {
+            s = (String) params.get(Proj4Keyword.lon_1);
+            if (s != null)
+                ((ObliqueMercatorProjection) projection).setProjectionLongitude1Degrees(parseAngle(s));
+            s = (String) params.get(Proj4Keyword.lon_2);
+            if (s != null)
+                ((ObliqueMercatorProjection) projection).setProjectionLongitude2Degrees(parseAngle(s));
+        }
+        if (projection instanceof LandsatProjection) {
+            s = (String) params.get(Proj4Keyword.lsat);
+            if (s != null)
+                ((LandsatProjection) projection).setLandsat(Integer.parseInt(s));
+            s = (String) params.get(Proj4Keyword.path);
+            if (s != null)
+                ((LandsatProjection) projection).setPath(Integer.parseInt(s));
+        }
+        if (projection instanceof KrovakProjection) {
+            if (params.containsKey(Proj4Keyword.czech))
+                ((KrovakProjection) projection).setCzech(true);
+        }
         if (projection instanceof ExtendedTransverseMercatorProjection) {
             s = (String) params.get(Proj4Keyword.zone);
             if (s != null)
@@ -245,9 +307,15 @@ private void parseEllipsoid(Map params, DatumParameters datumParam) {
         String s;
 
    /*
-    * // not supported by PROJ4 s = (String) params.get(Proj4Param.R); if (s !=
-    * null) a = Double.parseDouble(s);
+    * Radius of the sphere, given in meters. As in PROJ, "specifying R overrules
+    * everything": the ellipsoid is a sphere of that radius, and +ellps, +datum
+    * and the shape parameters are ignored.
     */
+        s = (String) params.get(Proj4Keyword.R);
+        if (s != null) {
+            datumParam.setR(Double.parseDouble(s));
+            return;
+        }
 
         String code = (String) params.get(Proj4Keyword.ellps);
         if (code != null) {

core/src/main/java/org/locationtech/proj4j/proj/AitoffProjection.java

@@ -22,6 +22,7 @@
 import java.util.Objects;
 
 import org.locationtech.proj4j.ProjCoordinate;
+import org.locationtech.proj4j.ProjectionException;
 
 
 public class AitoffProjection extends PseudoCylindricalProjection {
@@ -33,11 +34,18 @@ public class AitoffProjection extends PseudoCylindricalProjection {
 	private double cosphi1 = 0;
 
 	public AitoffProjection() {
+		// NaN marks "+lat_1 not given", so the Winkel Tripel default can be applied
+		projectionLatitude1 = Double.NaN;
+	}
+
+	protected AitoffProjection(int type) {
+		this();
+		winkel = type == WINKEL;
 	}
 
 	public AitoffProjection(int type, double projectionLatitude) {
+		this(type);
 		this.projectionLatitude = projectionLatitude;
-		winkel = type == WINKEL;
 	}
 
 	public ProjCoordinate project(double lplam, double lpphi, ProjCoordinate out) {
@@ -59,11 +67,10 @@ public ProjCoordinate project(double lplam, double lpphi, ProjCoordinate out) {
 	public void initialize() {
 		super.initialize();
 		if (winkel) {
-//FIXME
-//			if (pj_param(P->params, "tlat_1").i)
-//				if ((cosphi1 = Math.cos(pj_param(P->params, "rlat_1").f)) == 0.)
-//					throw new IllegalArgumentException("-22")
-//			else /* 50d28' or acos(2/pi) */
+			if (!Double.isNaN(projectionLatitude1)) {
+				if ((cosphi1 = Math.cos(projectionLatitude1)) == 0.)
+					throw new ProjectionException("Invalid value for lat_1: |lat_1| should be < 90");
+			} else /* 50d28' or acos(2/pi) */
 				cosphi1 = 0.636619772367581343;
 		}
 	}

core/src/main/java/org/locationtech/proj4j/proj/BipolarProjection.java

@@ -50,7 +50,8 @@ public class BipolarProjection extends Projection {
 	public BipolarProjection() {
 		minLatitude = Math.toRadians(-80);
 		maxLatitude = Math.toRadians(80);
-		projectionLongitude = Math.toRadians(-90);
+		// no lon_0 default: the projection constants already carry the bipolar
+		// orientation, and PROJ leaves lon_0 at 0 unless it is given
 		minLongitude = Math.toRadians(-90);
 		maxLongitude = Math.toRadians(90);
 	}

core/src/main/java/org/locationtech/proj4j/proj/BonneProjection.java

@@ -94,10 +94,9 @@ public void initialize() {
 
 		double c;
 
-//		phi1 = pj_param(params, "rlat_1").f;
-		phi1 = ProjectionMath.HALFPI;
+		phi1 = projectionLatitude1;
 		if (Math.abs(phi1) < EPS10)
-			throw new ProjectionException("-23");
+			throw new ProjectionException("Invalid value for lat_1: |lat_1| should be > 0");
 		if (!spherical) {
 			en = ProjectionMath.enfn(es);
 			m1 = ProjectionMath.mlfn(phi1, am1 = Math.sin(phi1),

core/src/main/java/org/locationtech/proj4j/proj/CassiniProjection.java

@@ -63,7 +63,7 @@ public ProjCoordinate project(double lplam, double lpphi, ProjCoordinate xy) {
 			c *= es * c / (1 - es);
 			a2 = a1 * a1;
 			xy.x = n * a1 * (1. - a2 * t *
-				(C1 - (8. - t + 8. * c) * a2 * C2));
+				(C1 + (8. - t + 8. * c) * a2 * C2));
 			xy.y -= m0 - n * tn * a2 *
 				(.5 + (5. - t + 6. * c) * a2 * C3);
 		}

core/src/main/java/org/locationtech/proj4j/proj/CylindricalEqualAreaProjection.java

@@ -29,7 +29,8 @@ public class CylindricalEqualAreaProjection extends Projection {
 	private double[] apa;
 
 	public CylindricalEqualAreaProjection() {
-		this(0.0, 0.0, 0.0);
+		// NaN marks "+lat_ts not given", so that a supplied +k_0 survives
+		this(0.0, 0.0, Double.NaN);
 	}
 
 	public CylindricalEqualAreaProjection(double projectionLatitude, double projectionLongitude, double trueScaleLatitude) {
@@ -41,9 +42,15 @@ public CylindricalEqualAreaProjection(double projectionLatitude, double projecti
 
 	public void initialize() {
 		super.initialize();
-		double t = trueScaleLatitude;
+		double t = 0.0;
 
-		scaleFactor = Math.cos(t);
+		// as in PROJ, lat_ts defines the scale factor; without it +k_0 is used as given
+		if (!Double.isNaN(trueScaleLatitude)) {
+			t = trueScaleLatitude;
+			scaleFactor = Math.cos(t);
+			if (scaleFactor < 0.)
+				throw new ProjectionException("Invalid value for lat_ts: |lat_ts| should be <= 90");
+		}
 		if (es != 0) {
 			t = Math.sin(t);
 			scaleFactor /= Math.sqrt(1. - es * t * t);