[Python] Why does my Kriging with External Drift (KED) fusion result deviate significantly...

I've encountered a tricky issue while using Kriging with External Drift (KED) to fuse satellite and rain gauge data, and I was hoping to get some advice.

The Setup:

• 
  Data Sources:
  
  1. Satellite rainfall data (as the external drift);
  2. Rain gauge measurements from ground stations.
• 
  Fusion Method: KED (Kriging with External Drift).
  

The Core Problem:

In most cases, the fusion results look reasonable.

However, at certain times—especially when a station observes localized, heavy rainfall, creating a large discrepancy with the satellite data—the KED result exhibits severe "overshooting".

Specifically, the peak value in the fused rainfall field becomes anomalously high, far exceeding the maximum value from either of the original input sources.

I've attached a typical image illustrating the problem, where:

• Max Satellite Rainfall: 7.8
  
• Max Station Rainfall (from IDW interpolation for comparison): 16.4

But the KED Fused Result Peak reaches 44.9

rainfall_comparison:



My Questions:

Is this "overshooting" an inherent flaw of the KED algorithm, especially when dealing with large local gradients between the primary variable (stations) and the external drift (satellite)?

Or, is it more likely an issue with my parameter settings (e.g., variogram model selection, nugget effect) or my code implementation?

Are there any recommended methods to suppress or avoid this numerical amplification effect during the fusion process?

This issue only occurs intermittently, but it significantly impacts the reliability of the results.

I would be very grateful to hear any analysis or advice from those with experience in this area.

My code: `

import numpy as np
from pykrige.uk import UniversalKriging

def ked_fusion(sat_data, sat_lat, sat_lon, gauge_values, station_coords,
variogram_model="spherical", variogram_parameters=None):
# Prepare station data
gauge_points = np.array(station_coords)
gauge_values = np.array(gauge_values) / 3600.0 # Convert mm/hr to mm/s

# Filter out invalid stations
valid_mask = ~np.isnan(gauge_values)
gauge_values = gauge_values[valid_mask]
gauge_points = gauge_points[valid_mask]

# Handle case with no valid stations
if gauge_values.size == 0:
return np.zeros_like(sat_data)

# Set default variogram parameters if not provided
if variogram_parameters is None:
variogram_parameters = {"sill": 0.1, "range": 1.0, "nugget": 1e-6}

# Perform KED fusion
uk = UniversalKriging(
x=gauge_points[:, 0], # Station longitudes
y=gauge_points[:, 1], # Station latitudes
z=gauge_values, # Station observations
variogram_model=variogram_model,
variogram_parameters=variogram_parameters,
drift_terms=["external_Z"], # External drift term
external_drift=sat_data, # Satellite data as external drift
external_drift_x=sat_lon[0, :], # Grid longitude vector
external_drift_y=sat_lat[:, 0], # Grid latitude vector
pseudo_inv=True, # Use pseudo-inverse for stability
exact_values=False # Allow non-exact matches
)

# Generate fused grid
z, _ = uk.execute("grid", sat_lon[0, :], sat_lat[:, 0])
return np.where(z < 0, 0, z) # Clip negative values to zero


`

Continue reading...