[pnuu]

Claude found an optimization on computing the projection coordinates. It now computes the destination projection coordinates only for every Nth pixel and then applies bilinear interpolation to get the dense grid.

Builds on top of #753.

Benchmark — SPEC-GS-004 and SPEC-GS-014 branches against main

Carried out on 2026-09-24.

Branch	Commit	Contents
main	aadaca1	Baseline
claude-per-block-gradient	062c72c	SPEC-GS-004: area source coordinates and gradients as 1D broadcast views, const kernel memoryviews
claude-coarse-dst-transform	1da2a61	SPEC-GS-014 prototype on top of GS-004: opt-in coarse target transform (precompute(transform_step=...))

Summary

• SPEC-GS-014 with transform_step=16 is 1.3–5.2× faster than main end to end, with
  indices within 0.003 px (LAEA) and 0.01 px (global lon/lat, equal to the tolerance) of
  main's.
• SPEC-GS-004 makes a single per-block call up to 6.8× faster and cuts its traced memory
  from 752 MiB to 32 MiB. In the full resampler its speed can't be distinguished from
  main's. It does cut peak memory by 26–45% for targets that cover the whole disc. All
  results are bit-identical to main.
• main silently drops about 0.2% of the valid disc pixels in the global case, because
  the search kernel does not converge there (SPEC-GS-006/GS-009). The GS-014 prototype
  happens to recover them.

Method

• Machine and software: 16-core Linux machine, pinned to 2 CPUs (taskset -c 0,1),
  OMP_NUM_THREADS=2, dask threaded scheduler with num_workers=2. Python 3.14.7,
  numpy 2.5.3, Cython 3.3.0, pyproj 3.8.0, PROJ 9.8.1.
• Builds: each branch was checked out in its own git worktree, with
  _gradient_search.pyx compiled in place, and selected with PYTHONPATH. Each run was a
  separate Python process.
• Source: SEVIRI full disc, geos, 3712×3712.
• Targets:
  • laea1024 / laea4096: EPSG:3035 over Europe, 1024² / 4096²;
  • geos1024: the source projection, 1024², part of the disc;
  • global4096: EPSG:4326, 4096×2048, 63% of it off the disc.
• Per-block case: gradient_resampler_indices(seviri, target) on the full-disc source,
  5 repeats. Memory is the tracemalloc peak of one call.
• End-to-end case: create_gradient_search_resampler(), precompute(), computing
  indices_xy, then compute() of float32 random data (1024² chunks), at CHUNK_SIZE
  4096 and 1024. Branches were interleaved over 4 rounds × 3 repeats, so each figure is the
  median of 12 runs. Memory is the peak RSS increase, sampled every 2 ms.
• Correctness: every branch's indices were compared with main's from the same case.

The machine is noisy. GS-004 and GS-014 without a step run effectively the same code, yet
their end-to-end medians differ by up to 20% and the interquartile ranges overlap widely.
Treat end-to-end differences below ~20% as noise.

Results

One per-block call, full-disc source

Target	main	GS-004	GS-014, exact	GS-014, step 16
LAEA 1024²	0.36 s, 752 MiB	0.19 s (1.9×), 32 MiB	0.19 s, 32 MiB	0.08 s (4.4×), 55 MiB
geos 1024² (same CRS)	0.19 s, 752 MiB	0.03 s (6.8×), 32 MiB	0.03 s, 32 MiB	0.04 s (5.2×), 51 MiB

GS-004 and GS-014-exact are bit-identical to main. At step 16 the maximum difference is
0.01 px (LAEA) and 0 px (geos, where the transform is the identity).

Full resampler, end to end

Median of 12 runs. "idx" is the time spent computing the indices; "resample" is the
bilinear resampling of the data.

Case	main	GS-004	GS-014, exact	GS-014, step 16
LAEA 4096², chunk 4096	3.93 s (idx 3.05)	4.20 s (idx 3.34)	3.90 s (idx 3.14)	1.50 s (2.6×) (idx 0.63)
LAEA 4096², chunk 1024	4.03 s (idx 3.40)	3.72 s (idx 3.13)	3.92 s (idx 3.27)	0.78 s (5.2×) (idx 0.33)
Global 4096×2048, chunk 4096	1.52 s (idx 1.08)	1.58 s (idx 1.05)	1.25 s (idx 0.82)	1.16 s (1.3×) (idx 0.68)
Global 4096×2048, chunk 1024	0.64 s (idx 0.44)	0.67 s (idx 0.43)	0.56 s (idx 0.36)	0.43 s (1.5×) (idx 0.25)

Peak RSS increase (median):

Case	main	GS-004	GS-014, exact	GS-014, step 16
LAEA 4096², chunk 4096	1510 MiB	1455 MiB	1475 MiB	1475 MiB
LAEA 4096², chunk 1024	598 MiB	614 MiB	604 MiB	586 MiB
Global 4096×2048, chunk 4096	1037 MiB	769 MiB	785 MiB	872 MiB
Global 4096×2048, chunk 1024	724 MiB	401 MiB	401 MiB	373 MiB

Differences from main's indices:

Case	GS-004	GS-014, exact	GS-014, step 16
LAEA 4096²	identical	identical	max 0.0032 px, same NaNs
Global 4096×2048	identical	identical	max 0.0100 px, NaNs differ (see below)

Findings

GS-004: large per-block win, small end-to-end effect

In the full resampler, resample_blocks gives each destination block only the source
crop it needs. For a European target that crop is small, so there was little gradient
work to remove, and the end-to-end time and memory don't change. For the global target
every block needs the whole disc, so the six source-sized float64 arrays that GS-004
removed were a real share of the peak: RSS drops by 26% (chunk 4096) and 45% (chunk 1024).
The large per-block gains only show up when a block's source crop is large.

GS-014: the destination transform was the bottleneck

With the target coordinates interpolated from a coarse grid, the index computation for
LAEA goes from ~3.1–3.4 s to 0.33–0.63 s. The global case gains less. The geos mapping
bends sharply near the edge of the disc, so the prototype has to refine that band down to
a step of 2, and the remaining cost there is mostly numpy bookkeeping over full-size masks.

main drops valid pixels in the global case

At step 16 the global case's NaN pattern differs from main's, but only in main's
disfavour. Step 16 returns 6,320 more valid pixels (of 3,095,692 valid ones) and loses
none. At those pixels the destination coordinates are bit-identical to the exact
transform, and they lie well inside the disc (e.g. source column 1835.6, row 51.5). main
returns NaN there because the search kernel gives up. The kernel warm-starts every pixel
from the previous pixel's solution, so tiny coordinate changes elsewhere on the path
decide whether it converges within its fixed budget of 5 steps. This is the behaviour
SPEC-GS-006 and SPEC-GS-009 describe, which until now had been established only by reading
the code: here main silently loses ~0.2% of the valid pixels.

The prototype recovering them is incidental, not a fix. A different target could just as
well lose pixels this way. The real fix belongs in the kernel.

Correction to SPEC-GS-004

SPEC-GS-004 states that np.broadcast_to views can't be passed to the non-const kernel
memoryviews. With numpy 2.5.3 and Cython 3.3.0, main's kernel does accept them; only
arrays explicitly marked read-only (flags.writeable = False) are rejected. The const
signatures are still the correct contract and needed for truly read-only input, but the
spec overstates what they were required for.

Reproducing

The benchmark was driven by two scripts: one runs a single case against the pyresample
on PYTHONPATH, and one loops over branches and cases pinned to 2 CPUs. For each branch:

git worktree add --detach <dir> <branch>
cp pyresample/ewa/*cpython-314*.so <dir>/pyresample/ewa/   # EWA extensions are unchanged
# compile <dir>/pyresample/gradient/_gradient_search.pyx in place (cythonize + numpy include)
PYTHONPATH=<dir> PYTROLL_CHUNK_SIZE=<chunk> taskset -c 0,1 python bench_branches.py <label> <case> [step]

Only main should save its index arrays for the comparison. A 4096² index array is
256 MB, and saving one per branch filled the 2 GB /tmp tmpfs on the first attempt.

• Closes #xxxx
• Tests added
• Tests passed
• Fully documented

[codecov]

Codecov Report

✅ All modified and coverable lines are covered by tests.
✅ Project coverage is 93.99%. Comparing base (4a46786) to head (b9a343b).
⚠️ Report is 17 commits behind head on main.

Additional details and impacted files

@@            Coverage Diff             @@
##             main     #754      +/-   ##
==========================================
+ Coverage   93.90%   93.99%   +0.09%     
==========================================
  Files          89       89              
  Lines       13924    14144     +220     
==========================================
+ Hits        13075    13295     +220     
  Misses        849      849              

Flag	Coverage Δ
unittests	93.99% <100.00%> (+0.09%)	⬆️

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.

[pnuu]

As far as I understand, the changes look reasonable. I'm still running tests with Satpy calls, but as Claude says the results are pretty noisy. I'll try by shutting down some other tasks and see if the results are more consistent.

[pnuu]

Here are some results with real FCI data, loading, creating, resampling and saving 16 composites to a 5624 x 5040 pixel EPSG:3035 (LAEA) area as COG geotiffs. First I ran the script three times for each of the three branches, building the cython extension when moving to the next one. As I noticed the first runs were always the fastest and the next runs got slower, and used more memory, I ran one additional run with each. These are the fourth number shown for each branch. Processing times are in seconds.

main
196.77108144760132
204.53654861450195
208.7838351726532
198.6278305053711

claude-per-block-gradient
166.24066615104675
188.39809799194336
194.8897659778595
154.08072543144226

claude-coarse-dst-transform
195.4351396560669
207.4542853832245
206.1188449859619
190.6227512359619

I don't understand what's going on. I also monitored the overall system CPU usage and disk read/write activity, and there was nothing notable happening there.

[pnuu]

I adjusted logging to emit only ERROR level and above, and fixed one (unused) import. The new results are mostly a bit more consistent. This time with 10 runs each. The script is run in a Bash loop, so a new interpreter is started for each run. I don't see any clear difference in maximum memory usage. Most likely it is dictated by the saving.

main
Elapsed time: 158.85512971878052
Elapsed time: 172.90649676322937
Elapsed time: 158.31317162513733
Elapsed time: 178.5110902786255
Elapsed time: 153.78597044944763
Elapsed time: 168.41702938079834
Elapsed time: 127.9962158203125
Elapsed time: 133.0429334640503
Elapsed time: 132.871568441391
Elapsed time: 147.755309343338

mean: 153.24549152851105
median: 156.04957103729248

claude-per-block-gradient
Elapsed time: 157.36019134521484
Elapsed time: 130.06735181808472
Elapsed time: 145.63860487937927
Elapsed time: 153.60946011543274
Elapsed time: 137.77070116996765
Elapsed time: 167.95250511169434
Elapsed time: 125.10808396339417
Elapsed time: 154.42376017570496
Elapsed time: 150.72444128990173
Elapsed time: 148.00337672233582

mean: 147.06584765911103
median: 149.36390900611877

claude-coarse-dst-transform
Elapsed time: 114.59825658798218
Elapsed time: 114.51854991912842
Elapsed time: 114.1517653465271
Elapsed time: 114.42341780662537
Elapsed time: 114.63398504257202
Elapsed time: 125.20915603637695
Elapsed time: 158.93143224716187
Elapsed time: 139.60028409957886
Elapsed time: 159.9323616027832
Elapsed time: 150.25194382667542

mean: 130.62511525154113
median: 119.92157053947449

There I ran some import timing tests that coincided with the increase in processing time. The tests were done during 2 or 3 of the runs, so I don't know if they had any actual effect as also the later runs had increased times 🙈

[pnuu]

Here's one dask profile, one of the fastest ones using this PR branch.

Another using main for reference:

The bouncy part is where the resampling happens. And maybe the memory usage in that part is also a bit lower.

[pnuu]

Now without scn.save_datasets() and instead doing da.compute(lcl_scn) after resampling. Note: this were running on battery power, so the computer might be running in some kind of powersaving mode resulting in slower times than the previous test.

main
Elapsed time: 166.5314645767212
Elapsed time: 163.60760140419006
Elapsed time: 166.1390540599823
Elapsed time: 166.32027912139893
Elapsed time: 164.83434295654297
Elapsed time: 165.5216464996338
Elapsed time: 165.65710830688477
Elapsed time: 165.02345633506775
Elapsed time: 165.60591316223145
Elapsed time: 165.80135226249695

Mean: 165.504221868515
Median: 165.6315107345581

claude-per-block-gradient
Elapsed time: 161.8393018245697
Elapsed time: 161.68171739578247
Elapsed time: 162.66494870185852
Elapsed time: 162.06241178512573
Elapsed time: 161.99148893356323
Elapsed time: 161.53293776512146
Elapsed time: 161.49555373191833
Elapsed time: 161.94949293136597
Elapsed time: 162.07496047019958
Elapsed time: 162.220276594162

Mean: 161.9513090133667
Median: 161.9704909324646

claude-coarse-dst-transform
Elapsed time: 161.61540293693542
Elapsed time: 160.37626123428345
Elapsed time: 161.98634028434753
Elapsed time: 161.99672389030457
Elapsed time: 161.77932143211365
Elapsed time: 162.20793223381042
Elapsed time: 162.8287959098816
Elapsed time: 162.2154746055603
Elapsed time: 161.93724155426025
Elapsed time: 162.87293100357056

Mean: 161.98164250850678
Median: 161.99153208732605

There is some variability in mem use between runs, but this branch has the lowest by my estimate looking through all the 30 graphs.

[djhoese]

I'm not so sure on this PR. The Claude summary talks about up to 5x timing improvement for some cases. Well, at least for the index generation step. So that should be ~3s of time savings for a calculation that your last profiles show take ~160s. It isn't easy to see the benefit in that case/profile. I started reviewing the code but it is a lot of code and I didn't understand the algorithm to start with so I think this will need to be a pure @mraspaud decision.

I do think any comparisons should be between your other PR and this PR, not with main since your other PR has so much memory savings.

[pnuu]

I agree, the comparison should between the two new PRs so that the benefits of this one can be seen. Frankly, the "real life" dask plots don't show much of a difference between the two.

[pnuu]

I asked for a new benchmark for end-to-end and individual public API funtions. TL;DR: the new feature is faster, but uses somewhat more memory.

Gradient search benchmarks: main vs claude-coarse-dst-transform

• main at 0291a40f, branch at b9a343bc (one commit on top of main: opt-in transform_step for precompute()).
• The Cython extensions were rebuilt with pip install -e . --no-deps after each checkout. _gradient_search.pyx is identical on both revisions.
• Machine: AMD Ryzen 7 PRO 8840U, 30 GB RAM, Linux; Python 3.14.7, numpy 2.5.3, dask 2026.8.0, pyproj 3.8.0 (PROJ 9.8.1), xarray 2026.7.0.
• dask threaded scheduler, 4 threads: the machine has 16 hardware threads, but the container has a cgroup CPU quota of 4 CPUs, which dask respects. Each scenario ran in a fresh process.
• Each timing is the median of 7 repeats after 1 warm-up, in seconds. Source data is random float32, persisted before timing. Ratios are main / branch, so above 1 means the branch is faster.
• Memory is measured in the same repeats, without a sampling thread. Before each repeat the script runs gc.collect() and malloc_trim(0), resets the kernel's peak-RSS counter (/proc/self/clear_refs = 5), and reads the current RSS. Peak is VmHWM after the step minus that baseline, i.e. the extra resident memory the step needed at its worst moment. Retained is the RSS right after the step returns minus the baseline. It includes the returned result (e.g. the persisted indices) and memory the allocator has not yet given back. Values are medians over the repeats, in MiB. Memory ratios are main / branch, so above 1 means the branch uses less.
• Benchmark script: bench_gradient.py (kept outside the repo, in the session scratchpad).

Summary

Time

• Default behaviour is unchanged. With the default precompute() the branch runs the same code as main.
  At the default chunk size every step longer than ~10 ms is within 0.96x-1.02x. At chunk size 1024 the
  sub-0.5 s steps are noisier, 0.90x-1.02x, and that includes the swath case, where the code paths are identical
  too. The indices and resampled data are bit-identical to main's.
• Index computation dominates. For area sources, indices_xy.persist() (the dask work behind precompute())
  is 70-80% of the end-to-end time. Creating the resampler and building the graphs take milliseconds. The
  interpolation step (compute() with ready indices) is the same on both revisions.
• transform_step speeds up the SEVIRI -> Europe 4096x4096 LAEA case a lot. The index computation is
  2.9-3.6x faster at the default chunk size and 2.9-4.7x faster at chunk size 1024. End to end it is
  2.1-3.4x faster (transform_step=16, bilinear: 4.31 s -> 1.87 s at chunk size 4096, 1.24 s -> 0.47 s at 1024).
• The gain is smaller when much of the target is off the source disc (global 0.1 deg: 1.1-1.3x end to end).
  The cells near the disc edge fall back to finer steps or to exact transforms. Here transform_step=4 is a bit
  faster than transform_step=16.
• Swath sources are unaffected. transform_step is ignored for them.

Memory

• Default behaviour: same memory as main. At the default chunk size the peaks are within 1%, except swath
  end-to-end nn (343 vs 373 MiB), where the code is identical. At chunk size 1024 the smaller peaks (75-370 MiB)
  vary by up to ~25% between the two revisions despite identical code. Treat differences of that size as
  allocator noise.

• transform_step raises the peak memory of the index computation:

  • Europe 4096x4096: 512 -> 698-746 MiB at chunk size 4096 (+36-46%), and 336 -> 400-406 MiB at 1024 (+19-21%).
  • Global target: 198 -> 421-461 MiB at chunk size 4096 (2.1-2.3x), and 100 -> 254-274 MiB at 1024 (2.5-2.7x).

  A full-resolution float64 array for one 4096x4096 chunk is 128 MiB, so the extra peak in the Europe case is
  about 1.5 such arrays per chunk. That fits the coarse path's per-chunk temporaries: the full-resolution
  interpolation result and the boolean pixel masks, held alongside dst_x/dst_y. The global target needs
  several refinement levels near the disc edge, which explains its larger relative increase. I have not
  profiled this line by line.

• End to end the extra memory is modest because the interpolation step peaks higher than the index
  computation. Europe: +4% (bilinear) and +7% (nn) at chunk size 4096, and +18% (bilinear) at 1024. Global:
  +13-15% (bilinear) and +24-29% (nn) at chunk size 4096.

• The interpolation step alone (compute() with persisted indices) also peaks 3-10% higher with transform_step
  results in the Europe case (both chunk sizes), even though the indices have the same shape and dtype. This is probably allocator
  state left behind by the earlier coarse-transform runs in the same process, not real extra work.

Accuracy

• The index error stays below the 0.01 source pixel tolerance: max 4e-4 px (transform_step=4) and
  6.3e-3 px (transform_step=16) for Europe, 1e-2 px for the global target. The large "max bilinear value diff"
  values come from the input being uniform random noise in [0, 100), whose pixel-to-pixel gradient is huge.
  Smooth real data would differ far less.
• A few NaN-pattern differences at the limb. For the global target, 14 pixels (28 x/y index elements; 16 at
  chunk size 1024) near +-80 deg latitude get valid indices on the branch but NaN on main. This is not off-disc
  leakage: the coarse transform's valid/invalid mask is identical to the exact transform's. On main, 4888 limb
  pixels have finite target coordinates but still come out NaN because the search fails right at the disc edge,
  and sub-0.01 px coordinate differences flip that outcome for a few of them. The branch indices there point at the
  outermost source row/column. For Europe at chunk size 1024 with transform_step=16, 2 pixels (4 elements)
  differ in the same way; I did not investigate them further. This contradicts the commit message's "same invalid
  pixels" for these targets, and is worth a note there or a test.

Timed steps (public API)

Step	What is timed
create_gradient_search_resampler	create_gradient_search_resampler(src, dst)
precompute (graph)	resampler.precompute(**kw); lazy, only builds the dask graph
indices (precompute + persist)	create + precompute() + indices_xy.persist()
indices persist only	indices_xy.persist(): the actual index computation
compute(method=...) (graph)	resampler.compute(data, method=...) with persisted indices; lazy
compute(method=...).compute() with persisted indices	Computing the interpolation alone
end-to-end resample(...).compute()	create + resample(data, method=..., **kw).compute()

PYTROLL_CHUNK_SIZE=4096 (default)

seviri_to_euro_laea_4k

SEVIRI full disc 3 km (geos, 3712x3712) -> Europe LAEA 4096x4096 (~1.5 km), target fully inside the disc

Time, seconds (median):

Step	main	branch default	ratio	branch transform_step=4	ratio	branch transform_step=16	ratio
create_gradient_search_resampler	0.000	0.000	0.97x	0.000	0.99x	0.000	1.06x
precompute (graph)	0.001	0.001	1.02x	0.001	1.01x	0.001	1.00x
indices (precompute + persist)	3.541	3.461	1.02x	1.213	2.92x	0.996	3.56x
indices persist only	3.468	3.463	1.00x	1.214	2.86x	1.001	3.46x
compute(method='bilinear') (graph)	0.002	0.002	0.99x	0.002	0.99x	0.002	0.98x
compute(method='bilinear').compute() with persisted indices	0.867	0.864	1.00x	0.856	1.01x	0.858	1.01x
compute(method='nn') (graph)	0.002	0.002	0.96x	0.002	0.96x	0.002	0.96x
compute(method='nn').compute() with persisted indices	0.242	0.241	1.00x	0.238	1.02x	0.239	1.02x
end-to-end resample(bilinear).compute()	4.308	4.324	1.00x	2.085	2.07x	1.872	2.30x
end-to-end resample(nn).compute()	3.706	3.700	1.00x	1.457	2.54x	1.249	2.97x

Memory, peak MiB above baseline (median), with retained MiB in parentheses:

Step	main	branch default	ratio	branch transform_step=4	ratio	branch transform_step=16	ratio
create_gradient_search_resampler	0 (0)	0 (0)	-	0 (0)	-	0 (0)	-
precompute (graph)	0 (0)	0 (0)	-	0 (0)	-	0 (0)	-
indices (precompute + persist)	512 (256)	512 (256)	1.00x	746 (256)	0.69x	698 (256)	0.73x
indices persist only	512 (256)	512 (256)	1.00x	746 (256)	0.69x	698 (256)	0.73x
compute(method='bilinear') (graph)	0 (0)	0 (0)	-	0 (0)	-	0 (0)	-
compute(method='bilinear').compute() with persisted indices	1136 (64)	1136 (64)	1.00x	1193 (170)	0.95x	1193 (170)	0.95x
compute(method='nn') (graph)	0 (0)	0 (0)	-	0 (0)	-	0 (0)	-
compute(method='nn').compute() with persisted indices	559 (64)	560 (64)	1.00x	617 (170)	0.91x	617 (170)	0.91x
end-to-end resample(bilinear).compute()	1392 (64)	1392 (64)	1.00x	1450 (170)	0.96x	1449 (170)	0.96x
end-to-end resample(nn).compute()	816 (64)	816 (64)	1.00x	874 (170)	0.93x	874 (170)	0.93x

Accuracy against main (exact transform), same chunk size:

Variant	index NaN-pattern mismatches (x and y planes counted separately)	max index diff (px)	p99 index diff (px)	bilinear NaN-pattern mismatches	max bilinear value diff
default	0	0.00e+00	0.00e+00	0	0.00e+00
transform_step=4	0	3.96e-04	3.89e-04	0	4.74e-02
transform_step=16	0	6.34e-03	6.08e-03	0	7.23e-01

seviri_to_globe_0.1deg

SEVIRI full disc 3 km (geos, 3712x3712) -> global EPSG:4326 0.1 deg (3600x1800), about half of the target outside the disc

Time, seconds (median):

Step	main	branch default	ratio	branch transform_step=4	ratio	branch transform_step=16	ratio
create_gradient_search_resampler	0.000	0.000	0.91x	0.000	0.98x	0.000	0.98x
precompute (graph)	0.001	0.001	0.99x	0.001	1.02x	0.001	1.02x
indices (precompute + persist)	0.703	0.695	1.01x	0.509	1.38x	0.556	1.26x
indices persist only	0.702	0.694	1.01x	0.511	1.37x	0.559	1.26x
compute(method='bilinear') (graph)	0.002	0.002	0.98x	0.002	1.00x	0.002	0.99x
compute(method='bilinear').compute() with persisted indices	0.276	0.282	0.98x	0.271	1.02x	0.282	0.98x
compute(method='nn') (graph)	0.002	0.002	1.00x	0.002	1.04x	0.002	1.01x
compute(method='nn').compute() with persisted indices	0.101	0.102	0.99x	0.107	0.94x	0.110	0.92x
end-to-end resample(bilinear).compute()	0.988	0.988	1.00x	0.795	1.24x	0.839	1.18x
end-to-end resample(nn).compute()	0.811	0.810	1.00x	0.630	1.29x	0.667	1.22x

Memory, peak MiB above baseline (median), with retained MiB in parentheses:

Step	main	branch default	ratio	branch transform_step=4	ratio	branch transform_step=16	ratio
create_gradient_search_resampler	0 (0)	0 (0)	-	0 (0)	-	0 (0)	-
precompute (graph)	0 (0)	0 (0)	-	0 (0)	-	0 (0)	-
indices (precompute + persist)	198 (99)	198 (99)	1.00x	421 (187)	0.47x	458 (250)	0.43x
indices persist only	198 (99)	197 (99)	1.00x	461 (250)	0.43x	458 (250)	0.43x
compute(method='bilinear') (graph)	0 (0)	0 (0)	-	0 (0)	-	0 (0)	-
compute(method='bilinear').compute() with persisted indices	432 (161)	432 (161)	1.00x	440 (70)	0.98x	451 (56)	0.96x
compute(method='nn') (graph)	0 (0)	0 (0)	-	0 (0)	-	0 (0)	-
compute(method='nn').compute() with persisted indices	210 (111)	210 (111)	1.00x	216 (70)	0.97x	228 (80)	0.92x
end-to-end resample(bilinear).compute()	588 (94)	588 (94)	1.00x	662 (168)	0.89x	676 (182)	0.87x
end-to-end resample(nn).compute()	371 (75)	371 (75)	1.00x	459 (193)	0.81x	478 (207)	0.77x

Accuracy against main (exact transform), same chunk size:

Variant	index NaN-pattern mismatches (x and y planes counted separately)	max index diff (px)	p99 index diff (px)	bilinear NaN-pattern mismatches	max bilinear value diff
default	0	0.00e+00	0.00e+00	0	0.00e+00
transform_step=4	28	1.00e-02	7.10e-03	14	1.32e+00
transform_step=16	28	1.00e-02	7.13e-03	14	1.32e+00

swath_to_nordic_laea_2k

Synthetic dask swath 2000x2000 (lons/lats) -> Nordic LAEA 2048x2048; transform_step does not apply to swath sources

Time, seconds (median):

Step	main	branch default	ratio
create_gradient_search_resampler	0.001	0.001	1.00x
precompute (graph)	0.001	0.001	0.98x
indices (precompute + persist)	0.445	0.451	0.99x
indices persist only	0.444	0.444	1.00x
compute(method='bilinear') (graph)	0.002	0.002	0.99x
compute(method='bilinear').compute() with persisted indices	0.207	0.217	0.96x
compute(method='nn') (graph)	0.002	0.002	1.00x
compute(method='nn').compute() with persisted indices	0.056	0.057	0.98x
end-to-end resample(bilinear).compute()	0.659	0.665	0.99x
end-to-end resample(nn).compute()	0.506	0.502	1.01x

Memory, peak MiB above baseline (median), with retained MiB in parentheses:

Step	main	branch default	ratio
create_gradient_search_resampler	0 (0)	0 (0)	-
precompute (graph)	0 (0)	0 (0)	-
indices (precompute + persist)	341 (217)	341 (217)	1.00x
indices persist only	341 (217)	341 (217)	1.00x
compute(method='bilinear') (graph)	0 (0)	0 (0)	-
compute(method='bilinear').compute() with persisted indices	287 (143)	285 (93)	1.01x
compute(method='nn') (graph)	0 (0)	0 (0)	-
compute(method='nn').compute() with persisted indices	157 (115)	155 (81)	1.01x
end-to-end resample(bilinear).compute()	503 (296)	504 (228)	1.00x
end-to-end resample(nn).compute()	343 (262)	373 (281)	0.92x

Accuracy against main (exact transform), same chunk size:

Variant	index NaN-pattern mismatches (x and y planes counted separately)	max index diff (px)	p99 index diff (px)	bilinear NaN-pattern mismatches	max bilinear value diff
default	0	0.00e+00	0.00e+00	0	0.00e+00

PYTROLL_CHUNK_SIZE=1024

seviri_to_euro_laea_4k

SEVIRI full disc 3 km (geos, 3712x3712) -> Europe LAEA 4096x4096 (~1.5 km), target fully inside the disc

Time, seconds (median):

Step	main	branch default	ratio	branch transform_step=4	ratio	branch transform_step=16	ratio
create_gradient_search_resampler	0.000	0.000	0.89x	0.000	0.99x	0.000	1.02x
precompute (graph)	0.008	0.008	0.99x	0.008	0.98x	0.008	1.00x
indices (precompute + persist)	0.958	1.015	0.94x	0.335	2.86x	0.209	4.59x
indices persist only	0.966	1.007	0.96x	0.324	2.98x	0.207	4.67x
compute(method='bilinear') (graph)	0.012	0.012	1.00x	0.012	1.01x	0.012	1.01x
compute(method='bilinear').compute() with persisted indices	0.280	0.283	0.99x	0.281	1.00x	0.280	1.00x
compute(method='nn') (graph)	0.012	0.012	1.01x	0.012	1.01x	0.012	1.02x
compute(method='nn').compute() with persisted indices	0.124	0.125	0.99x	0.127	0.98x	0.127	0.98x
end-to-end resample(bilinear).compute()	1.235	1.280	0.96x	0.592	2.09x	0.472	2.62x
end-to-end resample(nn).compute()	1.100	1.132	0.97x	0.437	2.51x	0.320	3.44x

Memory, peak MiB above baseline (median), with retained MiB in parentheses:

Step	main	branch default	ratio	branch transform_step=4	ratio	branch transform_step=16	ratio
create_gradient_search_resampler	0 (0)	0 (0)	-	0 (0)	-	0 (0)	-
precompute (graph)	0 (0)	0 (0)	-	0 (0)	-	0 (0)	-
indices (precompute + persist)	336 (336)	336 (336)	1.00x	406 (406)	0.83x	400 (400)	0.84x
indices persist only	352 (352)	352 (352)	1.00x	385 (385)	0.91x	373 (373)	0.94x
compute(method='bilinear') (graph)	0 (0)	0 (0)	-	0 (0)	-	0 (0)	-
compute(method='bilinear').compute() with persisted indices	375 (192)	383 (200)	0.98x	385 (190)	0.97x	389 (190)	0.96x
compute(method='nn') (graph)	0 (0)	0 (0)	-	0 (0)	-	0 (0)	-
compute(method='nn').compute() with persisted indices	192 (192)	200 (200)	0.96x	209 (209)	0.92x	208 (208)	0.92x
end-to-end resample(bilinear).compute()	392 (192)	391 (192)	1.00x	460 (210)	0.85x	462 (223)	0.85x
end-to-end resample(nn).compute()	344 (344)	355 (355)	0.97x	362 (362)	0.95x	368 (368)	0.94x

Accuracy against main (exact transform), same chunk size:

Variant	index NaN-pattern mismatches (x and y planes counted separately)	max index diff (px)	p99 index diff (px)	bilinear NaN-pattern mismatches	max bilinear value diff
default	0	0.00e+00	0.00e+00	0	0.00e+00
transform_step=4	0	3.96e-04	3.89e-04	0	4.74e-02
transform_step=16	4	6.34e-03	6.07e-03	2	7.23e-01

seviri_to_globe_0.1deg

SEVIRI full disc 3 km (geos, 3712x3712) -> global EPSG:4326 0.1 deg (3600x1800), about half of the target outside the disc

Time, seconds (median):

Step	main	branch default	ratio	branch transform_step=4	ratio	branch transform_step=16	ratio
create_gradient_search_resampler	0.000	0.000	1.02x	0.000	1.08x	0.000	1.06x
precompute (graph)	0.008	0.008	1.00x	0.008	1.02x	0.008	1.02x
indices (precompute + persist)	0.240	0.248	0.96x	0.196	1.22x	0.213	1.13x
indices persist only	0.232	0.244	0.95x	0.184	1.26x	0.206	1.13x
compute(method='bilinear') (graph)	0.011	0.010	1.02x	0.010	1.03x	0.010	1.04x
compute(method='bilinear').compute() with persisted indices	0.115	0.120	0.96x	0.118	0.98x	0.115	1.00x
compute(method='nn') (graph)	0.011	0.010	1.03x	0.010	1.03x	0.010	1.02x
compute(method='nn').compute() with persisted indices	0.054	0.058	0.94x	0.059	0.92x	0.056	0.97x
end-to-end resample(bilinear).compute()	0.322	0.352	0.91x	0.285	1.13x	0.292	1.10x
end-to-end resample(nn).compute()	0.268	0.285	0.94x	0.224	1.20x	0.244	1.10x

Memory, peak MiB above baseline (median), with retained MiB in parentheses:

Step	main	branch default	ratio	branch transform_step=4	ratio	branch transform_step=16	ratio
create_gradient_search_resampler	0 (0)	0 (0)	-	0 (0)	-	0 (0)	-
precompute (graph)	0 (0)	0 (0)	-	0 (0)	-	0 (0)	-
indices (precompute + persist)	100 (100)	100 (100)	1.00x	274 (274)	0.37x	254 (254)	0.40x
indices persist only	75 (44)	121 (90)	0.62x	275 (275)	0.27x	241 (241)	0.31x
compute(method='bilinear') (graph)	0 (0)	0 (0)	-	0 (0)	-	0 (0)	-
compute(method='bilinear').compute() with persisted indices	281 (91)	330 (169)	0.85x	325 (198)	0.86x	274 (232)	1.02x
compute(method='nn') (graph)	0 (0)	0 (0)	-	0 (0)	-	0 (0)	-
compute(method='nn').compute() with persisted indices	163 (153)	174 (145)	0.94x	199 (182)	0.82x	165 (165)	0.99x
end-to-end resample(bilinear).compute()	279 (187)	370 (178)	0.76x	354 (205)	0.79x	325 (228)	0.86x
end-to-end resample(nn).compute()	198 (198)	221 (213)	0.90x	274 (259)	0.72x	288 (282)	0.69x

Accuracy against main (exact transform), same chunk size:

Variant	index NaN-pattern mismatches (x and y planes counted separately)	max index diff (px)	p99 index diff (px)	bilinear NaN-pattern mismatches	max bilinear value diff
default	0	0.00e+00	0.00e+00	0	0.00e+00
transform_step=4	32	1.00e-02	7.13e-03	16	1.32e+00
transform_step=16	32	1.00e-02	7.16e-03	16	1.32e+00

swath_to_nordic_laea_2k

Synthetic dask swath 2000x2000 (lons/lats) -> Nordic LAEA 2048x2048; transform_step does not apply to swath sources

Time, seconds (median):

Step	main	branch default	ratio
create_gradient_search_resampler	0.001	0.001	0.96x
precompute (graph)	0.002	0.003	0.93x
indices (precompute + persist)	0.422	0.453	0.93x
indices persist only	0.420	0.457	0.92x
compute(method='bilinear') (graph)	0.004	0.004	0.95x
compute(method='bilinear').compute() with persisted indices	0.081	0.083	0.97x
compute(method='nn') (graph)	0.004	0.004	0.94x
compute(method='nn').compute() with persisted indices	0.040	0.041	0.97x
end-to-end resample(bilinear).compute()	0.467	0.520	0.90x
end-to-end resample(nn).compute()	0.448	0.487	0.92x

Memory, peak MiB above baseline (median), with retained MiB in parentheses:

Step	main	branch default	ratio
create_gradient_search_resampler	0 (0)	0 (0)	-
precompute (graph)	0 (0)	0 (0)	-
indices (precompute + persist)	727 (648)	746 (588)	0.97x
indices persist only	605 (587)	598 (560)	1.01x
compute(method='bilinear') (graph)	0 (0)	0 (0)	-
compute(method='bilinear').compute() with persisted indices	345 (345)	347 (305)	0.99x
compute(method='nn') (graph)	0 (0)	0 (0)	-
compute(method='nn').compute() with persisted indices	197 (197)	199 (198)	0.99x
end-to-end resample(bilinear).compute()	683 (642)	653 (559)	1.05x
end-to-end resample(nn).compute()	685 (674)	656 (580)	1.04x

Accuracy against main (exact transform), same chunk size:

Variant	index NaN-pattern mismatches (x and y planes counted separately)	max index diff (px)	p99 index diff (px)	bilinear NaN-pattern mismatches	max bilinear value diff
default	0	0.00e+00	0.00e+00	0	0.00e+00

[djhoese]

I think my main thing with this is that it is interpolating coordinates that we know the real values for. That fact that it is worse in memory makes me want it even less. Again, I'll let @mraspaud way in when he has time but I'm not a big fan so far. I'd rather find a way to speed up coordinate transformation for everything (better pyproj/PROJ performance somehow or better dask graph creation/tasks).

[pnuu]

I'm also leaning towards scrapping this. I guess I could give Claude some more cases to try and see if there are cases where the benefits are big, like FCI to a large target area. But looking at the large picture, saving a second or two from the resampling in a 3 minute process (like my operational FCI imager) doesn't make a big difference.

[pnuu]

I asked Claude to resample the SEVIRI example to FCI 500 m grid. Here the speed-up is noticeable. The values are not comparable to the earlier cases as the container has a 8 GB RAM limit, but might be more like what happens when we process data with Satpy and save as tiled geotiff.

If @mraspaud has time, I'd like him to test processing with SAR data or other very high resolution data, which I don't personally use.

SEVIRI -> MTG FCI FDSS 500 m target

SEVIRI full disc 3 km (geos, 3712x3712) -> mtg_fci_fdss_500m (geos, WGS84, h=35786400 m, 22272x22272, ~496 M pixels). The target corners are outside the Earth disc.

Same benchmark and machine as above: median of 7 repeats after 1 warm-up, memory as peak MiB above baseline with retained MiB in parentheses, ratios main / branch. Only timings and memory were recorded for this target, no accuracy comparison.

The container is limited to 8 GiB of memory (cgroup memory.max), and for this target the persisted indices alone are ~7.4 GiB (2 x 22272 x 22272 float64). A full .compute() of the 22272x22272 bilinear output also does not fit reliably: a probe at chunk size 4096 with transform_step=16 was OOM-killed. Every dask computation is therefore written chunk by chunk with da.store(..., lock=False) into a target that discards the data, so neither the indices nor the output is ever held in full. This is close to a pipeline that writes the result to disk chunk by chunk. As a result the steps differ from the tables above:

• indices (precompute + da.store to null): create + precompute() + computing all the indices.
• compute(method=...) (graph): resampler.compute(data, method=...) on the lazy, not persisted, indices.
• end-to-end resample(...), da.store to null: create + resample(data, method=..., **kw) + computing every output chunk.
• There are no "with persisted indices" steps. The interpolation cost is roughly the end-to-end time minus the index time.

PYTROLL_CHUNK_SIZE=4096 (default)

Time, seconds (median):

Step	main	branch default	ratio	branch transform_step=4	ratio	branch transform_step=16	ratio
create_gradient_search_resampler	0.000	0.000	0.78x	0.000	0.93x	0.000	0.86x
precompute (graph)	0.024	0.024	1.01x	0.024	0.99x	0.024	1.00x
indices (precompute + da.store to null)	22.625	27.401	0.83x	10.786	2.10x	8.729	2.59x
compute(method='bilinear') (graph)	0.030	0.031	0.99x	0.030	1.01x	0.030	1.01x
compute(method='nn') (graph)	0.031	0.031	1.00x	0.030	1.01x	0.030	1.01x
end-to-end resample(bilinear), da.store to null	29.979	37.229	0.81x	20.503	1.46x	18.583	1.61x
end-to-end resample(nn), da.store to null	27.142	30.902	0.88x	14.251	1.90x	12.326	2.20x

Memory, peak MiB above baseline (median), with retained MiB in parentheses:

Step	main	branch default	ratio	branch transform_step=4	ratio	branch transform_step=16	ratio
create_gradient_search_resampler	0 (0)	0 (0)	-	0 (0)	-	0 (0)	-
precompute (graph)	0 (0)	0 (0)	-	0 (0)	-	0 (0)	-
indices (precompute + da.store to null)	2047 (0)	2048 (0)	1.00x	2441 (475)	0.84x	2434 (488)	0.84x
compute(method='bilinear') (graph)	1 (1)	1 (1)	1.00x	1 (1)	0.97x	1 (1)	0.96x
compute(method='nn') (graph)	1 (1)	1 (1)	1.00x	1 (1)	1.06x	1 (1)	1.04x
end-to-end resample(bilinear), da.store to null	4590 (135)	4526 (165)	1.01x	5311 (490)	0.86x	5445 (488)	0.84x
end-to-end resample(nn), da.store to null	2701 (158)	2647 (162)	1.02x	2873 (490)	0.94x	2883 (551)	0.94x

PYTROLL_CHUNK_SIZE=1024

Time, seconds (median):

Step	main	branch default	ratio	branch transform_step=4	ratio	branch transform_step=16	ratio
create_gradient_search_resampler	0.000	0.000	1.04x	0.000	1.10x	0.000	1.20x
precompute (graph)	1.462	1.491	0.98x	1.531	0.95x	1.518	0.96x
indices (precompute + da.store to null)	26.294	27.360	0.96x	9.742	2.70x	6.417	4.10x
compute(method='bilinear') (graph)	1.603	1.644	0.98x	1.650	0.97x	1.623	0.99x
compute(method='nn') (graph)	1.597	1.653	0.97x	1.651	0.97x	1.625	0.98x
end-to-end resample(bilinear), da.store to null	33.473	35.092	0.95x	18.135	1.85x	15.083	2.22x
end-to-end resample(nn), da.store to null	30.606	31.144	0.98x	13.685	2.24x	10.677	2.87x

Memory, peak MiB above baseline (median), with retained MiB in parentheses:

Step	main	branch default	ratio	branch transform_step=4	ratio	branch transform_step=16	ratio
create_gradient_search_resampler	0 (0)	0 (0)	-	0 (0)	-	0 (0)	-
precompute (graph)	0 (0)	0 (0)	-	0 (0)	-	0 (0)	-
indices (precompute + da.store to null)	112 (1)	112 (1)	1.00x	210 (210)	0.54x	210 (203)	0.53x
compute(method='bilinear') (graph)	1 (1)	1 (1)	1.01x	1 (1)	0.93x	1 (1)	0.92x
compute(method='nn') (graph)	1 (1)	1 (1)	1.00x	1 (1)	1.02x	1 (1)	1.02x
end-to-end resample(bilinear), da.store to null	307 (204)	325 (233)	0.95x	336 (224)	0.91x	344 (208)	0.89x
end-to-end resample(nn), da.store to null	172 (172)	186 (186)	0.93x	215 (215)	0.80x	212 (205)	0.81x

Note on the branch default ratios: the default path runs the same code as main, but at chunk size 4096 the table shows it 0.81-0.88x as fast. The repeats within each run were tight (main 22.5-22.7 s and branch 27.4-27.5 s for the indices). A follow-up A/B run alternated main -> branch -> main -> branch, rebuilding the extensions every time. It timed the default index step at 28.0-28.5 s on both revisions (3 repeats each), so the gap came from the shared host getting slower between the two runs, not from the code. Comparisons between runs made at different times can be off by ~20%. Each branch variant ran in the same process right after branch default, so comparing within the branch columns is the more reliable view. At chunk size 4096 the indices take 27.4 s -> 10.8 s (transform_step=4) and 8.7 s (transform_step=16), and end-to-end bilinear takes 37.2 s -> 20.5 s and 18.6 s.

[pnuu]

This whole thing will be replaced with a caching solution that I PoC'd locally for benchmarking:

Index computation, in seconds (create + precompute() + persist()):

  ┌───────┬──────────────┬──────┬──────────────┬──────────────┬────────────────┐
  │ chunk │    target    │ main │  cache cold  │  cache warm  │ read from disk │
  ├───────┼──────────────┼──────┼──────────────┼──────────────┼────────────────┤
  │ 4096  │ Europe 4096² │ 3.47 │ 3.49 (0.99x) │ 0.24 (14.4x) │ 0.37 (9.5x)    │
  ├───────┼──────────────┼──────┼──────────────┼──────────────┼────────────────┤
  │ 4096  │ global 0.1°  │ 0.71 │ 0.73 (0.97x) │ 0.07 (10.1x) │ 0.13 (5.3x)    │
  ├───────┼──────────────┼──────┼──────────────┼──────────────┼────────────────┤
  │ 1024  │ Europe 4096² │ 0.99 │ 1.07 (0.92x) │ 0.11 (8.6x)  │ 0.17 (5.9x)    │
  ├───────┼──────────────┼──────┼──────────────┼──────────────┼────────────────┤
  │ 1024  │ global 0.1°  │ 0.25 │ 0.25 (0.97x) │ 0.04 (5.9x)  │ 0.06 (3.9x)    │
  └───────┴──────────────┴──────┴──────────────┴──────────────┴────────────────┘

I don't particularly like the implementation, so I'll start fresh (hopefully) tomorrow.