🔧 Fix Spatial machine learning Raster & rasterio Spatial Machine Learning · step 10 of 14 Updated

My Prediction Raster Is Striped, Blocky or Full of NoData

Problem statement

A prediction map comes out visibly wrong in a way the metrics did not show: horizontal stripes, hard rectangular edges, a flat region, or holes where there should be data.

Five distinct causes, each with a characteristic appearance:

symptom cause
horizontal or vertical stripes a covariate has stripes — usually flight lines or scene edges
hard rectangular blocks block-wise prediction with neighbourhood features and no buffer
a flat, featureless region the model extrapolating and saturating
holes one covariate has nodata there
a scrambled, noisy map features assembled in the wrong order, or a bad reshape

Only the last is a coding error. The others are the model or the data telling you something.

Quick answer

Check the inputs before the output:

import numpy as np

for name, layer in layers.items():
    finite = np.isfinite(layer)
    print(f"  {name:16} {1 - finite.mean():6.2%} nodata  "
          f"range {np.nanmin(layer):10.2f}..{np.nanmax(layer):10.2f}")

valid = np.stack([np.isfinite(l) for l in layers.values()]).all(axis=0)
print(f"  cells with every feature: {valid.mean():.1%}")

The prediction's holes are the union of the covariates' holes. If one layer is 20% nodata, the map is at most 80% complete no matter what the model does.

Five prediction artefacts: stripes from a covariate, blocks from unbuffered processing, a flat extrapolation region, holes from nodata, and a scrambled map.
Each artefact has a shape, and the shape identifies the cause.

Step-by-step solution

1. Stripes: look at the covariates

A model cannot invent structure. Stripes in the output are stripes in an input — a lidar flight line, a satellite scene boundary, a sensor artefact, or a mosaic seam.

Plot every covariate before modelling. A feature with visible stripes will produce a striped map, and the model will have learned the stripes as signal.

The fix is upstream: mosaic the covariate properly, normalise per flight line, or drop the feature.

2. Hard rectangular blocks: the processing tiles

Predicting block by block is correct for cell-wise features and wrong for neighbourhood features. A moving-window mean at a block edge sees only the cells inside the block.

buffer = window_size // 2 + 1
read_window = Window(col - buffer, row - buffer,
                     width + 2 * buffer, height + 2 * buffer)

Read buffered, compute, write the unbuffered part. The block size in the artefact matches the processing block, which is the tell.

3. Flat regions: the model is extrapolating

A random forest predicts the mean of training targets in a leaf. Where the covariates are outside anything it saw, it reuses its most extreme splits and the output converges on a constant.

The region is smooth, in range and information-free. Check it against the applicability mask — flat regions and low applicability usually coincide exactly.

4. Holes: one covariate has nodata

The valid mask is the intersection of all covariates' valid areas. A single layer with gaps punches those gaps through the whole prediction.

Options, in order of preference: fix the layer, drop the layer, or fit a second model without it and use that where the first cannot predict — recording which model produced each pixel.

5. A scrambled map: check the round trip 

direct = model.predict(stack[:, rows, cols].T)
from_raster = predicted[rows, cols]
assert np.allclose(direct, from_raster, atol=1e-4)

A reordered feature stack or a mismatched reshape produces a map that is spatially structured — because the inputs are — and completely wrong. Only a direct comparison catches it.

A neighbourhood feature computed per block without a buffer, producing discontinuities at every block boundary.
The artefact's block size matches the processing block size exactly, which identifies the cause.

Code examples

Example 1 — auditing the covariates before modelling 

import numpy as np


def covariate_audit(layers, stripe_threshold=3.0):
    """Nodata, ranges, and a stripe test on every input layer."""
    problems = []
    valid_stack = []

    for name, layer in layers.items():
        finite = np.isfinite(layer)
        valid_stack.append(finite)
        nodata = 1 - finite.mean()

        row_means = np.nanmean(layer, axis=1)
        col_means = np.nanmean(layer, axis=0)
        row_variation = np.nanstd(np.diff(row_means))
        col_variation = np.nanstd(np.diff(col_means))
        overall = np.nanstd(layer)

        flags = []
        if nodata > 0.05:
            flags.append(f"{nodata:.1%} nodata")
            problems.append(f"{name}: {nodata:.1%} nodata")
        if overall > 0 and row_variation / overall > stripe_threshold / 100:
            flags.append("row striping")
            problems.append(f"{name}: possible horizontal striping")
        if overall > 0 and col_variation / overall > stripe_threshold / 100:
            flags.append("column striping")
            problems.append(f"{name}: possible vertical striping")

        print(f"  {name:18} {nodata:6.2%} nodata  "
              f"{np.nanmin(layer):9.2f}..{np.nanmax(layer):9.2f}  "
              f"{' '.join(flags)}")

    valid = np.stack(valid_stack).all(axis=0)
    print(f"  cells with every feature: {valid.mean():.1%}")
    for name, layer in layers.items():
        unique = np.isfinite(layer) & ~valid
        if unique.mean() > 0.01:
            print(f"  ! {name} alone removes {unique.mean():.1%} more cells")
    return valid, problems

The last loop is the useful one: it identifies which single layer is responsible for the holes, which is the layer to fix or drop.

Example 2 — buffered block prediction 

import numpy as np
import rasterio
from rasterio.windows import Window
from scipy.ndimage import uniform_filter


def predict_blocks_buffered(model, layer_paths, feature_names, out_path,
                            block=512, neighbourhood=27):
    """Block-wise prediction that is correct for neighbourhood features."""
    buffer = neighbourhood // 2 + 1
    with rasterio.open(layer_paths[feature_names[0]]) as src:
        profile = src.profile
        height, width = src.height, src.width
    profile.update(count=1, dtype="float32", nodata=np.nan,
                   compress="deflate", tiled=True,
                   blockxsize=block, blockysize=block)

    sources = {n: rasterio.open(p) for n, p in layer_paths.items()}
    try:
        with rasterio.open(out_path, "w", **profile) as dst:
            for row in range(0, height, block):
                for col in range(0, width, block):
                    h = min(block, height - row)
                    w = min(block, width - col)

                    r0 = max(row - buffer, 0)
                    c0 = max(col - buffer, 0)
                    r1 = min(row + h + buffer, height)
                    c1 = min(col + w + buffer, width)
                    read = Window(c0, r0, c1 - c0, r1 - r0)

                    stack = np.stack([
                        sources[n].read(1, window=read).astype("float32")
                        for n in feature_names])
                    # neighbourhood features computed on the buffered block
                    extra = uniform_filter(np.nan_to_num(stack[0]),
                                           size=neighbourhood)
                    stack = np.concatenate([stack, extra[None, ...]])

                    valid = np.isfinite(stack).all(axis=0)
                    out = np.full(valid.shape, np.nan, dtype="float32")
                    if valid.any():
                        out[valid] = model.predict(stack[:, valid].T)

                    trim = out[row - r0:row - r0 + h, col - c0:col - c0 + w]
                    dst.write(trim, 1, window=Window(col, row, w, h))
    finally:
        for src in sources.values():
            src.close()
    print(f"  wrote {out_path} with a {buffer}-cell buffer")
    return out_path

Trimming the buffer before writing is what removes the block edges. Skipping the trim writes overlapping regions and produces a different artefact — visible seams where blocks overwrite each other.

Example 3 — locating the flat regions 

import numpy as np
from scipy.ndimage import uniform_filter


def find_flat_regions(prediction, applicability=None, window=9,
                      threshold_ratio=0.1):
    """Regions with far less local variation than the map as a whole."""
    filled = np.nan_to_num(prediction, nan=float(np.nanmean(prediction)))
    mean = uniform_filter(filled, size=window)
    squared = uniform_filter(filled ** 2, size=window)
    local_std = np.sqrt(np.maximum(squared - mean ** 2, 0))

    overall = float(np.nanstd(prediction))
    flat = (local_std < overall * threshold_ratio) & np.isfinite(prediction)

    print(f"  overall sd {overall:.4f}, "
          f"{flat.mean():.1%} of cells have local sd below "
          f"{threshold_ratio:.0%} of it")

    if applicability is not None:
        outside = ~applicability & np.isfinite(prediction)
        overlap = (flat & outside).sum() / max(flat.sum(), 1)
        print(f"  {overlap:.1%} of flat cells are outside the area of "
              "applicability")
        if overlap > 0.5:
            print("  -> the flatness is extrapolation, not real homogeneity")
    return flat

The overlap with the applicability mask is what distinguishes "the model is saturating" from "this area really is uniform". A lake genuinely is flat; a saturated forest prediction is not.

Explanation

Why a model cannot invent structure

Every prediction is a function of the covariates at that cell. If the output has stripes, some input has stripes at the same spacing and orientation.

That makes artefacts diagnostic. Their geometry points at the input responsible: flight-line stripes from lidar, scene-boundary steps from a mosaic, block edges from processing.

It also means the model has learned the artefact as signal. Its cross-validation score may be improved by it, because the artefact is consistent between training and test folds.

Why block edges appear only with neighbourhood features

Cell-wise features — elevation, reflectance, a distance raster — depend only on the cell. Splitting the raster into blocks changes nothing.

Neighbourhood features depend on surrounding cells. At a block boundary those cells are outside the block, so the filter sees a truncated window and returns a different value.

The result is a discontinuity at every block edge, whose spacing matches the processing block size exactly. That spacing is the identification.

Why extrapolation looks like homogeneity

A tree ensemble averages training targets in a leaf. Beyond the training feature space, every cell falls into the same extreme leaves, so every cell gets the same answer.

The output is smooth, within range, and free of the noise a real prediction has. It does not look like an error, which is why it survives review.

Comparing the flat regions against an applicability mask separates the two possibilities. Real homogeneity is inside the mask; saturation is outside it.

Why holes are an input problem

The valid mask is the intersection of every covariate's valid area. A single layer with 20% nodata caps the prediction at 80% coverage.

Filling that layer is a modelling choice with consequences: an interpolated covariate produces a prediction that inherits the interpolation's error, and nothing in the output records that.

The alternative — a second model without the problematic layer, used only where the first cannot predict — is more work and more honest, provided the output records which model produced each pixel.

Flat prediction regions compared against the applicability mask to separate real homogeneity from model saturation.
Real homogeneity is inside the mask. Saturation is outside it, and looks identical.

Edge cases or notes

  • Plot every covariate before modelling. Artefacts in the output are artefacts in the input.
  • Buffer blocks by half the largest neighbourhood, and trim before writing.
  • Flat regions plus low applicability means extrapolation, not homogeneity.
  • Holes are the union of the covariates' holes.
  • Verify the round trip on sampled cells; a reordering looks structured.
  • A striped covariate can improve the CV score while making the map wrong.
  • Record which model produced each pixel if you use a fallback model.
  • Build overviews so artefacts are visible at a glance in a viewer.

FAQ

Why does my prediction raster have stripes?

Because a covariate does. A model cannot invent structure, so stripes in the output come from flight lines, scene boundaries or mosaic seams in an input.

Why are there hard rectangular edges?

Block-wise processing with neighbourhood features and no buffer. Read each block with a buffer of half the window and trim before writing.

Why is one region completely flat?

The model is extrapolating and saturating. Check it against the applicability mask — flat regions outside the mask are saturation, not homogeneity.

Why does my prediction have holes?

One covariate has nodata there. The valid area is the intersection of all covariates, so a single gappy layer caps the coverage.

My map looks noisy and structured but wrong — what happened?

Probably a feature-order mismatch or a bad reshape. Predict a sample of cells directly and compare against the raster.

Can an artefact improve my cross-validation score?

Yes. A striped covariate is consistent between folds, so the model can exploit it and score well while producing a wrong map.

Should I fill the gappy covariate?

Only knowingly. The prediction then inherits the interpolation's error, and nothing in the output records it. A fallback model with a provenance band is more honest.