cs.CVApr 12, 2026

Neural Stochastic Processes for Satellite Precipitation Refinement

Authors: Shunya Nagashima, Takumi Bannai, Shuitsu Koyama, Tomoya Mitsui, Shuntaro Suzuki

Organizations: Neurogica Inc. · LTS, Inc. · ME-Lab Japan, Inc.

Read PDF
Ask AI
Cite

Copy citation

DOINot available for this paper

Abstract

Accurate precipitation estimation is critical for flood forecasting, water resource management, and disaster preparedness. Satellite products provide global hourly coverage but contain systematic biases; ground-based gauges are accurate at point locations but too sparse for direct gridded correction. Existing methods fuse these sources by interpolating gauge observations onto the satellite grid, but treat each time step independently and therefore discard temporal structure in precipitation fields. We propose Neural Stochastic Process (NSP), a model that pairs a Neural Process encoder conditioning on arbitrary sets of gauge observations with a latent Neural SDE on a 2D spatial representation. NSP is trained under a single variational objective with simulation-free cost. We also introduce QPEBench, a benchmark of 43{,}756 hourly samples over the Contiguous United States (2021--2025) with four aligned data sources and six evaluation metrics. On QPEBench, NSP outperforms 13 baselines across all six metrics and surpasses JAXA's operational gauge-calibrated product. An additional experiment on Kyushu, Japan confirms generalization to a different region with independent data sources.

Figures & tables

Figure 1: Motivation for the proposed Neural Stochastic Process (NSP). Satellite precipitation estimates provide gridded coverage but carry systematic biases, while gauges are sparse but locally accurate. NSP encodes satellite estimates, elevation, and gauges with the encoder EφE_{\varphi} , evolves a latent state through stochastic dynamics, and maps it to a calibrated precipitation estimate with the decoder DθD_{\theta} ; “Eval only” indicates that the radar reference is used only for comparison, not as a model input.
Figure 2: Architecture of Neural Stochastic Process (NSP). (a) During training, gauge observations are split into context and target sets: the black markers in the input maps denote context gauges encoded through the mask mt\mathbf{m}_{t} , while the right panel shows target gauges yt\mathbf{y}_{t} used by Lrec\mathcal{L}_{\mathrm{rec}} . The shared encoder maps consecutive inputs to zt\mathbf{z}_{t} and zt+1\mathbf{z}_{t+1} , and Ltrans\mathcal{L}_{\mathrm{trans}} aligns the Neural SDE prior pθ(zt+1∣zt)p_{\theta}(\mathbf{z}_{t+1}\mid\mathbf{z}_{t}) with the encoder posterior qφ(zt+1∣xt+1)q_{\varphi}(\mathbf{z}_{t+1}\mid\mathbf{x}_{t+1}) . (b) At inference, all available gauges are treated as context, and a single encoder–decoder pass predicts y^t\hat{\mathbf{y}}_{t} without SDE integration.
MethodRMSEr↓\text{RMSE}_{r}\downarrowMAEr↓\text{MAE}_{r}\downarrowRMSEg↓\text{RMSE}_{g}\downarrowMAEg↓\text{MAE}_{g}\downarrowrr,coll↑r_{r,\text{coll}}\uparrowFSSR↑\text{FSS}_{R}\uparrow
Quantile mapping Wood et al. (2004)4.073 ±\pm 0.1842.012 ±\pm 0.0630.885 ±\pm 0.0270.173 ±\pm 0.0060.288 ±\pm 0.0140.479 ±\pm 0.010
EMOS Gneiting et al. (2005)3.885 ±\pm 0.1271.844 ±\pm 0.0781.115 ±\pm 0.0960.179 ±\pm 0.0190.305 ±\pm 0.0130.478 ±\pm 0.029
XGBoost Chen and Guestrin (2016)3.741 ±\pm 0.1291.827 ±\pm 0.0771.062 ±\pm 0.0880.179 ±\pm 0.0200.284 ±\pm 0.0160.487 ±\pm 0.008
GWR Brunsdon et al. (1996)3.739 ±\pm 0.2301.627 ±\pm 0.0810.662 ±\pm 0.0210.155 ±\pm 0.0040.376 ±\pm 0.0160.229 ±\pm 0.017
Cokriging Cressie (1993)3.706 ±\pm 0.1211.789 ±\pm 0.0691.065 ±\pm 0.0810.220 ±\pm 0.0150.279 ±\pm 0.0160.461 ±\pm 0.009
GSMaP Kubota et al. (2020)3.638 ±\pm 0.1261.826 ±\pm 0.0761.017 ±\pm 0.0860.179 ±\pm 0.0200.288 ±\pm 0.0140.483 ±\pm 0.014
Table 1: Quantitative comparison of precipitation refinement methods. Values are mean ±\pm standard deviation over three-fold time-series cross-validation. The best score in each column is in bold .
Figure 3: Qualitative result at 03:00 UTC on March 16, 2025. NSP better matches the radar reference despite the satellite overestimation.
Figure 4: Regional zoom over the southeastern United States at 16:00 UTC on March 16, 2025.
MethodRMSEr ⁣↓\text{RMSE}_{r}\!\downarrowRMSEg ⁣↓\text{RMSE}_{g}\!\downarrowrr,coll ⁣↑r_{r,\text{coll}}\!\uparrow
GWR Brunsdon et al. (1996)1.7880.6550.389
CNP Garnelo et al. (2018a)1.7791.9530.437
EMOS Gneiting et al. (2005)1.7381.0820.169
Cokriging Cressie (1993)1.7151.0770.139
Linear regression1.6880.7680.356
Kriging Cressie (1993)1.6400.9730.114
Table 2: Kyushu (Japan) regional results.
 ModelRMSEr↓\text{RMSE}_{r}\downarrowMAEr↓\text{MAE}_{r}\downarrowRMSEg↓\text{RMSE}_{g}\downarrowMAEg↓\text{MAE}_{g}\downarrowrr,coll↑r_{r,\text{coll}}\uparrowFSSR↑\text{FSS}_{R}\uparrow
Architecture ablations
(1-i)w/o Neural SDE2.888 ±\pm 0.0391.487 ±\pm 0.0530.338 ±\pm 0.0660.062 ±\pm 0.0180.479 ±\pm 0.0200.509 ±\pm 0.022
(1-ii)Deterministic2.885 ±\pm 0.0801.492 ±\pm 0.0440.355 ±\pm 0.1470.060 ±\pm 0.0280.487 ±\pm 0.0330.495 ±\pm 0.012
(1-iii)w/o Satellite residual3.141 ±\pm 0.0891.613 ±\pm 0.0380.252 ±\pm 0.0240.041 ±\pm 0.0050.475 ±\pm 0.0180.153 ±\pm 0.047
(1-iv)Homoscedastic2.932 ±\pm 0.0621.492 ±\pm 0.0200.260 ±\pm 0.0160.052 ±\pm 0.0030.471 ±\pm 0.0230.445 ±\pm 0.012
(1-v)Log-space MSE3.554 ±\pm 0.1041.726 ±\pm 0.0480.405 ±\pm 0.0310.043 ±\pm 0.0060.475 ±\pm 0.0190.525 ±\pm 0.022
Table 3: Ablation study on CONUS. Mean ±\pm std over three folds; best in bold .
(a) RMSE by intensity
IntensityGSMaPGCNSPRel. improv.
(mm/h)   vs. GSMaPvs. GC
00.5920.3770.308++ 48.0%++ 18.3%
0–10.6930.4360.356++ 48.6%++ 18.3%
1–54.0112.5282.055++ 48.8%++ 18.7%
5–107.7255.5975.190++ 32.8%++ 7.3%
Table 4: Error analysis. (a) Pixel-level RMSE by radar precipitation intensity (fold 3, 8,759 samples). (b) Median spatial displacement (km) between predicted and radar exceedance masks, measured by the shift maximizing 2-D cross-correlation; lower is better.
Appendix figures & tables25 assets

Supplementary material from the paper’s appendix.

Appendix

PropertyValue
Temporal coverageJanuary 2021 – December 2025 (5 years)
Total hourly steps43,824 (43,756 after quality filtering)
Grid resolution0.1∘×0.1∘0.1^{\circ}\times 0.1^{\circ} ( ∼{\sim} 11 km)
Grid size260×590260\times 590 (153,400 pixels)
Spatial extent24.05∘24.05^{\circ} – 49.95∘49.95^{\circ} N, 124.95∘124.95^{\circ} – 66.05∘66.05^{\circ} W
Mean rainy fraction4.6% ±\pm 2.0% of radar pixels per sample
Appendix
Table 5: QPEBench dataset overview.
PercentileValue (mm/h)
P500.82
P903.55
P955.67
P9914.45
Max192.20
Appendix
Table 6: Distribution of non-zero radar precipitation across the full dataset.
ChannelTransform
Satellite (GSMaP)log⁡(1+x)\log(1+x)
Gauge observationslog⁡(1+x)\log(1+x)
Elevation (DEM)x / 2000x\,/\,2000
Context maskBinary {0,1}\{0,1\}
Radar (target)None (raw mm/h)
Appendix
Table 7: Input normalization.
HyperparameterValue
OptimizerAdamW ( β1=0.9\beta_{1}{=}0.9 , β2=0.999\beta_{2}{=}0.999 )
Learning rate3×10−33\times 10^{-3} (OneCycleLR)
Weight decay1.35×10−31.35\times 10^{-3}
Epochs4
Batch size4 per GPU ( ×8=32\times 8=32 effective)
Mixed precisionFP16
Appendix
Table 8: Hyperparameters of NSP.
FoldTrainValidationTest
1202120222023
22021–202220232024
32021–202320242025
Appendix
Table 9: Expanding-window time-series cross-validation splits over the 2021–2025 period.
MethodRMSEr↓\text{RMSE}_{r}\downarrowMAEr↓\text{MAE}_{r}\downarrowRMSEg↓\text{RMSE}_{g}\downarrowMAEg↓\text{MAE}_{g}\downarrowrr,coll↑r_{r,\text{coll}}\uparrow
ViT2.9792.0101.9931.2750.483
Quantile mapping1.8761.0731.3770.6230.216
GWR1.7880.8120.6550.2520.389
CNP1.7791.5591.9531.8150.437
EMOS1.7380.9541.0820.3260.169
Cokriging1.7150.9511.0770.3840.139
Appendix
Table 10: Full quantitative comparison on the Kyushu (Japan) regional dataset. Single train/validation/test split (2020/2021/2022); best in bold .
Threshold (mm/h)GSMaP GCNSPNote
0.133 km (3 px)33 km (3 px)Equal
1.033 km (3 px)55 km (5 px)GC better
5.0440 km (40 px)220 km (20 px)NSP 2×2\times better
10.0not reachednot reachedBoth below 0.5
Appendix
Table 11: Minimum useful scale: smallest neighborhood where FSS≥0.5\text{FSS}\geq 0.5 [ 52 ] . Lower indicates spatially sharper predictions.
ThresholdGCNSPImprov.NSP biasGC bias
(mm/h)cRMSEcRMSE   
0.13.663.54++ 3.3%−- 0.52−- 0.34
1.05.414.87++ 10.0%−- 0.92−- 0.41
5.09.729.24++ 4.9%−- 3.54−- 2.00
10.013.3811.88++ 11.2%−- 4.20−- 2.49
25.021.3214.70++ 31.1%−- 0.23++ 3.40
Appendix
Table 12: Conditional intensity RMSE (cRMSE) and bias, computed only on pixels where both prediction and radar exceed the threshold.
CategoryDefinitionMSE share
Hit bias (intensity)pred ≥0.1\geq 0.1 and radar ≥0.1\geq 0.167.0%
Miss (position)pred <0.1<0.1 and radar ≥0.1\geq 0.114.4%
False alarm (position)pred ≥0.1\geq 0.1 and radar <0.1<0.118.6%
Intensity-related 67.0%
Position-related 33.0%
Appendix
Table 13: MSE source decomposition into mutually exclusive pixel categories ( τ=0.1\tau{=}0.1 mm/h).
GC accuracyNNGC RMSENSP RMSEΔ\DeltaNSP wins
RMSE <1<1 (accurate)7,2360.5920.575−- 0.01752.7%
RMSE 1–21,4151.1931.014−- 0.17981.0%
RMSE 2–4342.3650.836−- 1.529100.0%
Overall8,686   57.5%
Appendix
Table 14: Per-sample comparison conditioned on GSMaP GC accuracy (fold 3). 73 samples with missing GC data are excluded; 1 sample with GC RMSE ≥4{}\geq 4 is included in Overall but omitted from bins.
 CountProportion
NSP better than GC4,99457.5%
NSP worse than GC3,69242.5%
Dominant error in 3,692 failure cases
Intensity error (hit bias)2,68172.6%
False precipitation57115.5%
Missed precipitation44011.9%
Appendix
Table 15: Per-sample comparison with GSMaP GC and failure mode decomposition ( τ=0.1\tau{=}0.1 mm/h).
ModeCountProportion
Intensity error (hit bias)7,12581.3%
False precipitation93010.6%
Missed precipitation7048.0%
Appendix
Table 16: Dominant error mode for all 8,759 test samples ( τ=0.1\tau{=}0.1 mm/h).
VariantRMSEr↓\text{RMSE}_{r}\downarrowMAEr↓\text{MAE}_{r}\downarrowRMSEg↓\text{RMSE}_{g}\downarrowMAEg↓\text{MAE}_{g}\downarrowrr,coll↑r_{r,\text{coll}}\uparrowFSSR↑\text{FSS}_{R}\uparrow
Matched variance2.883 ±\pm 0.1391.489 ±\pm 0.0760.341 ±\pm 0.0620.063 ±\pm 0.0150.478 ±\pm 0.0210.522 ±\pm 0.018
Naive temporal penalty2.853 ±\pm 0.1211.470 ±\pm 0.0890.366 ±\pm 0.0370.070 ±\pm 0.0110.478 ±\pm 0.0180.519 ±\pm 0.042
Full transition KL2.818 ±\pm 0.0621.444 ±\pm 0.0260.393 ±\pm 0.0470.076 ±\pm 0.0130.478 ±\pm 0.0210.527 ±\pm 0.022
Appendix
Table 17: Temporal objective variants on CONUS. The matched-variance row uses the Girsanov special case σφ2=σθ2Δt\sigma_{\varphi}^{2}=\sigma_{\theta}^{2}\Delta t ; the naive temporal penalty replaces the KL with mean-only temporal matching. Best values are in bold .
Figure 5: Hour-by-hour event tracking using the existing radar-evaluated metrics RMSEr↓\mathrm{RMSE}_{r}\downarrow and FSSR↑\mathrm{FSS}_{R}\uparrow on two 12-hour fold 3 events. NSP is compared against the no- Ltrans\mathcal{L}_{\mathrm{trans}} ablation and ConvCNP. In both events, removing Ltrans\mathcal{L}_{\mathrm{trans}} degrades RMSEr\mathrm{RMSE}_{r} and usually lowers FSSR\mathrm{FSS}_{R} , with the effect especially pronounced in the late-October event. ConvCNP improves over naive gauge-only interpolation but still fails to preserve thresholded spatial structure, yielding consistently low FSSR\mathrm{FSS}_{R} under dense-grid prediction.
EventModelRMSEr↓\mathrm{RMSE}_{r}\downarrowσ(RMSEr)↓\sigma(\mathrm{RMSE}_{r})\downarrowFSSR↑\mathrm{FSS}_{R}\uparrow
March 15, 2025NSP2.7700.3920.707
March 15, 2025NSP w/o Ltrans\mathcal{L}_{\mathrm{trans}}2.8280.4030.671
March 15, 2025ConvCNP3.1230.4270.069
October 30–31, 2025NSP1.8920.1910.539
October 30–31, 2025NSP w/o Ltrans\mathcal{L}_{\mathrm{trans}}2.2950.3230.472
October 30–31, 2025ConvCNP2.1660.2490.035
Appendix
Table 18: Event-wise temporal analysis on two 12-hour CONUS events in fold 3. Values are means over the 12 hourly steps; σ(RMSEr)\sigma(\mathrm{RMSE}_{r}) denotes the standard deviation of RMSEr\mathrm{RMSE}_{r} across the event window. ConvCNP denotes a reimplementation adapted to dense-grid prediction (Appendix D.3 ).
EventModeRMSEr↓\mathrm{RMSE}_{r}\downarrowFSSR↑\mathrm{FSS}_{R}\uparrowMean latency (ms)
March 15, 2025NSP (standard inference)2.7760.7132.60
March 15, 2025NSP (pure SDE rollout)3.2390.3512.22
March 15, 2025NSP (hybrid, α=0.5\alpha{=}0.5 )3.1020.4972.92
October 30–31, 2025NSP (standard inference)1.8910.5432.60
October 30–31, 2025NSP (pure SDE rollout)2.0430.3132.21
October 30–31, 2025NSP (hybrid, α=0.5\alpha{=}0.5 )1.8510.5032.93
Appendix
Table 19: SDE rollout analysis on the same two 12-hour fold 3 events. NSP (standard inference) denotes the deployed inference mode used throughout the paper, which re-encodes each hour independently. NSP (pure SDE rollout) propagates the latent state without re-encoding after initialization, and NSP (hybrid) linearly blends rolled-out and re-encoded latent states with α=0.5\alpha=0.5 . Values are means over the 12 hourly steps; latency is the mean per-hour inference time within each mode.
Figure 6: Effect of the inference-time context ratio on NSP performance (fold 3, model trained with 50% context). (a) Radar-evaluated metrics are best at 10–25% context and degrade with denser input, revealing a gauge over-interpolation effect. (b) Gauge-evaluated metrics improve monotonically with denser context, as expected.
Figure 7: Effect of the training-time context ratio on NSP performance (fold 3, seed 42). Each point retrains NSP from scratch with the indicated fixed context ratio; the 100% setting diverged and is omitted. At 75%, FSSR\text{FSS}_{R} collapses despite a slight RMSEr\text{RMSE}_{r} improvement, confirming that dense gauge supervision during training induces overfitting to point observations.
Figure 8: Effect of the fixed station-network ratio on NSP performance (fold 3, fixed 50% context split within each retained station subset). Unlike the earlier context-ratio sweep, this experiment first fixes the available gauge network and then partitions only that network into context and held-out targets. (a) Radar-evaluated metrics remain stable across 10–75% retained stations. (b) Gauge-evaluated metrics are shown for completeness, but their values are not directly comparable to the earlier sweeps because the evaluated station subset changes with the retained network.
Figure 9: Gauge context ablation for the widespread convective event shown in Figure 3 (03:00 UTC, March 16, 2025). Removing gauge context degrades RMSEr\text{RMSE}_{r} by 16.7%; the ∣Δ prediction∣|\Delta\,\text{prediction}| panel confirms that gauge corrections concentrate in active precipitation regions.
Figure 10: Gauge context ablation for a localized winter precipitation event (13:00 UTC, December 24, 2025). Removing gauge context changes RMSEr\text{RMSE}_{r} by only 1.5%; the spatial structure is predominantly satellite-driven.
MethodCategoryLatency (ms)Parameters
GSMaP GCOperational product0.04—
Linear regressionStatistical1.05—
IDWInterpolation1.88—
NSPDeep learning2.644.19M
Appendix
Table 20: Inference latency and model size. Measured on a single NVIDIA H200 GPU with batch size 1, excluding data loading and host–device transfer.
Figure 11: Two success cases, each showing a full CONUS view (top) and regional zoom (bottom). Case 1: A warm-front system over the southeastern US—GSMaP overestimates the spatial extent; GSMaP GC attenuates but introduces discontinuities; NSP produces a coherent field matching radar and preserves the frontal band orientation. Case 2: An intense squall line from Texas to the Great Lakes—GSMaP inflates the convective core far beyond the radar boundary; NSP recovers the sharp leading edge and trailing stratiform gradient.
Figure 12: Failure case at 13:00 UTC on January 22, 2024. Top: Full CONUS—radar detects moderate precipitation over the upper Midwest, but the satellite input is near-zero; without a spatial prior, NSP’s residual correction cannot reconstruct the field. Bottom: Regional zoom—GSMaP fails to detect shallow warm-rain precipitation; GSMaP GC produces diffuse, incoherent patches; NSP, dependent on the satellite residual, cannot recover the missing structure.
MethodRMSEr↓\text{RMSE}_{r}\downarrowMAEr↓\text{MAE}_{r}\downarrowRMSEg↓\text{RMSE}_{g}\downarrowMAEg↓\text{MAE}_{g}\downarrowrr,coll↑r_{r,\text{coll}}\uparrowFSSR↑\text{FSS}_{R}\uparrow
DDPM7.482 ±\pm 2.5346.441 ±\pm 2.3658.178 ±\pm 2.8147.605 ±\pm 2.5550.361 ±\pm 0.1610.060 ±\pm 0.039
NSP (Ours)2.818 ±\pm 0.0621.444 ±\pm 0.0260.393 ±\pm 0.0470.076 ±\pm 0.0130.478 ±\pm 0.0210.527 ±\pm 0.022
Appendix
Table 21: Conditional DDPM baseline vs. NSP on CONUS (mean ±\pm std over three-fold cross-validation). The DDPM performs far worse than NSP on every metric, with large cross-fold variance. All metrics follow the same protocol as Table 1 .
More in this topic

Neural Processes

Neural Processes learn distributions over functions from observed examples, enabling fast predictions with uncertainty for new inputs and tasks.

See all 56 papers

Explore similar work

CardsList

  1. Station2Radar: query conditioned gaussian splatting for precipitation field

    Feb 28, 2026Doyi Kim, Minseok Seo, Changick KimGaussian SplattingRemote Sensing

    Korea Advanced Institute of Science and Technology (KAIST)

  2. From Drops to Grid: Noise-Aware Spatio-Temporal Neural Process for Rainfall Estimation

    May 7, 2026Rafael Pablos Sarabia, Joachim Nyborg, Morten Birk +1Neural Processes

    Dept. Comp. Sc., Aarhus University & Cordulus · Cordulus · Dept. of Comp. Science, Aarhus University

  3. A plug-and-play generative framework for multi-satellite precipitation estimation

    May 14, 2026Yunfan Yang, Haofei Sun, Xiuyu Sun +7Remote SensingGenerative Modeling

    State Key Laboratory of Atmospheric Boundary Layer Physics and Atmospheric Chemistry (LAPC), Institute of Atmospheric Physics, Chinese Academy of Sciences, Beijing, China. · College of Earth and Planetary Sciences, University of Chinese Academy of Sciences, Beijing, China. · Shanghai Academy of Artificial Intelligence for Science (SAIS), Shanghai, China. +5