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.
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
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.
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.
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.
Float NoData, undeclared NoData and nansum
- Float NoData of β3.4e38 (the
float32minimum) behaves like-32768, only more so.masked=Truehandles it if the file declares it. - A file with no declared NoData gives
masked=Truenothing to mask. Check the minimum, and set the value withnp.where(band == -9999, np.nan, band)yourself. np.nansumover 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.asarrayand 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,
float320/0 gavenanwith a RuntimeWarning.
Internal links
- How to run a weighted site suitability analysis in Python β the pipeline these fixes belong in
- Site suitability explained: constraints, factors and weights β why factors must be standardised
- The raster data model explained: bands, dtype, NoData and the transform β where NoData and dtype live
- Rasterio returns the wrong values: NoData, scaling and dtype fixes β the read-side version of these problems
- Raster and vector do not line up in Python β the grid mismatch behind the NaN strips
- Slope values are wrong or absurdly steep β when the slope factor itself is broken
- How to rasterize a vector layer in Python β building constraint masks on the right grid
- My prediction raster is striped, blocky or full of NoData β the same symptoms from a model output
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.