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: stripes, hard rectangular edges, a flat region, or holes where there should be data.

Five distinct causes, each with a characteristic appearance:

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

Every artefact on this page is produced for real: a random forest predicting Sentinel-2 NDVI over a 10 km window on Snowdon from Copernicus terrain, with Landsat 7 land surface temperature from August 2022 as a sixth covariate. Landsat 7's scan-line corrector failed in 2003, so every scene since carries diagonal stripes of missing data — the covariate that breaks the map.

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:12} {1 - finite.mean():6.2%} nodata  "
          f"range {np.nanmin(layer):8.2f}..{np.nanmax(layer):8.2f}")

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

  elevation     0.00% nodata  range    53.68.. 1073.06
  slope         0.00% nodata  range     0.00..   67.00
  northness     0.00% nodata  range    -1.00..    1.00
  eastness      0.00% nodata  range    -1.00..    1.00
  tpi           0.00% nodata  range   -61.76..   76.79
  lst          23.73% nodata  range    15.78..   41.07
  cells with every feature: 76.3%

The prediction's holes are the union of the covariates' holes. The surface temperature layer is 23.7% nodata, so the map is at most 76.3% complete no matter what the model does — and because the gaps are scan lines, the holes come out as stripes.

Stripes, blocks, flat regions, holes, scrambles

1. Stripes: look at the covariates

A model cannot invent structure. Stripes in the output are stripes in an input — a satellite's scan-line gaps, a lidar flight line, a scene boundary, or a mosaic seam.

Plot every covariate before modelling. It is tempting to automate this with a row-and-column test, but on this stack a simple one cannot tell the striped layer from the terrain:

for name, layer in layers.items():
    row_step = np.nanstd(np.diff(np.nanmean(layer, axis=1))) / np.nanstd(layer)
    col_step = np.nanstd(np.diff(np.nanmean(layer, axis=0))) / np.nanstd(layer)
    print(f"  {name:12} row-to-row {row_step:.3f}  column-to-column {col_step:.3f}")

  elevation    row-to-row 0.009  column-to-column 0.009
  slope        row-to-row 0.024  column-to-column 0.026
  northness    row-to-row 0.034  column-to-column 0.015
  eastness     row-to-row 0.014  column-to-column 0.035
  tpi          row-to-row 0.043  column-to-column 0.038
  lst          row-to-row 0.041  column-to-column 0.035

The temperature layer's scan-line gaps run diagonally across the grid, so they barely move a row or column mean, and the topographic position index scores higher than it does. Look at the picture.

The fix is upstream: mosaic the covariate properly, use a sensor without the defect, 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.

from scipy import ndimage

def tpi_by_blocks(elevation, block, buffer, window=15):
    """Topographic position index computed one block at a time."""
    out = np.full_like(elevation, np.nan)
    height, width = elevation.shape
    for row in range(0, height, block):
        for col in range(0, width, block):
            h, w = min(block, height - row), min(block, width - col)
            r0, c0 = max(row - buffer, 0), max(col - buffer, 0)
            sub = elevation[r0:row + h + buffer, c0:col + w + buffer]
            tpi_sub = sub - ndimage.uniform_filter(sub, size=window)
            out[row:row + h, col:col + w] = tpi_sub[row - r0:row - r0 + h, col - c0:col - c0 + w]
    return out

for buffer in (0, 15 // 2 + 1):
    diff = np.abs(tpi_by_blocks(elev, 128, buffer) - tpi)
    print(f"  buffer {buffer}: {(diff > 1e-3).mean():.1%} of cells differ from the whole-raster TPI, "
          f"by up to {diff.max():.1f} m")

  buffer 0: 17.6% of cells differ from the whole-raster TPI, by up to 73.8 m
  buffer 8: 0.0% of cells differ from the whole-raster TPI, by up to 0.0 m

Read buffered, compute, write the unbuffered part. A buffer of half the window plus one makes the blocks identical to the whole-raster result. The block size in the artefact matches the processing block, which is the tell.

3. Flat regions: where the model stops varying

A random forest predicts the mean of training targets in a leaf. Where many different cells fall into the same leaves, the output converges on a constant.

The usual cause is extrapolation — covariates outside anything the model saw. On Snowdon, three-quarters of the flat cells in the eastern half are outside the area of applicability (Example 3): high ground, on average 260 m above the rest of the east, where the combination of terrain values is unlike anything in the training sample. The observed NDVI there is anything but flat; it varies more from cell to cell than the eastern half does on average.

The other quarter are inside the mask. Flatness can also mean the covariates barely change across a region, so the model has nothing to vary on. Either way the region is smooth, in range and uninformative; the applicability mask tells you which of the two you have.

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.

covered = np.isfinite(prediction_lst).mean()
print(f"  prediction with LST: {covered:.1%} of cells; without it: {np.isfinite(prediction).mean():.1%}")

  prediction with LST: 76.3% of cells; without it: 100.0%

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. Filling the gaps by interpolation is the fourth option, and it has a cost:

filled = ~np.isfinite(lst) & np.isfinite(lst_filled)
error = np.abs(prediction_filled - ndvi)
east_ok = east & target_ok
print(f"  MAE where LST was measured: {np.nanmean(error[~filled & east_ok]):.3f}")
print(f"  MAE where LST was filled:   {np.nanmean(error[filled & east_ok]):.3f}")

  MAE where LST was measured: 0.091
  MAE where LST was filled:   0.100

The map is complete, and worse where the covariate was invented — and nothing in the output says which cells those are.

5. A scrambled map: check the round trip

X_grid = np.stack([elev, slope, np.cos(aspect), np.sin(aspect), tpi]).reshape(5, -1).T
correct = model.predict(X_grid).reshape(elev.shape)
swapped = model.predict(X_grid[:, [1, 0, 2, 3, 4]]).reshape(elev.shape)   # slope and elevation swapped

ok = target_ok.ravel()
print(f"  correlation with observed NDVI: correct {np.corrcoef(correct.ravel()[ok], ndvi.ravel()[ok])[0, 1]:.3f}, "
      f"swapped {np.corrcoef(swapped.ravel()[ok], ndvi.ravel()[ok])[0, 1]:.3f}")

rows, cols = np.nonzero(target_ok)
sample = np.random.default_rng(1).choice(len(rows), 500, replace=False)
direct = model.predict(X_grid.reshape(elev.shape + (5,))[rows[sample], cols[sample]])
print("  round trip, correct map:", np.allclose(direct, correct[rows[sample], cols[sample]], atol=1e-4))
print("  round trip, swapped map:", np.allclose(direct, swapped[rows[sample], cols[sample]], atol=1e-4))

  correlation with observed NDVI: correct 0.489, swapped 0.003
  round trip, correct map: True
  round trip, swapped map: False

A reordered feature stack produces a map that is spatially structured — because the inputs are — and has almost nothing to do with the target. Only a direct comparison catches it.

A covariate audit and buffered block prediction

The data used on this page

A 10 km window on Snowdon, 4.16–4.00° W, 53.03–53.12° N, on the 20 m grid of the least cloudy Sentinel-2 L2A scene of June–September 2025 from Earth Search. The target is that scene's NDVI; the terrain covariates come from the Copernicus DEM GLO-30 — elevation, slope, northness, eastness and a topographic position index over 300 m. The sixth covariate is land surface temperature from the cloud-free Landsat 7 scene of 10 August 2022 (Collection 2 Level-2, from the Planetary Computer), with its scan-line gaps and flagged cells left as NaN, and lst_filled is the same layer gap-filled with GDAL's inverse-distance fillnodata. model is a random forest on the five terrain features and model_lst one on all six, both trained on 3,000 cells of the western half. The terrain layers are also written to GeoTIFFs for Example 2.

import numpy as np
import planetary_computer
import rasterio
import rasterio.merge
from pystac_client import Client
from rasterio.enums import Resampling
from rasterio.fill import fillnodata
from rasterio.warp import reproject, transform_bounds
from scipy import ndimage
from sklearn.ensemble import RandomForestRegressor

AOI = (-4.16, 53.03, -4.00, 53.12)
items = Client.open("https://earth-search.aws.element84.com/v1").search(
    collections=["sentinel-2-l2a"], bbox=list(AOI), datetime="2025-06-01/2025-10-01",
    query={"eo:cloud_cover": {"lt": 60}}).items()
item = min(items, key=lambda i: i.properties["eo:cloud_cover"])

with rasterio.open(item.assets["scl"].href) as ds:                 # the 20 m grid
    win = rasterio.windows.from_bounds(*transform_bounds("EPSG:4326", ds.crs, *AOI),
                                       transform=ds.transform).round_offsets().round_lengths()
    scl, transform, crs = ds.read(1, window=win), ds.window_transform(win), ds.crs
    bounds = rasterio.windows.bounds(win, ds.transform)
shape = scl.shape
band = {}
for name in ("red", "nir"):
    with rasterio.open(item.assets[name].href) as ds:
        band[name] = ds.read(1, window=rasterio.windows.from_bounds(*bounds, transform=ds.transform),
                             out_shape=shape, resampling=Resampling.average).astype("float32") / 1e4
ndvi = (band["nir"] - band["red"]) / (band["nir"] + band["red"])

tiles = [rasterio.open(f"https://copernicus-dem-30m.s3.amazonaws.com/Copernicus_DSM_COG_10_N53_00_{t}_00_DEM/"
                       f"Copernicus_DSM_COG_10_N53_00_{t}_00_DEM.tif") for t in ("W005", "W004")]
mosaic, mosaic_transform = rasterio.merge.merge(tiles, bounds=(-4.18, 53.02, -3.98, 53.13))
elev = np.full(shape, np.nan, dtype="float32")
reproject(mosaic[0], elev, src_transform=mosaic_transform, src_crs=tiles[0].crs,
          dst_transform=transform, dst_crs=crs, resampling=Resampling.bilinear)
dz_dy, dz_dx = np.gradient(elev, -20, 20)
slope = np.degrees(np.arctan(np.hypot(dz_dx, dz_dy)))
aspect = np.arctan2(-dz_dx, -dz_dy)
tpi = elev - ndimage.uniform_filter(elev, size=15)                  # a 300 m neighbourhood

# Landsat 7 surface temperature, SLC-off, onto the same grid
pc = Client.open("https://planetarycomputer.microsoft.com/api/stac/v1", modifier=planetary_computer.sign_inplace)
scene = pc.get_collection("landsat-c2-l2").get_item("LE07_L2SP_204023_20220810_02_T1")
st, qa = np.zeros(shape, "float32"), np.ones(shape, "uint16")
with rasterio.open(scene.assets["lwir"].href) as ds:
    reproject(rasterio.band(ds, 1), st, dst_transform=transform, dst_crs=crs,
              resampling=Resampling.nearest, src_nodata=0, dst_nodata=0)
with rasterio.open(scene.assets["qa_pixel"].href) as ds:
    reproject(rasterio.band(ds, 1), qa, dst_transform=transform, dst_crs=crs,
              resampling=Resampling.nearest, src_nodata=1, dst_nodata=1)
flagged = (qa & 0b11001) > 0                                        # fill, cloud, cloud shadow
lst = np.where((st > 0) & ~flagged, st * 0.00341802 + 149.0 - 273.15, np.nan).astype("float32")
lst_filled = fillnodata(np.nan_to_num(lst), mask=np.isfinite(lst).astype("uint8"), max_search_distance=20)
lst_filled[lst_filled == 0] = np.nan

layers = {"elevation": elev, "slope": slope, "northness": np.cos(aspect),
          "eastness": np.sin(aspect), "tpi": tpi, "lst": lst}
target_ok = np.isin(scl, (2, 4, 5, 7)) & np.isfinite(ndvi)
east = np.tile(np.arange(shape[1]) >= shape[1] // 2, (shape[0], 1))

grid6 = np.stack(list(layers.values())).reshape(6, -1).T
grid6_filled = grid6.copy()
grid6_filled[:, 5] = lst_filled.ravel()
usable = (target_ok & ~east).ravel() & np.isfinite(grid6).all(axis=1)
pick = np.random.default_rng(0).choice(np.flatnonzero(usable), 3000, replace=False)
model = RandomForestRegressor(n_estimators=200, min_samples_leaf=2, random_state=0, n_jobs=4)
model.fit(grid6[pick, :5], ndvi.ravel()[pick])
model_lst = RandomForestRegressor(n_estimators=200, min_samples_leaf=2, random_state=0, n_jobs=4)
model_lst.fit(grid6[pick], ndvi.ravel()[pick])

def predict_grid(fitted, features):
    out = np.full(len(features), np.nan, dtype="float32")
    ok = np.isfinite(features).all(axis=1)
    out[ok] = fitted.predict(features[ok])
    return out.reshape(shape)

prediction = predict_grid(model, grid6[:, :5])
prediction_lst = predict_grid(model_lst, grid6)
prediction_filled = predict_grid(model_lst, grid6_filled)

profile = {"driver": "GTiff", "height": shape[0], "width": shape[1], "count": 1, "dtype": "float32",
           "crs": crs, "transform": transform, "nodata": np.nan}
layer_paths = {}
for name in ("elevation", "slope", "northness", "eastness"):
    layer_paths[name] = f"{name}.tif"
    with rasterio.open(layer_paths[name], "w", **profile) as dst:
        dst.write(layers[name].astype("float32"), 1)

print(f"{item.id}: {shape[1]} x {shape[0]} cells of 20 m; LST from {scene.id}; "
      f"{len(pick):,} training cells in the west")

S2C_30UVD_20250712_0_L2A: 543 x 509 cells of 20 m; LST from LE07_L2SP_204023_20220810_02_T1; 3,000 training cells in the west

Example 1 — auditing the covariates before modelling

import numpy as np


def covariate_audit(layers, max_nodata=0.05):
    """Nodata and ranges for every input layer, and which layer causes the holes."""
    problems = []
    valid_stack = []

    for name, layer in layers.items():
        finite = np.isfinite(layer)
        valid_stack.append(finite)
        nodata = 1 - finite.mean()
        if nodata > max_nodata:
            problems.append(f"{name}: {nodata:.1%} nodata")
        print(f"  {name:12} {nodata:6.2%} nodata  {np.nanmin(layer):9.2f}..{np.nanmax(layer):9.2f}")

    valid = np.stack(valid_stack).all(axis=0)
    print(f"  cells with every feature: {valid.mean():.1%}")
    for name, finite in zip(layers, valid_stack):
        others = np.stack([f for n, f in zip(layers, valid_stack) if n != name]).all(axis=0)
        alone = others & ~finite
        if alone.mean() > 0.01:
            print(f"  ! {name} alone removes {alone.mean():.1%} of cells")
    return valid, problems


valid, problems = covariate_audit(layers)
print(problems)

  elevation     0.00% nodata      53.68..  1073.06
  slope         0.00% nodata       0.00..    67.00
  northness     0.00% nodata      -1.00..     1.00
  eastness      0.00% nodata      -1.00..     1.00
  tpi           0.00% nodata     -61.76..    76.79
  lst          23.73% nodata      15.78..    41.07
  cells with every feature: 76.3%
  ! lst alone removes 23.7% of cells
['lst: 23.7% nodata']

The last loop is the useful one: it names the single layer responsible for the holes, which is the layer to fix or drop. It compares each layer with all the others, so a hole that two layers share is not blamed on either.

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=128, neighbourhood=15, buffer=None):
    """Block-wise prediction that is correct for neighbourhood features."""
    if buffer is None:
        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])
                    # the neighbourhood feature, computed on the buffered block
                    tpi_block = stack[0] - uniform_filter(stack[0], size=neighbourhood)
                    stack = np.concatenate([stack, tpi_block[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()
    return out_path


names = ["elevation", "slope", "northness", "eastness"]
for buffer in (0, None):
    path = predict_blocks_buffered(model, layer_paths, names, f"prediction_buffer_{buffer}.tif", buffer=buffer)
    with rasterio.open(path) as src:
        blocks = src.read(1)
    diff = np.abs(blocks - prediction)
    step = np.abs(np.diff(blocks, axis=1))
    at_edges = np.zeros(step.shape[1], bool)
    at_edges[127::128] = True
    print(f"  buffer {buffer if buffer is not None else 8}: {(diff > 0.01).mean():.1%} of cells off by more than 0.01 "
          f"(max {diff.max():.3f}); mean step across block edges {step[:, at_edges].mean():.3f}, "
          f"elsewhere {step[:, ~at_edges].mean():.3f}")

  buffer 0: 6.4% of cells off by more than 0.01 (max 0.239); mean step across block edges 0.044, elsewhere 0.023
  buffer 8: 0.0% of cells off by more than 0.01 (max 0.000); mean step across block edges 0.022, elsewhere 0.023

Trimming the buffer before writing is what removes the block edges. Without the buffer, the step in predicted NDVI across a block boundary is about twice the step between any other pair of neighbouring cells — the seams you see on the map.

Example 3 — locating the flat regions

import numpy as np
from scipy.ndimage import uniform_filter
from scipy.spatial import cKDTree
from sklearn.preprocessing import StandardScaler


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 {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("  -> mostly extrapolation, not real homogeneity")
    return flat


# a simple applicability mask: distance to the nearest training cell in scaled feature space
X_grid = grid6[:, :5]
scaler = StandardScaler().fit(X_grid[pick])
tree = cKDTree(scaler.transform(X_grid[pick]))
threshold = np.quantile(tree.query(scaler.transform(X_grid[pick]), k=2)[0][:, 1], 0.95)
applicable = (tree.query(scaler.transform(X_grid), k=1)[0] <= threshold).reshape(shape)

east_prediction = np.where(east, prediction, np.nan)
flat = find_flat_regions(east_prediction, applicability=applicable)
print(f"  {flat[east].mean():.1%} of the eastern half is flat")

nd = np.nan_to_num(ndvi)
real_sd = np.sqrt(np.maximum(uniform_filter(nd ** 2, 9) - uniform_filter(nd, 9) ** 2, 0))
print(f"  observed NDVI local sd: {real_sd[flat].mean():.3f} in the flat cells, {real_sd[east].mean():.3f} across the east")
print(f"  flat cells: mean slope {slope[flat].mean():.1f}° against {slope[east].mean():.1f}°, "
      f"mean elevation {elev[flat].mean():.0f} m against {elev[east].mean():.0f} m")

  overall sd 0.1011, 0.5% of cells have local sd below 10% of it
  77.4% of flat cells are outside the area of applicability
  -> mostly extrapolation, not real homogeneity
  1.0% of the eastern half is flat
  observed NDVI local sd: 0.095 in the flat cells, 0.075 across the east
  flat cells: mean slope 15.8° against 21.4°, mean elevation 707 m against 446 m

The overlap with the applicability mask is what distinguishes "the model is extrapolating" from "the model has nothing to vary on". Here most of the flat area is outside the mask, on high ground whose terrain the western training sample does not cover, so the fix is training data from terrain like it — not a different model. The flat cells inside the mask are the other case, where a covariate that varies there — soil, land cover — is what is missing.

A model cannot invent structure

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: scan-line gaps from Landsat 7, flight-line stripes from lidar, scene-boundary steps from a mosaic, block edges from processing.

It also means a model can learn an 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 a flat region does not look like an error

A tree ensemble averages training targets in a leaf. Every cell whose features land in the same leaves gets the same answer — whether that is because the features are outside the training range, or because they barely change across the region.

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. Flatness outside the mask is extrapolation; flatness inside it means the covariates cannot see what is changing there.

Why holes are an input problem

The valid mask is the intersection of every covariate's valid area. A single layer with 23.7% nodata caps the prediction at 76.3% 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 — on Snowdon the error rose where the temperature was filled.

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.

Block buffers, flat regions and the union of holes

• Plot every covariate before modelling. Artefacts in the output are artefacts in the input, and a row-and-column stripe test missed diagonal scan-line gaps here.
• Buffer blocks by half the largest neighbourhood plus one, and trim before writing.
• Flat regions outside the applicability mask are extrapolation; inside it, the covariates are too uniform to explain the target there. On Snowdon about three-quarters were the first kind.
• 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, and which cells were gap-filled if you fill.
• Build overviews so artefacts are visible at a glance in a viewer.

Internal links

• How to turn model predictions back into a raster — the reshape and mask ordering
• Extrapolation in space explained: the area of applicability — identifying the flat regions
• How to build spatial features for a machine learning model — where the neighbourhood features come from
• How to evaluate a spatial model honestly — why metrics missed this
• My model scores 0.95 in testing and fails in the field — the same failure in numbers
• How to clip and tile a point cloud in Python — the buffering pattern
• Raster and vector do not line up in Python — alignment artefacts
• My composite has holes, stripes or grey patches — the same symptoms in imagery

FAQ

Why does my prediction raster have stripes?

Because a covariate does. A model cannot invent structure, so stripes in the output come from scan-line gaps, flight lines, scene boundaries or mosaic seams in an input — Landsat 7's since 2003 being the best-known case.

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 plus one cell and trim before writing.

Why is one region completely flat?

Many cells are landing in the same leaves. Check it against the applicability mask: flat regions outside it are extrapolation; flat regions inside it mean the covariates barely change there, so the model has nothing to vary on.

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 — on this page the error was higher where Landsat 7's gaps were filled — and nothing in the output records it. A fallback model with a provenance band is more honest.