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

On this page

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,219 Γ— 1,610 cells, a Copernicus DSM and a USGS 3DEP bare-earth DTM, 2020 census blocks and TIGER roads β€” 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())
slope 0.0 56.0 0.5005922873881413

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.9983 and 1.0000, a standard deviation of 0.000206. The masked read gave the full 0–1 range with a standard deviation of 0.1208. Same file, one argument.

Symptoms of a broken suitability overlay, their causes, and the fix for each. the symptom points at the causesymptom β€” causefixall NaN β€” .min() on a layer withNaNnp.nanmin and np.nanmaxone value β€” βˆ’32768 read as datasrc.read(1, masked=True)NaN strip β€” factor on anotherextentcompare transform, not only shapebest sites scrambled β€” uint8 wrapscast to float32 before summingmap equals one factor β€” raw unitsstandardise each factor to 0–1Only a shape mismatch raises an exception; the rest render a map with a legend spatialworkflow.io
Only a misaligned layer raises an exception; the other six causes produce a map.

Print every range, mask NoData, use NaN-aware statistics

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,572,590 of 3,572,590 cells NaN
lo, hi = np.nanmin(dist), np.nanmax(dist) # 0.0, 4781.7

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,572,590 cells. The NaN-aware versions returned 0 and 4,782 m, leaving NaN only in the 1,788,411 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.275% 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.219% 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, all 1,788,411 cells outside the county, which should have been NoData, came back as a perfectly valid-looking score of 0.

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

What happens to a weighted overlay when the factors are on different scales. the weights said 50/50 and the arithmetic said 100/0distance, metresslope, degreesrange over the countyβœ•4,782 units56 unitsthe weight you wrote0.500.50the weight it hadβœ•β‰ˆ 1.00βœ•β‰ˆ 0.00The score correlated with distance alone at 0.9999, and with slope alone at 0.04 spatialworkflow.io
Standardise with explicit thresholds before weighting, or the weights are decoration.

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.04. The weights say 50/50 and the arithmetic says 100/0, because distance ranged over 4,782 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 299

Measured, with lower ground scoring higher on elevation, 44.4% 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 β€” 5, 10, 15 and 25% β€” 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.

Why an unmasked βˆ’32768 flattens a min–max normalisation: real slopes occupy 0.17% of the range. the same int16 slope layer, normalised two waysread without a mask: range βˆ’32768 … 56βˆ’32768 (NoData)0–56Β° here: 0.17%every county cell scores 0.9983 … 1.0000read with masked=True: range 0 … 560Β°15Β°30Β°45Β°56Β°Standard deviation inside the county: 0.000206 unmasked, 0.1208 masked spatialworkflow.io
The normalisation is correct; its input range is not.

Pre-flight checks and a safe weighted overlay

The data used on this page

Chittenden County, Vermont, on a 30 m grid in NAD83 / Vermont (EPSG:32145) covering the county plus 2 km. The county comes from the TIGER/Line 2020 blocks; distance is to the nearest road in the county's TIGER 2020 road layer; slope comes from the Copernicus GLO-30 surface model and elevation from the USGS 3DEP 1 arc-second bare-earth DTM. The slope is also written out as an int16 GeoTIFF with -32768 as NoData outside the county β€” the file the quick answer reads.

import os
import urllib.request

import geopandas as gpd
import numpy as np
import rasterio
from rasterio.enums import Resampling
from rasterio.features import rasterize
from rasterio.merge import merge
from rasterio.transform import from_origin
from rasterio.warp import reproject
from scipy import ndimage

CRS, CELL = 32145, 30.0                                   # NAD83 / Vermont, 30 m cells
downloads = {
    "tl_2020_50_tabblock20.zip": "https://www2.census.gov/geo/tiger/TIGER2020/TABBLOCK20/tl_2020_50_tabblock20.zip",
    "tl_2020_50007_roads.zip": "https://www2.census.gov/geo/tiger/TIGER2020/ROADS/tl_2020_50007_roads.zip",
}
for tile in ("n45w074", "n45w073"):                     # USGS 3DEP 1 arc-second bare-earth DTM
    downloads[f"USGS_1_{tile}.tif"] = ("https://prd-tnm.s3.amazonaws.com/StagedProducts/Elevation/1/TIFF/"
                                       f"current/{tile}/USGS_1_{tile}.tif")
for tile in ("N44_00_W074_00", "N44_00_W073_00"):       # Copernicus GLO-30 DSM, 1 degree tiles
    name = f"Copernicus_DSM_COG_10_{tile}_DEM"
    downloads[name + ".tif"] = f"https://copernicus-dem-30m.s3.amazonaws.com/{name}/{name}.tif"
for name, url in downloads.items():
    if not os.path.exists(name):
        urllib.request.urlretrieve(url, name)

blocks = gpd.read_file("tl_2020_50_tabblock20.zip").to_crs(CRS)
blocks = blocks[blocks.COUNTYFP20 == "007"]                           # Chittenden County
west, south, east, north = blocks.total_bounds + np.array([-2000, -2000, 2000, 2000])   # 2 km pad
width, height = int(np.ceil((east - west) / CELL)), int(np.ceil((north - south) / CELL))
transform = from_origin(west, north, CELL, CELL)
grid = dict(out_shape=(height, width), transform=transform)

in_county = rasterize(blocks.geometry, fill=0, default_value=1, dtype="uint8", **grid).astype(bool)
water = rasterize(blocks[blocks.ALAND20 == 0].geometry, fill=0, default_value=1, dtype="uint8", **grid).astype(bool)

roads = gpd.read_file("tl_2020_50007_roads.zip").to_crs(CRS)             # every road in the county
on_road = rasterize(roads.geometry, fill=0, default_value=1, dtype="uint8", all_touched=True, **grid).astype(bool)
dist_m = (ndimage.distance_transform_edt(~on_road) * CELL).astype("float32")



def to_grid(names):
    sources = [rasterio.open(n) for n in names]
    mosaic, mosaic_transform = merge(sources)
    out = np.full((height, width), np.nan, dtype="float32")
    reproject(mosaic[0], out, src_transform=mosaic_transform, src_crs=sources[0].crs, src_nodata=sources[0].nodata,
              dst_transform=transform, dst_crs=CRS, resampling=Resampling.bilinear, dst_nodata=np.nan)
    return out


dsm = to_grid([n for n in downloads if n.endswith("_DEM.tif")])
dtm = to_grid([n for n in downloads if n.startswith("USGS_1_")])
dz_dy, dz_dx = np.gradient(dsm, CELL)
slope_deg = np.degrees(np.arctan(np.hypot(dz_dx, dz_dy)))

# a slope layer as it is often delivered: int16 whole degrees, -32768 outside the county
slope_int = np.where(in_county, np.round(slope_deg), -32768).astype("int16")
with rasterio.open("slope_int16.tif", "w", driver="GTiff", height=height, width=width, count=1, dtype="int16",
                   crs=CRS, transform=transform, nodata=-32768) as dst:
    dst.write(slope_int, 1)
print(f"grid {height:,} x {width:,} = {height * width:,} cells, {in_county.sum():,} in the county")
grid 2,219 x 1,610 = 3,572,590 cells, 1,784,179 in the county

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


slope_raw = rasterio.open("slope_int16.tif").read(1)
print(check_layers({"slope": slope_raw, "distance": dist_m}, in_county))
slope        int16    min          0 max         56 NaN inside 0.00%
distance     float32  min          0 max       4782 NaN inside 0.00%
  ! slope: integer dtype β€” cast to float32 first
  ! slope: range 0..56 is not standardised
  ! distance: range 0..4782 is not standardised
False

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,782 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")


with rasterio.open("slope_int16.tif") as ref:
    slope_f = read_float("slope_int16.tif", reference=ref)
factors = {"slope": decreasing(slope_f, good=3, bad=15), "distance": decreasing(dist_m, good=0, bad=3000)}
print(check_layers(factors, in_county))
slope        float32  min          0 max          1 NaN inside 0.00%
distance     float32  min          0 max          1 NaN inside 0.00%
True

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


score = weighted_overlay(factors, {"slope": 0.5, "distance": 0.5}, in_county & (slope_f <= 15), in_county)
scored 100.0%, excluded 13.1%, NoData 0.00%

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.

Slope classes written for percent slope, applied to slope in degrees and in percent, Chittenden County. the same breaks, 5 / 10 / 15 / 25, on two unitsslope layer inflattest classsteepest classdegrees (wrong)44.6%1.9%percent (right)31.1%16.8%share of the county in each class; a degree is a larger step than a percentNothing fails; the map just leans towards "suitable" spatialworkflow.io
Reclassification breaks carry a unit even when the table does not say so.

Example 4 β€” every symptom on this page, measured

import warnings

inside = in_county
with rasterio.open("slope_int16.tif") as src:
    raw, masked = src.read(1), src.read(1, masked=True)
print(f"raw read {raw.min()}..{raw.max()}; masked {masked.min()}..{masked.max()}; "
      f"np.asarray on the masked array: min {np.asarray(masked).min()}")
r = raw.astype("float32")                                # int16 arithmetic would overflow as well
naive = (r - r.min()) / (r.max() - r.min())
good = (masked - masked.min()) / (masked.max() - masked.min())
print(f"min-max with -32768 inside: county cells {naive[inside].min():.4f}..{naive[inside].max():.4f}, "
      f"sd {naive[inside].std():.6f}, distinct values {np.unique(np.round(naive[inside], 4)).size}; "
      f"masked sd {good[inside].std():.4f}")

dist = np.where(inside, dist_m, np.nan)
with warnings.catch_warnings(record=True) as caught:
    warnings.simplefilter("always")
    lo, hi = dist.min(), dist.max()
f = (hi - dist) / (hi - lo)
print(f"dist.min() {lo}, warnings raised: {len(caught)}; NaN after normalising: {np.isnan(f).sum():,} of {f.size:,}; "
      f"nanmin {np.nanmin(dist):.1f}, nanmax {np.nanmax(dist):,.1f}; cells outside the county {(~inside).sum():,}")

fs = decreasing(slope_deg, good=3, bad=15)
as_nan = np.where(slope_deg > 15, np.nan, fs) * decreasing(dist_m, 0, 3000)
print(f"constraint as NaN: {np.isnan(as_nan[inside]).mean():.3%} of the county NoData")
shifted = np.full_like(dist_m, np.nan)
shifted[:, 34:] = dist[:, :-34]                          # dist is NaN outside the county
print(f"one factor shifted 34 cells ({34 * 30 / 1000:.2f} km): {np.isnan(shifted[inside]).mean():.3%} of the county NaN")
both = np.nansum(np.stack([np.where(inside, fs, np.nan), dist]), axis=0)
print(f"np.nansum turns {(~inside).sum():,} all-NaN cells into 0: {(both[~inside] == 0).sum():,}")

raw_score = 0.5 * slope_deg + 0.5 * dist_m
print(f"raw 50/50 score: corr with distance {np.corrcoef(raw_score[inside], dist_m[inside])[0, 1]:.4f}, "
      f"with slope {np.corrcoef(raw_score[inside], slope_deg[inside])[0, 1]:.2f}; "
      f"ranges {np.ptp(dist_m[inside]):,.0f} m and {np.ptp(slope_deg[inside]):.0f} degrees")

to100 = lambda x: np.round(100 * x).astype("uint8")
slope_u8, dist_u8 = to100(fs), to100(decreasing(dist_m, 0, 3000))
elev_u8 = to100(decreasing(dtm, np.nanmin(dtm[inside]), np.nanmax(dtm[inside])))   # lower ground scores higher
wrapped = (slope_u8 + dist_u8 + elev_u8)[inside]
true = (slope_u8.astype("int32") + dist_u8 + elev_u8)[inside]
top_w, top_t = wrapped >= np.percentile(wrapped, 95), true >= np.percentile(true, 95)
print(f"uint8 sum: {(true > 255).mean():.1%} of cells wrapped; max true {true.max()}; "
      f"top-5% overlap (Jaccard) {(top_w & top_t).sum() / (top_w | top_t).sum():.3f}")

breaks = [5, 10, 15, 25]                                 # classes written for percent slope
percent = 100 * np.tan(np.radians(slope_deg))
for label, x in (("degrees", slope_deg), ("percent", percent)):
    cls = np.digitize(x[inside], breaks)
    print(f"percent breaks on slope in {label}: flattest {np.mean(cls == 0):.1%}, steepest {np.mean(cls == 4):.1%}")
print(f"mask built from slope > 15: zeroes {(~(slope_deg > 15))[inside].mean():.1%} of the county")
raw read -32768..56; masked 0..56; np.asarray on the masked array: min -32768
min-max with -32768 inside: county cells 0.9983..1.0000, sd 0.000206, distinct values 18; masked sd 0.1208
dist.min() nan, warnings raised: 0; NaN after normalising: 3,572,590 of 3,572,590; nanmin 0.0, nanmax 4,781.7; cells outside the county 1,788,411
constraint as NaN: 14.275% of the county NoData
one factor shifted 34 cells (1.02 km): 4.219% of the county NaN
np.nansum turns 1,788,411 all-NaN cells into 0: 1,788,411
raw 50/50 score: corr with distance 0.9999, with slope 0.04; ranges 4,782 m and 56 degrees
uint8 sum: 44.4% of cells wrapped; max true 299; top-5% overlap (Jaccard) 0.000
percent breaks on slope in degrees: flattest 44.6%, steepest 1.9%
percent breaks on slope in percent: flattest 31.1%, steepest 16.8%
mask built from slope > 15: zeroes 85.7% of the county

How one NoData value flattens the map

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.

Measured share of cells affected by five suitability-overlay failures. share of cells affected, Chittenden County 30 m grid.min() with NaN100% NaNmask inverted85.7% zeroeduint8 sum wrapped44.4% wrappedconstraint as NaN14.3% NaNlayer 1 km offset4.2% NaN strip0100%After the uint8 wrap, the top 5% of cells shared none with the correct top 5% spatialworkflow.io
Each bar is a map that rendered, had a legend, and was wrong.

Float NoData, undeclared NoData and nansum

  • 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,572,590 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, 44.4% 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.