Fix various bugs and docs related to MODIS interpolation

AGENTS.md

@@ -106,7 +106,7 @@ graph without computing.
   `# cython: language_level=3, boundscheck=False, cdivision=True, wraparound=False, initializedcheck=False, nonecheck=False`
 - **`wraparound=False` means negative indexing is banned** and silently wrong if used. There are past
   bugfix commits for exactly this. The sole exception is `_get_coords_5km`, explicitly annotated
-  `@cython.wraparound(True)` (`_modis_interpolator.pyx:231`) so it can use `x[-2]`, `x[-7]`.
+  `@cython.wraparound(True)` (`_modis_interpolator.pyx:237`) so it can use `x[-2]`, `x[-7]`.
 - **Two different `floating` fused types.** `_modis_utils.pxd` defines its own (float32/float64) and the
   MODIS `.pyx` files `cimport` it; `multilinear_cython.pyx` uses `from cython cimport floating`, the
   builtin, which *also* includes `long double`.
@@ -119,6 +119,9 @@ graph without computing.
   and a `Free Threading :: 1 - Unstable` classifier.
 - `_simple_modis_interpolator.pyx` redundantly re-applies the file-level directives as per-function
   decorators; `_modis_interpolator.pyx` does not.
+- **Compile-time constants use `cdef extern from *` with a C `#define`**, not `DEF` (deprecated in
+  Cython 3): `EARTH_RADIUS` in `_modis_utils.pyx` and `R`/`H` in `_modis_interpolator.pyx`. This keeps
+  them true compile-time constants inside the `nogil` hot loops rather than module-level global reads.
 
 ## Build, test, lint
 
@@ -151,24 +154,15 @@ make -C doc doctest
 
 Verified as of this writing; fix them only when the task calls for it.
 
-- `_modis_utils.pyx:103` — `if first_arr.ndim != 2 or first_arr.ndim != 2:` tests the same condition
-  twice; the second was presumably meant to validate another array.
-- `_modis_utils.pyx:173` — `good_col_chunks = len(col_chunks) == 1 and col_chunks[0] != num_cols` is
-  always `False` (a single column chunk must equal `num_cols`), so the rechunk branch always runs.
-  The `!=` looks like it should be `==`.
-- `_modis_utils.pyx:178` — the error message hardcodes "(10 rows per scan)" regardless of resolution.
-- `_modis_interpolator.pyx:4` imports `scanline_mapblocks` from `.simple_modis_interpolator` rather than
-  from `._modis_utils` where it is defined, coupling the two MODIS front-ends for no reason.
 - `multilinear_cython.pyx` uses `prange`, but `setup.py` compiles with only `-O3` and no
   `-fopenmp`/`/openmp`, so those loops are serial. `multilinear_interpolation_5d` is unreachable — the
   dispatcher raises for `d > 4`.
-- **Three Earth radii**: `6370997.0` (`__init__.py`, `geointerpolator.py`, `_modis_utils.pyx`),
-  `6370.997` km (`_modis_interpolator.pyx`), `6371008.7714` (`viiinterpolator.py`).
+- **Three Earth radii, accepted as-is**: `6370997.0` (`__init__.py`, `geointerpolator.py`,
+  `_modis_utils.pyx`) and `6370.997` km (`_modis_interpolator.pyx`) are the same value in different
+  units; `6371008.7714` (`viiinterpolator.py`) is the IUGG mean radius, a genuinely different number,
+  so unifying it would change `viiinterpolator` results. Deliberate, not an oversight.
 - `AbstractMultipleInterpolator.interpolate` (inherited by both `Multiple*Interpolator` classes)
   returns a **generator**, not a tuple.
-- `simple_modis_interpolator.interpolate_geolocation_cartesian`'s docstring documents a `res_factor`
-  argument that no longer exists.
-- `DEF` compile-time constants (`DEF R`, `DEF H`, `DEF EARTH_RADIUS`) are deprecated in Cython 3.
 - `interpolator.py`'s `Interpolator` docstring describes `kx_`/`ky_` as orders "in x and y", but
   `_interp` passes `kx=self.kx_` to the *row* (first) axis of `RectBivariateSpline`.
 
@@ -183,7 +177,7 @@ start it unprompted.
   constants; reconcile module naming (`modisinterpolator.py`/`viiinterpolator.py` vs
   `simple_modis_interpolator.py`/`basic_interpolator.py`); document or rename the undocumented `h`
   prefix (`hrow_indices`/`hcol_indices` = high-resolution) and the `kx_`/`x__` trailing/double
-  underscore habits; fix the `course_col_idx` typo (should be `coarse_`); remove dead code.
+  underscore habits; remove dead code.
 - **Interpolator performance.** The MODIS Cython kernels and the `scanline_mapblocks` chunking are the
   hot spots. Any change must hold the existing geodetic-distance tolerances and must not introduce dask
   computes (`CustomScheduler(0)` will fail the tests if it does).

geotiepoints/_modis_interpolator.pyx

@@ -1,15 +1,21 @@
 # cython: language_level=3, boundscheck=False, cdivision=True, wraparound=False, initializedcheck=False, nonecheck=False
 cimport cython
 from ._modis_utils cimport lonlat2xyz, xyz2lonlat, floating, deg2rad
-from .simple_modis_interpolator import scanline_mapblocks
+from ._modis_utils import scanline_mapblocks
 
 from libc.math cimport asin, sin, cos, sqrt
 cimport numpy as np
 import numpy as np
 
-DEF R = 6370.997
-# Aqua altitude in km
-DEF H = 709.0
+cdef extern from *:
+    """
+    #define R 6370.997
+    /* Aqua altitude in km */
+    #define H 709.0
+    """
+    const double R
+    # Aqua altitude in km
+    const double H
 
 np.import_array()
 
@@ -563,14 +569,14 @@ cdef class MODISInterpolator:
             self,
             floating[:, :] input_arr,
             floating[:, ::1] expanded_arr,
-            Py_ssize_t course_col_idx,
+            Py_ssize_t coarse_col_idx,
             Py_ssize_t fine_col_idx,
     ) noexcept nogil:
         cdef floating tiepoint_value
         cdef Py_ssize_t row_idx, row_offset
         for row_idx in range(input_arr.shape[0]):
             row_offset = row_idx * self._fine_pixels_per_coarse_pixel * 2
-            tiepoint_value = input_arr[row_idx, course_col_idx]
+            tiepoint_value = input_arr[row_idx, coarse_col_idx]
             self._expand_tiepoint_array_5km_with_repeat(
                 tiepoint_value,
                 expanded_arr,

geotiepoints/_modis_utils.pyx

@@ -20,7 +20,11 @@ except ImportError:
     xr = None
 
 
-DEF EARTH_RADIUS = 6370997.0
+cdef extern from *:
+    """
+    #define EARTH_RADIUS 6370997.0
+    """
+    const double EARTH_RADIUS
 
 
 cdef void lonlat2xyz(
@@ -100,7 +104,7 @@ def scanline_mapblocks(func):
         if coarse_resolution is None or fine_resolution is None:
             raise ValueError("'coarse_resolution' and 'fine_resolution' are required keyword arguments.")
         first_arr = [arr for arr in args if hasattr(arr, "ndim")][0]
-        if first_arr.ndim != 2 or first_arr.ndim != 2:
+        if any(arr.ndim != 2 for arr in args if hasattr(arr, "ndim")):
             raise ValueError("Expected 2D input arrays.")
         if hasattr(first_arr, "compute"):
             # assume it is dask or xarray with dask, ensure proper chunk size
@@ -170,20 +174,20 @@ def _rechunk_dask_arrays_if_needed(args, rows_per_scan: int):
     num_rows = first_arr.shape[0]
     num_cols = first_arr.shape[1]
     good_row_chunks = all(x % rows_per_scan == 0 for x in row_chunks)
-    good_col_chunks = len(col_chunks) == 1 and col_chunks[0] != num_cols
+    good_col_chunks = len(col_chunks) == 1 and col_chunks[0] == num_cols
     all_orig_chunks = [arr.chunks for arr in args if hasattr(arr, "chunks")]
 
     if num_rows % rows_per_scan != 0:
         raise ValueError("Input longitude/latitude data does not consist of "
-                         "whole scans (10 rows per scan).")
+                         f"whole scans ({rows_per_scan} rows per scan).")
     all_same_chunks = all(
         all_orig_chunks[0] == some_chunks
         for some_chunks in all_orig_chunks[1:]
     )
     if good_row_chunks and good_col_chunks and all_same_chunks:
         return args
 
-    new_row_chunks = (row_chunks[0] // rows_per_scan) * rows_per_scan
+    new_row_chunks = max(1, row_chunks[0] // rows_per_scan) * rows_per_scan
     new_args = [arr.rechunk((new_row_chunks, -1)) if hasattr(arr, "chunks") else arr for arr in args]
     return new_args
 

geotiepoints/simple_modis_interpolator.py

@@ -30,8 +30,10 @@ def interpolate_geolocation_cartesian(lon_array, lat_array, coarse_resolution, f
             The input data is expected to represent 1000m geolocation.
         lat_array: Latitude data as a 2D numpy, dask, or xarray DataArray object.
             The input data is expected to represent 1000m geolocation.
-        res_factor (int): Expansion factor for the function. Should be 2 for
-            500m output or 4 for 250m output.
+        coarse_resolution (int): Resolution in meters of the input arrays.
+            Keyword-only; consumed by the ``scanline_mapblocks`` decorator.
+        fine_resolution (int): Resolution in meters of the output arrays.
+            Keyword-only; consumed by the ``scanline_mapblocks`` decorator.
 
     Returns:
         A two-element tuple (lon, lat).

geotiepoints/tests/test_simple_modis_interpolator.py

@@ -5,7 +5,12 @@
 import dask
 import dask.array as da
 
-from geotiepoints.simple_modis_interpolator import modis_1km_to_250m, modis_1km_to_500m
+from geotiepoints._modis_utils import _rechunk_dask_arrays_if_needed
+from geotiepoints.simple_modis_interpolator import (
+    interpolate_geolocation_cartesian,
+    modis_1km_to_250m,
+    modis_1km_to_500m,
+)
 from .test_modisinterpolator import (
     assert_geodetic_distance,
     load_1km_lonlat_as_xarray_dask,
@@ -51,3 +56,44 @@ def test_nonstandard_scan_size():
     lat1 = lat1[:-1]
 
     pytest.raises(ValueError, modis_1km_to_250m, lon1, lat1)
+
+
+def test_3d_array_raises():
+    """Every array argument must be 2D, not just the first one."""
+    lon1, lat1 = load_1km_lonlat_as_numpy()
+    lat1 = lat1[np.newaxis]
+
+    with pytest.raises(ValueError, match="Expected 2D input arrays"):
+        modis_1km_to_250m(lon1, lat1)
+
+
+def test_missing_resolutions_raises():
+    """The resolutions are required keyword arguments of the decorator."""
+    lon1, lat1 = load_1km_lonlat_as_numpy()
+
+    with pytest.raises(ValueError, match="required keyword arguments"):
+        interpolate_geolocation_cartesian(lon1, lat1)
+
+
+def test_aligned_chunks_are_not_rechunked():
+    """Scan-aligned, full-width, identically-chunked arrays are passed through."""
+    lon1, lat1 = load_1km_lonlat_as_dask()
+    assert lon1.chunks == ((20,), (1354,))
+
+    with dask.config.set(scheduler=CustomScheduler(0)):
+        result = _rechunk_dask_arrays_if_needed([lon1, lat1], 10)
+
+    assert result[0] is lon1
+    assert result[1] is lat1
+
+
+def test_nonstandard_scan_size_error_names_rows_per_scan():
+    """The whole scans error message uses the actual rows per scan."""
+    lon1, lat1 = load_1km_lonlat_as_xarray_dask()
+    # remove 1 row from the end so 5km's 2 rows per scan doesn't divide evenly
+    lon1 = lon1[:-1]
+    lat1 = lat1[:-1]
+
+    with pytest.raises(ValueError, match="2 rows per scan"):
+        interpolate_geolocation_cartesian(
+            lon1, lat1, coarse_resolution=5000, fine_resolution=1000)