Filling Raster Gaps with Cokriging

tool

geostatistics

A cloud-covered satellite band, a sparse survey rasterized with mostly-empty cells: fill the gaps using a second, fully-populated band as a covariate, via cokriging instead of single-variable kriging. Runs on R; Python shown as reference only.

Published

September 20, 2026

The kriging article on this site predicts a value everywhere from scattered points, using only how far apart those points are. But a raster with gaps is a different problem: you already have a dense grid, most of it valid, and a second band of the same file that’s fully populated and correlated with the one that isn’t (elevation next to a cloud-gapped spectral band, say). Ignoring that second band and kriging the gaps from the target’s own neighbours alone throws away real information. Cokriging is kriging that uses it: fit a model of how the two bands vary together, not just how one varies with distance, then predict the missing values from wherever the covariate is available.

Prefer to skip straight to a result? The standalone tool runs the same computation: upload your own multi-band raster, pick the two bands, and it fills the gaps for you.

1. Upload your raster

A multi-band GeoTIFF (.tif/.tiff) with at least two bands: one you want filled (has real gaps, NoData cells), one that doesn’t. That second band is the covariate: elevation, another spectral band, anything correlated with the first.

mutable uploadedFiles = null

mutable uploadStatus = WebGeoDS.Upload.defaultStatus

2. Why a covariate band helps

If you’ve read the kriging article, you’ve already met the variogram: bin every pair of points by distance, average their squared difference within each bin (half of it, actually, which is why the axis says semivariance, not variance), fit a smooth curve to those bins. That curve’s nugget, sill, and range describe how one variable varies with itself over distance. Cokriging needs one more piece: the cross-variogram, the same idea applied to two variables at once. It asks how much the target’s value at one location tends to differ from the covariate’s value at another location, as a function of the distance between them.

Fitting a target’s variogram, a covariate’s variogram, and their cross-variogram together, consistently enough that the combined system stays mathematically valid (technically: positive definite), is a linear model of coregionalization (LMC). Once fit, predicting the target at a gap uses both: the target’s own nearby known values, and the covariate’s value right there, weighted by how strongly the two actually move together. Where the covariate is missing too, there’s nothing extra to weight in. Cokriging can’t fill that pixel any differently than plain kriging could, which is why this article (and the standalone tool) leaves those cells empty rather than guessing.

3. Shared map

sharedMapEl = sharedMap.element

legendEl = {
  const p = predictResult;
  const f = fitResult;
  const range = p ? p.raster : (f ? f.raster : null);
  if (!range) {
    const wrap = document.createElement("div");
    wrap.className = "webgeods-legend";
    wrap.hidden = true;
    return wrap;
  }
  return WebGeoDS.legend({ gradient: FILL_RAMP, minLabel: `${range.min}`, maxLabel: `${range.max}` });
}

base64ToBytes = (base64Str) => {
  const binary = atob(base64Str);
  const bytes = new Uint8Array(binary.length);
  for (let i = 0; i < binary.length; i++) bytes[i] = binary.charCodeAt(i);
  return bytes;
}

mutable fitResult = null

mutable predictResult = null

4. Compute: fit the model in R — runnable here; Python — reference only

R — runnable in this browser
Python — reference (a different technique, not this one)

Warning

This code doesn’t run on this page, and it isn’t the same algorithm. See why Python takes a different approach here.

# Reference only -- not executable in this browser tool, and NOT a
# port of the R cell to the left. GSTools' "collocated cokriging"
# assumes the covariate is known at every point you're predicting
# (true here -- the covariate band has no NoData either, at any pixel
# the target needs filling), and uses that single co-located value
# directly instead of fitting a full cross-variogram between the two
# variables. Simpler, and it asks less of the data, but it's a
# genuinely different method from the fit.lmc() coregionalization
# model the R cell above fits -- see this article's own note on why
# that distinction matters before treating the two as interchangeable.

aFit = WebGeoDS.getCellValue("cokriging-fit-r", Generators)
{
  if (aFit && aFit.raster) {
    const r = aFit.raster;
    await sharedMap.setRasterImage("cokriging-target-raw", {
      bounds: r.wgs84Bounds,
      width: r.width,
      height: r.height,
      values: new Float32Array(base64ToBytes(r.data).buffer),
      min: r.min,
      max: r.max,
      colorRamp: FILL_RAMP
    });
    mutable fitResult = aFit;
    // A freshly fit model invalidates any fill predicted from a
    // PREVIOUS one -- same reasoning the kriging article's own upload
    // handler uses when new data arrives.
    await sharedMap.removeRasterImage("cokriging-target-filled");
    mutable predictResult = null;
  }
}

statCardPart1 = {
  const f = fitResult;
  const box = document.createElement("div");
  box.className = "webgeods-stat-grid";
  const row = (label, value) => {
    const dt = document.createElement("div");
    dt.className = "webgeods-stat-label";
    dt.textContent = label;
    const dd = document.createElement("div");
    dd.className = "webgeods-stat-value";
    dd.textContent = value;
    box.append(dt, dd);
  };
  if (!f) {
    row("Result", "—");
    return box;
  }
  const s = f.summary;
  row("Target band", s.targetBand);
  row("Covariate band", s.covariateBand);
  row("Pixels used to fit", s.coFitPixels.toLocaleString());
  row("Nugget", s.nugget);
  row("Partial sill", s.partialSill);
  row("Range", s.range);
  row("Sample mean", s.sampleMean);
  row("Sample variance", s.sampleVariance);
  return box;
}

5. The target’s own variogram

This is the target band’s direct variogram, the same shape as the kriging article’s own chart: half the average squared difference within each distance bin, with the fitted curve’s nugget, sill, and range marked on it. It’s not the whole story. The actual cokriging fit above also needs the covariate’s own variogram and the cross-variogram between the two, neither shown separately here. But it’s the piece that carries over directly from ordinary kriging, and the one place this article’s chart can show something the reader already knows how to read.

variogramChart = {
  const f = fitResult;
  const width = 640, height = 220, margin = { top: 20, right: 20, bottom: 36, left: 46 };

  const svg = d3.create("svg")
    .attr("viewBox", [0, 0, width, height])
    .attr("width", "100%")
    .attr("height", height);

  if (!f) {
    svg.append("text")
      .attr("x", width / 2).attr("y", height / 2)
      .attr("text-anchor", "middle")
      .attr("fill", "#8a8378")
      .text("Run the R cell above to see the variogram.");
    return svg.node();
  }

  const s = f.summary;
  const nugget = s.nugget;
  // Both operands already rounded to 4dp on the R side -- re-rounding
  // the sum catches the odd binary-float artifact (e.g. "14.9639999999999999")
  // that adding two such numbers in JS can still produce.
  const sill = Math.round((s.nugget + s.partialSill) * 10000) / 10000;
  const range = s.range;

  const emp = f.variogram.empirical.dist.map((d, i) => ({ dist: d, gamma: f.variogram.empirical.gamma[i], npairs: f.variogram.empirical.npairs[i] }));
  const fit = f.variogram.fitted.dist.map((d, i) => ({ dist: d, gamma: f.variogram.fitted.gamma[i] }));

  const x = d3.scaleLinear()
    .domain([0, d3.max([...emp, ...fit, { dist: range }], (d) => d.dist)])
    .range([margin.left, width - margin.right]);
  const y = d3.scaleLinear()
    .domain([0, d3.max([...emp, ...fit, { gamma: sill }], (d) => d.gamma) * 1.1])
    .range([height - margin.bottom, margin.top]);

  svg.append("g")
    .attr("transform", `translate(0,${height - margin.bottom})`)
    .call(d3.axisBottom(x).ticks(6))
    .append("text")
    .attr("x", width - margin.right).attr("y", 32)
    .attr("fill", "#2a2117").attr("text-anchor", "end")
    .text("Distance");

  svg.append("g")
    .attr("transform", `translate(${margin.left},0)`)
    .call(d3.axisLeft(y).ticks(5))
    .append("text")
    .attr("x", -margin.left).attr("y", margin.top - 6)
    .attr("fill", "#2a2117").attr("text-anchor", "start")
    .text("Semivariance");

  const refLine = (x1, y1, x2, y2) => {
    svg.append("line")
      .attr("x1", x1).attr("y1", y1).attr("x2", x2).attr("y2", y2)
      .attr("stroke", "#8a8378").attr("stroke-dasharray", "4,3");
  };
  refLine(margin.left, y(nugget), width - margin.right, y(nugget));
  refLine(margin.left, y(sill), width - margin.right, y(sill));
  refLine(x(range), margin.top, x(range), height - margin.bottom);

  svg.append("path")
    .datum(fit)
    .attr("fill", "none")
    .attr("stroke", "#ab502b")
    .attr("stroke-width", 2)
    .attr("d", d3.line().x((d) => x(d.dist)).y((d) => y(d.gamma)));

  svg.append("g")
    .selectAll("circle")
    .data(emp)
    .join("circle")
    .attr("cx", (d) => x(d.dist))
    .attr("cy", (d) => y(d.gamma))
    .attr("r", 5)
    .attr("fill", "#42583c")
    .append("title")
    .text((d) => `distance ${d.dist}, semivariance ${d.gamma}, ${d.npairs} pairs`);

  const pillLabel = (px, py, anchor, text) => {
    const t = svg.append("text")
      .attr("x", px).attr("y", py)
      .attr("text-anchor", anchor)
      .attr("dominant-baseline", "middle")
      .attr("fill", "#5b5347").attr("font-size", 11).attr("font-weight", 600)
      .text(text);
    const box = t.node().getBBox();
    svg.insert("rect", () => t.node())
      .attr("x", box.x - 4).attr("y", box.y - 2)
      .attr("width", box.width + 8).attr("height", box.height + 4)
      .attr("rx", 3)
      .attr("fill", "#fffdf8").attr("stroke", "#d8cdb8").attr("stroke-width", 1);
  };

  // d3.create() built a DETACHED svg -- getBBox() inside pillLabel()
  // returns all zeros in that state, so it's attached off-screen for
  // these three calls only, then detached again immediately after.
  document.body.appendChild(svg.node());
  svg.style("position", "absolute").style("visibility", "hidden").style("left", "-9999px");

  pillLabel(margin.left + 6, y(nugget), "start", `nugget ${nugget}`);
  pillLabel(width - margin.right - 6, y(sill), "end", `sill ${sill}`);
  const rangeAnchor = x(range) > width - margin.right - 60 ? "end" : "start";
  const rangeDx = rangeAnchor === "end" ? -6 : 6;
  pillLabel(x(range) + rangeDx, height - margin.bottom - 12, rangeAnchor, `range ${range}`);

  svg.node().remove();
  svg.style("position", null).style("visibility", null).style("left", null);

  return svg.node();
}

Once the model above looks right, predict the gaps themselves:

R — runnable in this browser
Python — reference (a different technique, not this one)

Warning

Still a different method, continued. If GSTools’ collocated cokriging were used instead of the R cell’s LMC, this is the step that would actually predict. Same caveat as above.

# Reference only, for the same reason as the fit step above: GSTools'
# collocated cokriging predicts directly from the fitted model and the
# covariate's own values at the target locations, no separate
# "prediction dataset" the way gstat's predict() needs one built by
# hand -- a real structural difference from the R cell to the left,
# not just a syntax one.

bPredict = WebGeoDS.getCellValue("cokriging-predict-r", Generators)
{
  if (bPredict && bPredict.raster) {
    const r = bPredict.raster;
    await sharedMap.setRasterImage("cokriging-target-filled", {
      bounds: r.wgs84Bounds,
      width: r.width,
      height: r.height,
      values: new Float32Array(base64ToBytes(r.data).buffer),
      min: r.min,
      max: r.max,
      colorRamp: FILL_RAMP
    });
    mutable predictResult = bPredict;
  }
}

statCardPart2 = {
  const p = predictResult;
  const box = document.createElement("div");
  box.className = "webgeods-stat-grid";
  if (!p) {
    box.hidden = true;
    return box;
  }
  const row = (label, value) => {
    const dt = document.createElement("div");
    dt.className = "webgeods-stat-label";
    dt.textContent = label;
    const dd = document.createElement("div");
    dd.className = "webgeods-stat-value";
    dd.textContent = value;
    box.append(dt, dd);
  };
  const s = p.summary;
  row("Pixels filled", s.pixelsFilled.toLocaleString());
  if (s.crossValidationRmse !== null && s.crossValidationRmse !== undefined) {
    row("Cross-validation RMSE", s.crossValidationRmse);
  }
  if (s.pixelsStillEmpty > 0) {
    row("Still empty", `${s.pixelsStillEmpty.toLocaleString()} (no covariate data there either)`);
  }
  return box;
}

6. Why Python takes a different approach here

R’s gstat (built on sp/sf) has no WebAssembly build of its closest Python equivalent, pykrige. That’s the same gap the kriging article already documents, and it’s true again here for the same reason: no pykrige wheel published for Pyodide. But cokriging specifically has a second, more interesting wrinkle worth being precise about instead of glossed over: Python’s own GSTools library does offer cokriging (Simple and Intrinsic Collocated Cokriging), but that’s a genuinely different technique from the fit.lmc() linear model of coregionalization the R cell above fits, not just a different implementation of the same idea. Collocated cokriging assumes the covariate is known at every location you’re predicting, and uses that one co-located value directly, without fitting a full cross-variogram between the two variables the way an LMC does. It asks less of the data and skips a real modeling step, which is reasonable when its own assumption holds but not a drop-in substitute for what this article and the standalone tool actually compute. No verified reference code for it sits on this page: writing some without running it against a real GSTools install would risk teaching a technique this project hasn’t actually checked.

7. Where to next

Just needed the map? The standalone Raster Gap Filler: upload your own multi-band file, pick the two bands, done.
Only have one band, no covariate? Ordinary kriging on the raster’s own valid pixels still works: the standalone tool falls back to it automatically when you leave the covariate on “None.” See the kriging article for how that half of the pipeline works on its own.

Second tool in the GeoStatistics family.