reproject: LAEA inverse off by up to 2.6 km, AEA/CEA inverse off by ~5 m vs PROJ

[brendancol]

Describe the bug

Two numerical bugs in the inverse projection kernels in xrspatial/reproject/_projections.py (and their CUDA twins in _projections_cuda.py) make reprojections into LAEA, AEA, and CEA targets sample the wrong source location.

1. _laea_inv_point (oblique/equatorial modes) normalizes the input by x / (a * xmf) and y / (a * ymf), which strips the rq factor that PROJ keeps in rho. The angular distance is then computed as sce = 2 * asin(0.5 * rho / rq), dividing by rq a second time. The result is a radial error of about 0.11% of the distance from the projection origin. For EPSG:3035 over Europe that is up to 2.6 km versus pyproj. The CUDA kernel _d_laea_inv has the same formula.

2. _authalic_apa claims sub-centimetre accuracy from a 6-term inverse authalic latitude series, but the coefficients do not invert _authalic_q. Measured against exact numerical inversion of q(phi), the series has a max error of 4.79 m (7.5e-7 rad). PROJ's standard 3-term series (517/5040, 23/360 + 251/3780, 761/45360) lands at 1.6 mm. The wrong coefficient is visible at leading order in term 2: e417/360 vs PROJ's e423/360. This series feeds the AEA and CEA inverses (_aea_inv_point, _cea_inv_point, plus the inline 3-term copy in _laea_inv_point and the CUDA _d_authalic_inv call sites).

Measured parity vs pyproj (random points in each CRS's region of use):

pair	direction	max error
EPSG:3035 -> 4326	laea inverse	2643 m
EPSG:5070 -> 4326	aea inverse	4.86 m
EPSG:6933 -> 4326	cea inverse	4.81 m
4326 -> each of the above	forward	<= 1.1e-5 m

The inverse direction is what reproject() uses to build the sampling grid when the target CRS is LAEA/AEA/CEA, so reproject(raster_4326, 'EPSG:3035') georeferences the output wrong by the amounts above. With the LAEA normalization fixed and PROJ's series coefficients in place, the 3035 inverse parity drops to ~1 mm.

Expected behavior

Inverse fast-path kernels should match pyproj within the same tolerance the forward kernels already achieve (sub-mm to ~1e-5 m).

Reproduction

import numpy as np, pyproj
from xrspatial.reproject._projections import transform_points

tr = pyproj.Transformer.from_crs(4326, 3035, always_xy=True)
lons = np.random.uniform(-10, 30, 500); lats = np.random.uniform(35, 70, 500)
px, py = tr.transform(lons, lats)
tx, ty = transform_points(pyproj.CRS(3035), pyproj.CRS(4326), np.asarray(px), np.asarray(py))
inv = pyproj.Transformer.from_crs(3035, 4326, always_xy=True)
rx, ry = inv.transform(px, py)
err = np.hypot((tx - rx) * np.cos(np.radians(ry)), ty - ry) * 111320
print(err.max())   # ~2600 m

Additional context

Found by /sweep-accuracy against the reproject module. End-to-end check: reproject(..., 'EPSG:3035') with the default fast path vs transform_precision=0 (exact pyproj) differs accordingly; numpy and cupy backends agree with each other because both kernels carry the same bugs.