Lightweight CRS parser to remove hard pyproj dependency (#1072)

* Add design spec for fused_overlap and multi_overlap utilities

* Add implementation plan for fused_overlap and multi_overlap

* Add lightweight CRS class with embedded EPSG table (#1057)

* Fix ValueError message to include pyproj install hint (#1057)

* Two-tier CRS resolution: lite CRS first, pyproj fallback (#1057)

* Add scatter-point transform_points() for boundary estimation (#1057)

* Use Numba scatter-point transform for grid boundary estimation (#1057)

* Wire chunk functions and merge helper to use lite CRS (#1057)

* Add integration tests for pyproj-free CRS resolution (#1057)

Author: brendancol

xrspatial/reproject/__init__.py

@@ -192,11 +192,10 @@ def _reproject_chunk_numpy(
     Called inside ``dask.delayed`` for the dask path, or directly for numpy.
     CRS objects are passed as WKT strings for pickle safety.
     """
-    from ._crs_utils import _require_pyproj
+    from ._crs_utils import _crs_from_wkt
 
-    pyproj = _require_pyproj()
-    src_crs = pyproj.CRS.from_wkt(src_wkt)
-    tgt_crs = pyproj.CRS.from_wkt(tgt_wkt)
+    src_crs = _crs_from_wkt(src_wkt)
+    tgt_crs = _crs_from_wkt(tgt_wkt)
 
     # Try Numba fast path first (avoids creating pyproj Transformer)
     numba_result = None
@@ -212,6 +211,8 @@ def _reproject_chunk_numpy(
         src_y, src_x = numba_result
     else:
         # Fallback: create pyproj Transformer (expensive)
+        from ._crs_utils import _require_pyproj
+        pyproj = _require_pyproj()
         transformer = pyproj.Transformer.from_crs(
             tgt_crs, src_crs, always_xy=True
         )
@@ -321,15 +322,10 @@ def _reproject_chunk_cupy(
     """CuPy variant of ``_reproject_chunk_numpy``."""
     import cupy as cp
 
-    from ._crs_utils import _require_pyproj
+    from ._crs_utils import _crs_from_wkt
 
-    pyproj = _require_pyproj()
-    src_crs = pyproj.CRS.from_wkt(src_wkt)
-    tgt_crs = pyproj.CRS.from_wkt(tgt_wkt)
-
-    transformer = pyproj.Transformer.from_crs(
-        tgt_crs, src_crs, always_xy=True
-    )
+    src_crs = _crs_from_wkt(src_wkt)
+    tgt_crs = _crs_from_wkt(tgt_wkt)
 
     # Try CUDA transform first (keeps coordinates on-device)
     cuda_result = None
@@ -371,6 +367,11 @@ def _reproject_chunk_cupy(
         _use_native_cuda = True
     else:
         # CPU fallback (Numba JIT or pyproj)
+        from ._crs_utils import _require_pyproj
+        pyproj = _require_pyproj()
+        transformer = pyproj.Transformer.from_crs(
+            tgt_crs, src_crs, always_xy=True
+        )
         src_y, src_x = _transform_coords(
             transformer, chunk_bounds_tuple, chunk_shape, transform_precision,
             src_crs=src_crs, tgt_crs=tgt_crs,
@@ -513,16 +514,13 @@ def reproject(
         If vertical transformation was applied, ``attrs['vertical_crs']``
         records the target vertical datum.
     """
-    from ._crs_utils import _require_pyproj
-
     if not isinstance(raster, xr.DataArray):
         raise TypeError(
             f"reproject(): raster must be an xr.DataArray, "
             f"got {type(raster).__name__}"
         )
 
     _validate_resampling(resampling)
-    _require_pyproj()
 
     # Resolve CRS
     src_crs = _resolve_crs(source_crs)
@@ -984,11 +982,10 @@ def _reproject_dask_cupy(
     """
     import cupy as cp
 
-    from ._crs_utils import _require_pyproj
+    from ._crs_utils import _crs_from_wkt
 
-    pyproj = _require_pyproj()
-    src_crs = pyproj.CRS.from_wkt(src_wkt)
-    tgt_crs = pyproj.CRS.from_wkt(tgt_wkt)
+    src_crs = _crs_from_wkt(src_wkt)
+    tgt_crs = _crs_from_wkt(tgt_wkt)
 
     # Use larger chunks for GPU to amortize kernel launch overhead
     gpu_chunk = chunk_size or 2048
@@ -1048,6 +1045,8 @@ def _reproject_dask_cupy(
                 c_max = int(np.ceil(c_max_val)) + 3
             else:
                 # CPU fallback for this chunk
+                from ._crs_utils import _require_pyproj
+                pyproj = _require_pyproj()
                 transformer = pyproj.Transformer.from_crs(
                     tgt_crs, src_crs, always_xy=True
                 )
@@ -1120,30 +1119,44 @@ def _reproject_dask_cupy(
 
 
 def _source_footprint_in_target(src_bounds, src_wkt, tgt_wkt):
-    """Compute an approximate bounding box of the source raster in target CRS.
-
-    Transforms corners and edge midpoints (12 points) to handle non-linear
-    projections.  Returns ``(left, bottom, right, top)`` in target CRS, or
-    *None* if the transform fails (e.g. out-of-domain).
-    """
+    """Compute approximate bounding box of source raster in target CRS."""
     try:
-        from ._crs_utils import _require_pyproj
-        pyproj = _require_pyproj()
-        src_crs = pyproj.CRS(src_wkt)
-        tgt_crs = pyproj.CRS(tgt_wkt)
-        transformer = pyproj.Transformer.from_crs(
-            src_crs, tgt_crs, always_xy=True
-        )
+        from ._crs_utils import _crs_from_wkt, _resolve_crs
+        try:
+            src_crs = _crs_from_wkt(src_wkt)
+        except Exception:
+            src_crs = _resolve_crs(src_wkt)
+        try:
+            tgt_crs = _crs_from_wkt(tgt_wkt)
+        except Exception:
+            tgt_crs = _resolve_crs(tgt_wkt)
     except Exception:
         return None
 
     sl, sb, sr, st = src_bounds
     mx = (sl + sr) / 2
     my = (sb + st) / 2
-    xs = [sl, mx, sr, sl, mx, sr, sl, mx, sr, sl, sr, mx]
-    ys = [sb, sb, sb, my, my, my, st, st, st, mx, mx, sb]
+    xs = np.array([sl, mx, sr, sl, mx, sr, sl, mx, sr, sl, sr, mx])
+    ys = np.array([sb, sb, sb, my, my, my, st, st, st, mx, mx, sb])
+
     try:
-        tx, ty = transformer.transform(xs, ys)
+        from ._projections import transform_points
+        result = transform_points(src_crs, tgt_crs, xs, ys)
+        if result is not None:
+            tx, ty = result
+            tx = [v for v in tx if np.isfinite(v)]
+            ty = [v for v in ty if np.isfinite(v)]
+            if not tx or not ty:
+                return None
+            return (min(tx), min(ty), max(tx), max(ty))
+    except (ImportError, ModuleNotFoundError):
+        pass
+
+    try:
+        from ._crs_utils import _require_pyproj
+        pyproj = _require_pyproj()
+        transformer = pyproj.Transformer.from_crs(src_crs, tgt_crs, always_xy=True)
+        tx, ty = transformer.transform(xs.tolist(), ys.tolist())
         tx = [v for v in tx if np.isfinite(v)]
         ty = [v for v in ty if np.isfinite(v)]
         if not tx or not ty:
@@ -1298,14 +1311,11 @@ def merge(
     -------
     xr.DataArray
     """
-    from ._crs_utils import _require_pyproj
-
     if not rasters:
         raise ValueError("merge(): rasters list must not be empty")
 
     _validate_resampling(resampling)
     _validate_strategy(strategy)
-    pyproj = _require_pyproj()
 
     # Resolve target CRS
     tgt_crs = _resolve_crs(target_crs)
@@ -1485,9 +1495,8 @@ def _merge_inmemory(
     Detects same-CRS tiles and uses fast direct placement instead
     of reprojection.
     """
-    from ._crs_utils import _require_pyproj
-    pyproj = _require_pyproj()
-    tgt_crs = pyproj.CRS.from_wkt(tgt_wkt)
+    from ._crs_utils import _crs_from_wkt
+    tgt_crs = _crs_from_wkt(tgt_wkt)
 
     arrays = []
     for info in raster_infos:

xrspatial/reproject/_crs_utils.py

@@ -1,36 +1,86 @@
-"""CRS detection utilities and optional pyproj import guard."""
+"""CRS detection utilities and optional pyproj import guard.
+
+Uses a two-tier strategy: try the lightweight built-in CRS first,
+then fall back to pyproj for codes/formats not in the built-in table.
+"""
 from __future__ import annotations
 
+from xrspatial.reproject._lite_crs import CRS as LiteCRS
 
-def _require_pyproj():
-    """Import and return the pyproj module, raising a clear error if missing."""
+
+def _try_import_pyproj():
+    """Try to import pyproj, returning the module or None."""
     try:
         import pyproj
         return pyproj
     except ImportError:
+        return None
+
+
+def _require_pyproj():
+    """Import and return the pyproj module, raising a clear error if missing."""
+    pyproj = _try_import_pyproj()
+    if pyproj is None:
         raise ImportError(
             "pyproj is required for CRS reprojection. "
             "Install it with:  pip install pyproj  "
             "or:  pip install xarray-spatial[reproject]"
         )
+    return pyproj
 
 
 def _resolve_crs(crs_input):
-    """Convert *crs_input* to a ``pyproj.CRS`` object.
-
-    Accepts anything ``pyproj.CRS()`` accepts: EPSG int, authority string,
-    WKT, proj4 dict, or an existing ``pyproj.CRS`` instance.
-
-    Returns None if *crs_input* is None.
+    """Convert *crs_input* to a CRS object.
+
+    Resolution order:
+
+    1. ``None`` passes through as ``None``.
+    2. An existing ``LiteCRS`` instance passes through unchanged.
+    3. An existing ``pyproj.CRS`` instance passes through unchanged
+       (only checked when pyproj is importable).
+    4. Try ``LiteCRS(crs_input)`` -- covers EPSG ints and ``"EPSG:XXXX"``
+       strings for codes in the built-in table.
+    5. Fall back to ``pyproj.CRS(crs_input)`` -- raises ``ImportError``
+       if pyproj is not installed.
     """
     if crs_input is None:
         return None
-    pyproj = _require_pyproj()
-    if isinstance(crs_input, pyproj.CRS):
+
+    # Pass through existing LiteCRS
+    if isinstance(crs_input, LiteCRS):
+        return crs_input
+
+    # Pass through existing pyproj.CRS (if pyproj available)
+    pyproj = _try_import_pyproj()
+    if pyproj is not None and isinstance(crs_input, pyproj.CRS):
         return crs_input
+
+    # Try lite CRS first
+    try:
+        return LiteCRS(crs_input)
+    except (ValueError, TypeError):
+        pass
+
+    # Fall back to pyproj
+    pyproj = _require_pyproj()
     return pyproj.CRS(crs_input)
 
 
+def _crs_from_wkt(wkt):
+    """Build a CRS from an OGC WKT string.
+
+    Tries ``LiteCRS.from_wkt`` first (extracts the AUTHORITY tag),
+    then falls back to ``pyproj.CRS.from_wkt``.
+    """
+    try:
+        return LiteCRS.from_wkt(wkt)
+    except (ValueError, TypeError):
+        pass
+
+    pyproj = _require_pyproj()
+    return pyproj.CRS.from_wkt(wkt)
+
+
 def _detect_source_crs(raster):
     """Auto-detect the CRS of a DataArray.
 
@@ -47,7 +97,7 @@ def _detect_source_crs(raster):
 
     crs_wkt = raster.attrs.get('crs_wkt')
     if crs_wkt is not None:
-        return _resolve_crs(crs_wkt)
+        return _crs_from_wkt(crs_wkt)
 
     # rioxarray fallback
     try:

xrspatial/reproject/_grid.py

@@ -4,6 +4,38 @@
 import numpy as np
 
 
+def _transform_boundary(source_crs, target_crs, xs, ys):
+    """Transform coordinate arrays, preferring Numba fast path over pyproj.
+
+    Parameters
+    ----------
+    source_crs, target_crs : CRS-like
+        Source and target coordinate reference systems.
+    xs, ys : ndarray
+        1-D arrays of x and y coordinates in *source_crs*.
+
+    Returns
+    -------
+    tx, ty : ndarray
+        Transformed coordinates as numpy arrays.
+    """
+    from ._projections import transform_points
+
+    result = transform_points(source_crs, target_crs, xs, ys)
+    if result is not None:
+        return result
+
+    # Fall back to pyproj
+    from ._crs_utils import _require_pyproj
+
+    pyproj = _require_pyproj()
+    transformer = pyproj.Transformer.from_crs(
+        source_crs, target_crs, always_xy=True
+    )
+    tx, ty = transformer.transform(xs, ys)
+    return np.asarray(tx), np.asarray(ty)
+
+
 def _compute_output_grid(source_bounds, source_shape, source_crs, target_crs,
                          resolution=None, bounds=None, width=None, height=None):
     """Compute the output raster grid parameters.
@@ -14,7 +46,7 @@ def _compute_output_grid(source_bounds, source_shape, source_crs, target_crs,
         (left, bottom, right, top) in source CRS.
     source_shape : tuple
         (height, width) of source raster.
-    source_crs, target_crs : pyproj.CRS
+    source_crs, target_crs : CRS-like
         Source and target coordinate reference systems.
     resolution : float or tuple or None
         Target resolution. If tuple, (x_res, y_res).
@@ -27,13 +59,6 @@ def _compute_output_grid(source_bounds, source_shape, source_crs, target_crs,
     -------
     dict with keys: bounds, shape, res_x, res_y
     """
-    from ._crs_utils import _require_pyproj
-
-    pyproj = _require_pyproj()
-    transformer = pyproj.Transformer.from_crs(
-        source_crs, target_crs, always_xy=True
-    )
-
     if bounds is None:
         # Transform source corners and edges to target CRS
         src_left, src_bottom, src_right, src_top = source_bounds
@@ -76,7 +101,7 @@ def _compute_output_grid(source_bounds, source_shape, source_crs, target_crs,
         ixx, iyy = np.meshgrid(ix, iy)
         xs = np.concatenate([edge_xs, ixx.ravel()])
         ys = np.concatenate([edge_ys, iyy.ravel()])
-        tx, ty = transformer.transform(xs, ys)
+        tx, ty = _transform_boundary(source_crs, target_crs, xs, ys)
         tx = np.asarray(tx)
         ty = np.asarray(ty)
         # Filter out inf/nan from failed transforms
@@ -110,7 +135,9 @@ def _compute_output_grid(source_bounds, source_shape, source_crs, target_crs,
             ix = np.linspace(src_left, src_right, n_dense)
             iy = np.linspace(src_bottom, src_top, n_dense)
             ixx, iyy = np.meshgrid(ix, iy)
-            itx, ity = transformer.transform(ixx.ravel(), iyy.ravel())
+            itx, ity = _transform_boundary(
+                source_crs, target_crs, ixx.ravel(), iyy.ravel()
+            )
             itx = np.asarray(itx)
             ity = np.asarray(ity)
             ivalid = np.isfinite(itx) & np.isfinite(ity)
@@ -150,13 +177,15 @@ def _compute_output_grid(source_bounds, source_shape, source_crs, target_crs,
         src_res_y = (src_top - src_bottom) / src_h
         center_x = (src_left + src_right) / 2
         center_y = (src_bottom + src_top) / 2
-        tc_x, tc_y = transformer.transform(center_x, center_y)
-        # Step along x only
-        tx_x, tx_y = transformer.transform(center_x + src_res_x, center_y)
-        dx = np.hypot(float(tx_x) - float(tc_x), float(tx_y) - float(tc_y))
-        # Step along y only
-        ty_x, ty_y = transformer.transform(center_x, center_y + src_res_y)
-        dy = np.hypot(float(ty_x) - float(tc_x), float(ty_y) - float(tc_y))
+        # Batch the three resolution-estimation points into one call
+        pts_x = np.array([center_x, center_x + src_res_x, center_x])
+        pts_y = np.array([center_y, center_y, center_y + src_res_y])
+        tp_x, tp_y = _transform_boundary(source_crs, target_crs, pts_x, pts_y)
+        tc_x, tc_y = float(tp_x[0]), float(tp_y[0])
+        tx_x, tx_y = float(tp_x[1]), float(tp_y[1])
+        ty_x, ty_y = float(tp_x[2]), float(tp_y[2])
+        dx = np.hypot(tx_x - tc_x, tx_y - tc_y)
+        dy = np.hypot(ty_x - tc_x, ty_y - tc_y)
         if dx == 0 or dy == 0:
             res_x = (right - left) / src_w
             res_y = (top - bottom) / src_h