Rewrite VII interpolator to use dask map_blocks

[djhoese]

This is an evolution of #134 but basically a whole rewrite for better performance and better results at the edge of scans. This is entirely Claude written with some intervention by me. Claude basically got rid of the scipy interp1d usage and wrote it manually. This concerned me, but I asked Claude to justify it (see bottom section).

CC @ameraner @sjoro @pepephillips

Here is Claude's summary:

Interpolate VII tie points per scan with dask map_blocks

This replaces the two xarray.DataArray.interp() passes in viiinterpolator with a small numpy kernel that bilinearly interpolates one scan at a time. Dask input goes through da.map_blocks. It builds on bugfix-vii-c-contiguous and removes that branch's _as_c_contiguous copy, which is no longer needed.

Nothing changes in the public API. tie_points_interpolation and tie_points_geo_interpolation keep their signatures and still return DataArrays. Those keep the tie points' dimension names, name, attrs and non-dimension coordinates.

How the interpolation works

A METimage scan is 24 pixel rows with 4 tie rows. The tie rows sit at pixel positions 0, 8 and 16 of the scan, plus a trailing edge row at 24. Scans are not contiguous because of the bowtie effect, so every pixel should depend only on the tie rows of its own scan. The old code got this by putting all tie rows on one evenly spaced grid and skipping the gap between scans. The new code reshapes the tie points into scans instead, so the interpolation cannot reach outside a scan.

Problems with the old implementation that this fixes

1. Missing (NaN) tie points leaked into valid pixels, including across scans.

   • Linear interpolation computes a + (b - a) * w. For a pixel exactly on a tie point, w is 0 or 1, so the other tie point should drop out. But NaN * 0 is still NaN.
   • scipy's interp1d puts a pixel on tie point i in the interval (i−1, i) with w = 1. A NaN at i−1 therefore blanks the pixel.
   • For the first row of a scan, i−1 is the previous scan's edge tie row.
   • The new kernel copies tie-aligned pixels straight from the tie point. Only the pixels actually interpolated from a missing tie point become NaN.
   • This also removes ≤3.6e-14° differences caused by a + (b - a) not always rounding back to exactly b.

   Real-world impact: on a CONUS granule whose first scan is missing, the old code blanked the first row of scan 1. pyresample's EWA builds each scan's footprint ellipses from the scan's first and last rows. With a NaN there, it falls back to a flat box footprint for the whole scan. In a Polar2Grid run (vii_865, EWA, fixed grid):

   • 242,103 grid pixels changed: median 1 level, 99th percentile 11 levels, max 210.
   • 26,533 grid pixels that were empty are now filled.
   • 488 grid pixels smeared into by the fallback boxes are now masked.
   • All of these are within 55 grid cells of scan 1. Everything else is bit-identical.

2. Chunking did not follow scans.

   • .interp() merges the interpolated dimension into one chunk and splits the other. The old output came back as rows of 2052 pixels (85.5 scans), whatever the input chunks were.
   • The new output has one chunk per tie-point chunk, of whole scans and the full swath width. With Satpy's reader that is exactly its pixel grid, so _rechunk_to_pixel_grid becomes a no-op: 0 rechunk layers instead of 5.
   • Input that isn't chunked in whole scans and the full width is rechunked first, the same way as scanline_mapblocks does for MODIS.

3. Blocks were not C-contiguous.

   • Interpolating along the last dimension produced Fortran-ordered blocks.
   • pyresample ≤ 1.35.0 then fails in EWA ll2cr with ValueError: ndarray is not C-contiguous. This hit Polar2Grid users with orthorectification off.
   • Blocks are now C-contiguous by construction, with no extra copy.

4. tie_points_geo_interpolation computed dask input while building the graph, twice, to choose between geodetic and cartesian interpolation. (a or b on a 0-d DataArray evaluated the first reduction twice.) Satpy paid for this in the METimage file handler's __init__ for every granule.

   • Now the two reductions (max |lat| and longitude range) stay lazy. They are passed to map_blocks as 0-d dask arrays.
   • Dask computes them once per compute and hands the same two numbers to every task, so the choice stays global and chunking can't affect it.
   • Neither function computes anything while building the graph. Creating a Satpy Scene for this granule goes from 6 dask computes to 4, and from 0.093 s to 0.052 s. The remaining 4 are Satpy's own reads of the calibration coefficients.
   • The tie points are also read once instead of twice.
   • Trade-off: no pixel chunk starts until every tie-point chunk has been reduced. The tie points are about 64× smaller than the pixels, and this didn't show up in the timings.

5. Output dtype was inconsistent.

   • float32 tie points came back as float64 for numpy input but as float32 for dask input, and the geo function returned float64.
   • Now floating-point input keeps its dtype and integer input becomes float64, for numpy and dask alike. This is a behavior change for float32 numpy input; see the decisions below.

6. Smaller fixes:

   • The ValueError message about scan multiples was never formatted, because the format argument was passed as a second argument.
   • The geo docstring had duplicated Args/Returns sections.
   • The module docstring described a "per granule, not per scan" approximation that no longer applies.

Results are unchanged otherwise. The geodetic path is bit-identical to before, and the cartesian path is within 3.6e-14°. Results were already independent of chunking; that is now enforced by tests.

Performance

These numbers come from one real granule: 1368×394 tie points interpolated to 8208×3144 pixels. Test setup:

• 22 cores; xarray 2026.7.0, dask 2026.7.1, numpy 2.5.2.
• Best of 3 runs; "1 thread" means dask's synchronous scheduler.
• Peak memory is from tracemalloc.
• The tie-point chunks are the ones Satpy's reader produces at Polar2Grid's 75MiB and 32MiB settings.

tie_points_geo_interpolation	Wall time, old → new	1-thread time, old → new	Peak MiB, old → new	Dask tasks
lon/lat, 75MiB chunks	0.39 → 0.31 s	0.70 → 0.42 s	1675 → 890	43 → 22
lon/lat, 32MiB chunks	0.47 → 0.23 s	0.79 → 0.42 s	1479 → 847	111 → 44
cartesian, 75MiB chunks	1.05 → 1.04 s	2.11 → 1.27 s	2513 → 1242	114 → 22
cartesian, 32MiB chunks	1.02 → 0.67 s	2.26 → 1.27 s	1823 → 1361	284 → 44
lon/lat, numpy input	0.52 → 0.25 s	–	1009 → 443	–
cartesian, numpy input	1.88 → 1.11 s	–	1809 → 1430	–

Times include building the graph, so the old version's decision compute is counted. Of the new version's tasks, 12 (75MiB) or 26 (32MiB) are the small lazy reductions.

• Satpy metimage_l1b_nc file handler lon/lat (75MiB, cartesian because the granule reaches 68°N): wall time 0.97 → 0.99 s, 1-thread time 2.01 → 1.28 s, peak memory 2512 → 1242 MiB, 118 → 22 tasks.
• Full Polar2Grid run (2 bands, EWA): 11.5 s → 11.3 s. Geolocation is a small part of the total, so there's no measurable change.
• Why wall time barely improves at 75MiB: Satpy's reader chunks this granule into 6240 and 1968 rows, so one task carries 76% of the work. The old graph did its heavy work on even 4104-row chunks before Satpy rechunked. With several granules or more even chunks, the lower CPU time turns into lower wall time, as in the 32MiB rows.
• Where the gains come from:
  • Along-track interpolation runs on tie-point-sized arrays.
  • Everything is written in place into the output.
  • The cartesian back-conversion no longer evaluates both latitude formulas for every pixel; np.where did. The high-latitude formula is now only applied where it's needed.

Decisions still to make

1. Geodetic or cartesian interpolation.
   • The old global rule is kept: cartesian if any |lat| > 60° or the longitude span exceeds 180°. It is now evaluated lazily, so it no longer costs a compute.
   • Deciding per chunk was rejected: results would depend on chunk size. The chunk-independence tests catch that.
   • The remaining question is whether the rule itself is right.
     • This CONUS granule (28.6–68.1°N) goes cartesian for every pixel because of its northern end. The rule is per call, which in Satpy means per granule, so different granules in one scene can use different methods (as before).
     • Always using cartesian would be simpler and consistent across granules. It costs about 3× the CPU of the lon/lat path (1.27 s vs 0.42 s per granule on one thread).
   • On this granule the two methods differ by at most 22 m (99.9th percentile 16 m), about 4% of a 500 m pixel. Which one is closer to the truth has not been checked; neither is exact great-circle interpolation.
2. dtype policy. Is keeping float32 output for float32 input OK for all users? Previously numpy float32 input came back as float64. Satpy's METimage tie points are float64, so Satpy output is unaffected.
3. Chunk size policy. Each task holds a whole chunk of pixels: x, y, z and the outputs for the cartesian path. The function only rechunks input that isn't whole scans and the full width, and otherwise leaves chunk size to the caller. Should it also split very large chunks for more parallelism and less memory per task?
4. Rechunking non-conforming input. It currently mirrors MODIS: the first chunk's row count rounded down to whole scans (at least one scan), full width. An alternative is to size chunks with dask's automatic chunk size applied to the pixel output.

Follow-ups elsewhere

• Satpy metimage_nc.py:
  • After a release, _rechunk_to_pixel_grid and its docstring (about xarray.interp dividing chunks evenly) can go.
  • The reader could use more even chunk sizes; 6240 + 1968 rows limits parallelism for single granules.
• pyresample EWA: a NaN in the first or last row of a scan silently turns the whole scan into box footprints. Worth documenting, or warning about.

Testing

• test_viiinterpolator.py has 52 tests, each run with numpy input and with dask input. They cover:
  • C-contiguous blocks, and one scan-aligned, full-width chunk per tie-point chunk.
  • Laziness: neither function may compute while building the graph.
  • Bit-identical numpy and dask results across 6 chunkings, including partial scans and partial width that need rechunking. In the cartesian case only some chunks reach 60° on their own, so a per-chunk decision would fail (checked).
  • NaN handling: a missing scan, edge and middle tie rows, the first tie row of the next scan, and a tie column.
  • Exact values at tie points, and dtype and metadata behavior.
• 31 of the 52 fail against the previous implementation: the NaN cases, float32 dtype, and compute counts.
• The full suite passes, including doctests. flake8 is clean.
• Checked on a real METimage granule against the old implementation: through Satpy's reader, and end to end in Polar2Grid (results above).

🤖 Generated with Claude Code

Rewrite of interpolation kernel

Where both are valid, the values are bit-identical.

How to justify it to the original authors:

1. It isn't a new algorithm. It's the same linear formula scipy uses, y0 + (y1 − y0) × (x − x0)/(x1 − x0). Here every pixel sits at a fixed fraction (k/8) between two tie points, so the weights are constants and no interval search is needed. The output is bit-identical to scipy's, and the authors' own expected-value tables (TEST_LON_1…TEST_LAT_3) still pass.
2. Calling scipy per block brings back both bugs.
   • The cross-scan NaN leak comes back: scipy assigns a pixel sitting exactly on a tie point to the interval on its left.
   • The output comes back Fortran-ordered, which needs an extra copy.
   • Avoiding those means special handling around scipy anyway.
3. np.interp doesn't fit. It only works on 1-D arrays, so it would need a Python loop over thousands of rows and columns.
4. It's faster and uses less memory. It skips the interval search and writes straight into the output instead of building temporary arrays: 1.4× faster than scipy per block, or 1.8× once scipy's output is copied to C order.
5. It's what they originally wanted. Their V2 rewrite (bfe48f4, "without internal loop over scans") gave up per-scan interpolation for speed. Reshaping into scans gets exact per-scan interpolation back without a Python loop.

Suggested paragraph for the PR:

Why not scipy/numpy per block? Tie points and pixels are on fixed regular grids, so every pixel sits at a known fraction (k/tie_points_factor) between two tie points, and linear interpolation reduces to a constant-weight blend. _interpolate_intervals applies exactly the formula interp1d(kind="linear") uses (y0 + (y1 - y0) * w) and gives bit-identical values, with no interval search. Calling interp1d per block instead is about 1.4× slower (1.8× with the copy it would need). It returns Fortran-ordered arrays, and it keeps the zero-weight NaN leak across scans, because it assigns a pixel sitting exactly on a tie point to the interval on its left. np.interp only handles 1-D arrays.

• Closes #xxxx
• Tests added
• Tests passed
• Passes git diff origin/main **/*py | flake8 --diff
• Fully documented

[codecov]

Codecov Report

✅ All modified and coverable lines are covered by tests.
✅ Project coverage is 90.54%. Comparing base (0a01996) to head (523e6ed).
⚠️ Report is 1 commits behind head on main.

Additional details and impacted files

@@            Coverage Diff             @@
##             main     #135      +/-   ##
==========================================
+ Coverage   89.74%   90.54%   +0.79%     
==========================================
  Files          20       20              
  Lines        1541     1671     +130     
==========================================
+ Hits         1383     1513     +130     
  Misses        158      158              

Flag	Coverage Δ
unittests	90.54% <100.00%> (+0.79%)	⬆️

Flags with carried forward coverage won't be shown. Click here to find out more.

☔ View full report in Codecov by Harness.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

[coveralls]

coverage: 89.691% (+0.8%) from 88.861% — djhoese:feat-vii-map-blocks into pytroll:main