🔧 Fix Location & accessibility Raster & rasterio Location & Accessibility · step 20 of 20 Updated

Fixing a Suitability Map That Comes Out All One Value or All NoData

Problem statement

The weighted overlay ran without an error, and the map is useless. One of two things happened:

>>> print(score.min(), score.max())
nan nan

or, just as common and harder to notice:

>>> np.unique(np.round(score[in_county], 4)).size
18

Eighteen distinct values over 1.78 million cells is a flat map with a little noise on it. Both symptoms were reproduced on a real 30 m grid of Chittenden County, Vermont — 2,220 × 1,611 cells, a Copernicus DSM and a USGS 3DEP bare-earth DTM, 2020 census blocks — and every one of them came from a single line that looked correct.

Seven causes account for nearly all broken suitability maps: an integer NoData value read as data, NaN reaching min() or max(), NaN travelling through the weighted sum, factors left in their own units, an unsigned integer that wraps, reclassification breaks in the wrong units, and a constraint mask applied the wrong way round. Only one of them raises an exception.

Quick answer

Read with a mask, normalise with the NaN-aware functions, work in float32, and check the range of every layer before combining them:

import numpy as np
import rasterio

def read_factor(path):
    with rasterio.open(path) as src:
        band = src.read(1, masked=True)            # NoData becomes a mask, not -32768
    return band.astype("float32").filled(np.nan)

def standardise(x, good, bad):
    return np.clip((bad - x) / (bad - good), 0, 1)  # explicit thresholds, NaN stays NaN

slope = read_factor("slope_int16.tif")
for name, layer in {"slope": slope}.items():
    print(name, np.nanmin(layer), np.nanmax(layer), np.isnan(layer).mean())

Measured on the county grid, reading an int16 slope layer without masked=True let -32768 into a min–max normalisation. Every cell inside the county then scored between −0.0017 and 0.0000, a standard deviation of 0.000207. The masked read gave the full 0–1 range with a standard deviation of 0.1208. Same file, one argument.

Triage table matching two symptoms, all NoData and all one value, to their causes and fixes.
Only a misaligned layer raises an exception; the other six causes produce a map.

Step-by-step solution

1. Print the range of every layer before you combine anything

Most of these failures are visible in one line per input:

for name, layer in factors.items():
    finite = np.isfinite(layer)
    print(f"{name:10s} {layer.dtype}  min {np.nanmin(layer):10.3f}  "
          f"max {np.nanmax(layer):10.3f}  NaN {1 - finite.mean():.1%}")

A minimum of -32768, -9999 or -3.4e38 is a NoData value being treated as data. A layer that is 100% NaN will not normalise. A layer whose maximum is in the thousands, next to one whose maximum is 1, is unnormalised. An integer dtype on a factor that will be multiplied by weights is an overflow waiting to happen.

2. Mask integer NoData at read time 

with rasterio.open("slope_int16.tif") as src:
    raw = src.read(1)                 # min -32768, max 56
    masked = src.read(1, masked=True) # min 0,      max 56

Measured, the raw read reported a range of −32768 to 56. Min–max scaling divides by that range of 32,824, so the real slopes — 0 to 56 degrees — occupy the last 0.17% of it, and the whole county collapses onto one value.

The fix is masked=True, then .filled(np.nan) on a float copy. Do not use np.asarray() on a masked array: measured, it hands back the original -32768 underneath the mask.

3. Use np.nanmin and np.nanmax, never .min() and .max()

A distance layer clipped to the county has NaN outside it. The ordinary methods propagate NaN, silently:

lo, hi = dist.min(), dist.max()           # nan, nan — no warning
f = (hi - dist) / (hi - lo)               # 3,576,420 of 3,576,420 cells NaN
lo, hi = np.nanmin(dist), np.nanmax(dist) # 0.0, 4481.0

Measured with NumPy 2.5.3: ndarray.min() on an array containing NaN returned nan with no warning at all, and the normalised layer was NaN in every one of 3,576,420 cells. The NaN-aware versions returned 0 and 4,481 m, leaving NaN only in the 1,792,247 cells outside the county, where it belongs.

4. Find where NaN enters the weighted sum

a + b is NaN wherever either input is NaN, so the output's NoData is the union of every input's NoData:

for name, layer in factors.items():
    print(name, f"{np.isnan(layer[in_county]).mean():.3%} NaN inside the study area")

Two measured ways it happens:

  • A constraint encoded as NaN instead of 0. Marking slopes over 15° as NaN rather than as excluded made 14.282% of the county NoData in the product, with no way to tell excluded cells from missing data.
  • A factor on a slightly different extent. Shifting one factor by 34 cells (1.02 km) left 4.220% of the county NaN along one edge.

np.nansum is not the fix. It returns 0 for a cell where every input is NaN — measured, 1,723,461 cells that should have been NoData came back as a perfectly valid-looking score of 0.

5. Put every factor on the same 0–1 scale

Adding slope in degrees to distance in metres produces a map of distance:

raw_score = 0.5 * slope_deg + 0.5 * distance_m
np.corrcoef(raw_score[in_county], distance_m[in_county])[0, 1]   # 0.9999

Measured, the correlation between the score and distance alone was 0.9999; with slope it was −0.19. The weights say 50/50 and the arithmetic says 100/0, because distance ranged over 4,481 units and slope over 56. Standardise each factor with explicit thresholds before weighting.

6. Work in float32, not uint8

Scaling factors to 0–100 and storing them as uint8 is a common way to save memory, and summing three of them wraps:

total = slope_u8 + dist_u8 + elev_u8          # dtype uint8, max 255
true = slope_u8.astype("int32") + dist_u8 + elev_u8   # max 297

Measured, 38.3% of the county's cells wrapped past 255, and the overlap between the top 5% of cells in the wrapped score and in the correct score was 0.000 — not one cell in common. There is no warning: NumPy wraps array arithmetic silently and only warns for scalars.

7. Check the units of reclassification breaks

Breaks written for percent slope applied to a slope in degrees put 44.6% of the county in the flattest class and 1.9% in the steepest. The same breaks on percent slope gave 31.1% and 16.8%. Not one value, but a map skewed towards "suitable" for no reason in the data.

8. Check which way round the constraint mask is

score * mask needs True where a site is allowed. Built from slope > 15 rather than slope <= 15, it zeroed 85.7% of the county and kept only the steep ground.

Two number lines showing real slope values squashed into a tiny fraction of the range when -32768 is included.
The normalisation is correct; its input range is not.

Code examples

Example 1 — a pre-flight check for every factor 

import numpy as np


def check_layers(layers, study_mask, expect=(0.0, 1.0)):
    """Print the facts that reveal a broken overlay before it is computed."""
    shapes = {name: layer.shape for name, layer in layers.items()}
    if len(set(shapes.values())) > 1:
        raise ValueError(f"layers are on different grids: {shapes}")

    problems = []
    for name, layer in layers.items():
        inside = layer[study_mask]
        nan_share = np.isnan(inside).mean() if inside.dtype.kind == "f" else 0.0
        lo, hi = np.nanmin(inside), np.nanmax(inside)
        print(f"{name:12s} {str(layer.dtype):8s} min {lo:10.4g} max {hi:10.4g} "
              f"NaN inside {nan_share:.2%}")
        if layer.dtype.kind in "iu":
            problems.append(f"{name}: integer dtype — cast to float32 first")
        if lo <= -9999:
            problems.append(f"{name}: minimum {lo} looks like NoData")
        if nan_share > 0:
            problems.append(f"{name}: {nan_share:.2%} NaN inside the study area")
        if lo < expect[0] - 1e-6 or hi > expect[1] + 1e-6:
            problems.append(f"{name}: range {lo:.4g}..{hi:.4g} is not standardised")
    for p in problems:
        print("  !", p)
    return not problems

Run on the raw inputs from this article, it returned False with three flags: the int16 dtype of the slope layer, and two unstandardised ranges, 0–56 for slope and 0–4,481 for distance. Run on the standardised layers from Example 2, it returned True. Note that the checks look inside the study mask only. A NoData value outside it, as here, is harmless until a layer's extent turns out to be smaller than the study area.

Example 2 — reading and standardising factors safely 

import numpy as np
import rasterio


def read_float(path, reference=None):
    """Read band 1 as float32 with NoData as NaN, and refuse a mismatched grid."""
    with rasterio.open(path) as src:
        if reference is not None and (src.transform != reference.transform
                                      or src.shape != reference.shape
                                      or src.crs != reference.crs):
            raise ValueError(f"{path} is not on the reference grid")
        band = src.read(1, masked=True)
    return band.astype("float32").filled(np.nan)


def decreasing(x, good, bad):
    """1 at or below `good`, 0 at or beyond `bad`, linear between; NaN stays NaN."""
    return np.clip((bad - x) / (bad - good), 0.0, 1.0).astype("float32")


def increasing(x, bad, good):
    return np.clip((x - bad) / (good - bad), 0.0, 1.0).astype("float32")

Thresholds chosen from the problem — "flat enough below 3°, unusable above 15°" — do not depend on the extremes of the data, so an outlier or a stray NoData value cannot compress the scale.

Example 3 — a weighted overlay that keeps NoData and exclusion apart 

import numpy as np


def weighted_overlay(factors, weights, allowed, study_mask):
    """Weighted sum with three distinct outcomes: a score, excluded (0), or NoData (NaN)."""
    total = sum(weights.values())
    if not np.isclose(total, 1.0):
        raise ValueError(f"weights sum to {total}, not 1")

    score = np.zeros(study_mask.shape, dtype="float32")
    missing = np.zeros(study_mask.shape, dtype=bool)
    for name, w in weights.items():
        layer = factors[name]
        missing |= np.isnan(layer)
        score += np.float32(w) * np.nan_to_num(layer, nan=0.0)

    score[~allowed] = 0.0                # excluded by a constraint
    score[missing & allowed] = np.nan    # genuinely unknown
    score[~study_mask] = np.nan          # outside the study area

    inside = study_mask
    print(f"scored {np.isfinite(score[inside]).mean():.1%}, "
          f"excluded {(~allowed[inside]).mean():.1%}, "
          f"NoData {(missing & allowed)[inside].mean():.2%}")
    return score

Keeping "excluded" and "unknown" as different values is what lets a reviewer see a 4% NaN strip along one edge instead of a map with an unexplained hole.

Explanation

Why an integer NoData value flattens the whole map

Min–max normalisation maps the smallest value to 0 and the largest to 1. It is linear, so the result is only as useful as the range it is given. With -32768 as the minimum and 56 as the maximum, the county's real slopes occupy 56 of 32,824 units, and every one of them lands within 0.0017 of the same number. The map is not wrong so much as rendered with a colour ramp stretched over a value that does not exist.

Explicit thresholds avoid the problem entirely, because they never look at the data's extremes.

Why NaN spreads without a warning

IEEE floating point defines every arithmetic operation involving NaN to return NaN, and ndarray.min() follows the same rule. NumPy warns when it creates a NaN — 0/0 raised invalid value encountered in scalar divide — but not when it passes one along. So a single NaN in a 3.6-million-cell array makes the minimum NaN, the normalisation NaN, and the map NaN, all without a message.

np.nanmin does warn, but only in the opposite case: an array that is entirely NaN gives RuntimeWarning: All-NaN slice encountered. If you see that, the layer did not overlap the study area at all.

Why unsigned integers wrap instead of failing

Integer arrays in NumPy use the machine's fixed-width arithmetic: 200 + 100 in uint8 is 44, because 300 modulo 256 is 44. Checking every element of a 3.6-million-cell operation for overflow would be slow, so array arithmetic does not check. Scalars are handled differently — NumPy 2 raises OverflowError: Python integer 300 out of bounds for uint8 for an out-of-range Python integer, and warns on scalar overflow — which is why a quick test in the console can behave while the full raster wraps.

A wrapped cell does not move slightly; it moves from the top of the scale to the bottom. That is why the top-5% sets of the wrapped and correct scores shared nothing at all.

Why misalignment is the only loud failure

Arrays of different shapes cannot be added, so a layer one column short fails immediately with ValueError: operands could not be broadcast together with shapes (2220,1611) (2220,1610). A layer with the same shape but a different transform — clipped, shifted, or resampled onto another origin — adds without complaint and produces NaN strips or a quietly offset score. Compare the transform and CRS, not only the shape.

Bar chart of the share of cells affected by five measured suitability failures, from all NaN to a 4% edge strip.
Each bar is a map that rendered, had a legend, and was wrong.

Edge cases or notes

  • Float NoData of −3.4e38 (the float32 minimum) behaves like -32768, only more so. masked=True handles it if the file declares it.
  • A file with no declared NoData gives masked=True nothing to mask. Check the minimum, and set the value with np.where(band == -9999, np.nan, band) yourself.
  • np.nansum over factors turns all-NaN cells into 0. That is a valid-looking score, not NoData.
  • Weights that do not sum to 1 do not flatten the map, but they change the scale, so a fixed "suitable above 0.7" threshold stops meaning anything.
  • Resampling a categorical layer bilinearly creates in-between classes that reclassification tables do not cover, which become NoData or 0.
  • Masked arrays lose their mask through np.asarray and some SciPy functions. Convert to float with NaN once, at read time.
  • A slope raster computed in degrees of longitude is a different failure — every cell comes out nearly flat or absurdly steep — and needs a projected DEM, not a nodata fix.
  • A constant layer (standard deviation 0) divides by zero in min–max scaling. Measured, float32 0/0 gave nan with a RuntimeWarning.

FAQ

Why is my whole suitability raster NaN?

Usually because .min() or .max() was used to normalise a layer containing NaN. Measured, that made all 3,576,420 cells NaN with no warning; np.nanmin and np.nanmax fixed it.

Why does my suitability map look like one colour?

A NoData value such as −32768 was read as data and stretched the normalisation range. Measured, it squeezed every county cell into a band 0.0017 wide. Read with masked=True.

Should I use np.nansum to combine factors?

No. It turns cells where every factor is NaN into a score of 0, which looks valid. Combine with ordinary addition and handle NoData and exclusions explicitly.

Why do my best sites change completely when I store scores as integers?

Unsigned 8-bit sums wrap past 255. Measured, 38.3% of cells wrapped and the top 5% of the wrapped score shared no cells with the correct top 5%. Use float32.

Why does my weighted map look exactly like my distance layer?

Because distance in metres has a range thousands of times larger than the other factors. Measured, the score correlated 0.9999 with distance. Standardise every factor to 0–1 first.

How do I tell excluded cells from missing data?

Give them different values: 0 for cells a constraint excludes and NaN for cells with no data. A map where both are NaN hides edge strips and misaligned layers.