Rasterio Output GeoTIFF Is Black or Empty: How to Fix It

Problem statement

The script runs, the file appears, and it is useless. Here is NDVI computed from a clear Sentinel-2 scene over Manchester (the setup is under The data used on this page), written with the profile of the band it came from:

import numpy as np
import rasterio

with rasterio.open("red.tif") as r, rasterio.open("nir.tif") as n:
    red, nir = r.read(1).astype("float32"), n.read(1).astype("float32")
    profile = r.profile                                  # inherited from a uint16 source!
ndvi = (nir - red) / (nir + red)
print(f"computed: {ndvi.min():.2f} … {ndvi.max():.2f}")

with rasterio.open("ndvi.tif", "w", **profile) as dst:
    dst.write(ndvi, 1)

computed: -0.36 … 1.00

Open it in QGIS and there is nothing there. Open it in Python and every cell is missing.

The frustrating part is that there is no error. dst.write accepted the array, GDAL wrote a valid GeoTIFF, and every piece of software that opens it agrees it is a well-formed file containing nothing useful.

There are six causes. Two are about the data, two about the dtype, and two about how the file is displayed — and that last pair means the file is sometimes fine and only looks broken.

Quick answer

Read the file you just wrote and look at the numbers, not the picture:

def inspect(path):
    with rasterio.open(path) as src:
        a = src.read(1, masked=True)
        print(f"{path}: dtype {src.dtypes[0]}, nodata {src.nodata}, "
              f"valid {100 * (1 - np.ma.getmaskarray(a).mean()):.1f}%", end="")
        if a.count():
            print(f", range {a.min():.3f} … {a.max():.3f}, {len(np.unique(a.compressed())):,} distinct")
        else:
            print(", nothing valid to show")


inspect("ndvi.tif")

ndvi.tif: dtype uint16, nodata 0.0, valid 0.0%, nothing valid to show

The most useful single line: valid and range. If those look right, the file is fine and the problem is display.

Check the array before writing, dtype, NoData

1. Check the array before writing it

Most "the file is black" problems are "the array was black". Assert on the way out:

def check_before_write(arr, name="array"):
    a = np.ma.masked_invalid(arr) if np.issubdtype(arr.dtype, np.floating) else arr
    finite = a.compressed() if np.ma.isMaskedArray(a) else a.ravel()
    if finite.size == 0:
        raise ValueError(f"{name} is entirely NoData or NaN — nothing to write")
    if finite.min() == finite.max():
        raise ValueError(f"{name} is constant at {finite.min()} — check the computation")
    print(f"{name}: dtype {arr.dtype}, shape {arr.shape}, range {finite.min():.3f} … {finite.max():.3f}, "
          f"{len(np.unique(finite)):,} distinct")
    return arr


check_before_write(ndvi, "ndvi")

ndvi: dtype float32, shape (898, 787), range -0.364 … 0.999, 552,331 distinct

The array was fine. That is the point of checking here: it moves the search from "what did I compute" to "what happened on the way into the file".

2. dtype truncation — the commonest real cause

print("profile dtype:", profile["dtype"], "| profile nodata:", profile["nodata"])
with rasterio.open("ndvi.tif") as src:
    stored = src.read(1)
print("what the file holds:", np.unique(stored))

profile dtype: uint16 | profile nodata: 0.0
what the file holds: [0]

Every NDVI value between −0.36 and 1.0 became 0 on its way into a uint16 file — negatives clamped to zero, positives truncated towards it. And because the profile was copied from a Sentinel-2 band whose NoData value is 0, every one of those zeros is now declared missing. The file is a valid uint16 raster in which nothing is valid.

The fix is to update the profile to match the data, not the source:

profile = profile.copy()                        # copy, then update
profile.update(dtype="float32", nodata=np.nan, count=1)
with rasterio.open("ndvi.tif", "w", **profile) as dst:
    dst.write(ndvi.astype("float32"), 1)
inspect("ndvi.tif")

ndvi.tif: dtype float32, nodata nan, valid 100.0%, range -0.364 … 0.999, 552,331 distinct

Or scale deliberately if an integer dtype is required:

scaled = profile | {"dtype": "int16", "nodata": -32768}
with rasterio.open("ndvi_int16.tif", "w", **scaled) as dst:
    dst.write(np.round(ndvi * 10_000).astype("int16"), 1)
inspect("ndvi_int16.tif")

ndvi_int16.tif: dtype int16, nodata -32768.0, valid 100.0%, range -3637.000 … 9992.000, 10,367 distinct

Guard against it generically:

def assert_dtype_fits(arr, dtype):
    info = np.iinfo(dtype) if np.issubdtype(np.dtype(dtype), np.integer) else np.finfo(dtype)
    lo, hi = float(np.nanmin(arr)), float(np.nanmax(arr))
    if lo < info.min or hi > info.max:
        raise ValueError(f"data {lo:.3f} … {hi:.3f} does not fit in {dtype} ({info.min} … {info.max})")
    if np.issubdtype(arr.dtype, np.floating) and np.issubdtype(np.dtype(dtype), np.integer):
        if not np.allclose(arr[np.isfinite(arr)] % 1, 0):
            raise ValueError(f"float data would be truncated by {dtype}")


try:
    assert_dtype_fits(ndvi, "uint16")
except ValueError as err:
    print(f"ValueError: {err}")

ValueError: data -0.364 … 0.999 does not fit in uint16 (0 … 65535)

3. NoData that swallows the whole raster

Setting a nodata value that is also a legitimate value marks every such cell as missing. Step 2 was the extreme case — every cell became the declared NoData — but the same thing happens quietly whenever 0 is both a real value and the declared fill: an NDVI of exactly zero, a count of zero, sea level. Pick a sentinel outside the data's range, or use NaN for float data:

profile.update(nodata=-9999)              # a sentinel outside the data range
# or, for float data:
profile.update(dtype="float32", nodata=np.nan)

The reverse also happens: an operation that produced NaN everywhere writes a legitimately empty raster. valid 0.0% cannot distinguish the two — check what produced the array. Here is the same NDVI for a scene three weeks earlier, 99.9% cloud, masked with the scene classification band:

with rasterio.open("cloudy_red.tif") as r, rasterio.open("cloudy_nir.tif") as n, \
        rasterio.open("cloudy_scl.tif") as s:
    cred, cnir, scl = r.read(1).astype("float32"), n.read(1).astype("float32"), s.read(1)
cloudy_ndvi = np.where(np.isin(scl, [4, 5, 6, 7]), (cnir - cred) / (cnir + cred), np.nan)
print(f"clear cells: {np.isin(scl, [4, 5, 6, 7]).sum():,} of {scl.size:,}")
try:
    check_before_write(cloudy_ndvi, "cloudy_ndvi")
except ValueError as err:
    print(f"ValueError: {err}")

clear cells: 0 of 706,726
ValueError: cloudy_ndvi is entirely NoData or NaN — nothing to write

Nothing was wrong with the code; the sky was the problem. A job that writes this file anyway produces something indistinguishable from a bug.

4. The array was never populated

from rasterio.windows import Window

height, width = ndvi.shape
out = np.zeros((height, width), dtype="float32")
for row in range(0, height, 256):
    window = Window(0, row, width, min(256, height - row))
    result = ndvi[row:row + 256]                       # computed…
                                                       # …and never assigned to out
print(f"range {out.min()} … {out.max()}")

out = np.full((height, width), np.nan, dtype="float32")   # initialise with something impossible
print("any window wrote data:", bool(np.isfinite(out).any()))

range 0.0 … 0.0
any window wrote data: False

The tell is range 0 … 0 with a plausible file size: every cell holds the initialisation value. Initialise with something impossible instead of zero, and an assert np.isfinite(out).any() before writing turns the mistake into an exception.

5. Band index and shape mistakes

count in the profile must match what you write: a profile with count=3 and one dst.write(arr, 1) leaves bands 2 and 3 filled with the default, and a three-band viewer renders that as a nearly black RGB image. And mask() returns a 3-D array even for a single band:

from rasterio.mask import mask
from shapely import box

with rasterio.open("ndvi.tif") as src:
    out_image, out_transform = mask(src, [box(*src.bounds).buffer(-500)], crop=True)
    one_band = src.profile | {"height": out_image.shape[1], "width": out_image.shape[2], "transform": out_transform}
print("mask() returned", out_image.shape)
with rasterio.open("clip.tif", "w", **one_band) as dst:
    dst.write(out_image)                               # 3-D into a count=1 dataset: fine
    dst.write(out_image[0], 1)                         # 2-D into band 1: fine
    try:
        dst.write(out_image, 1)                        # 3-D into band 1
    except ValueError as err:
        print(f"ValueError: {err}")

mask() returned (1, 798, 687)
ValueError: Source shape (1, 1, 798, 687) is inconsistent with given indexes 1

6. When the file is fine and it only looks wrong

Check the numbers before changing any code — if valid and range are sensible, the data is fine. In QGIS 4.2 a single-band float raster opens with a grey ramp stretched from its actual minimum to its actual maximum, so a correct NDVI shows as a grey gradient from −0.36 to 1.0: dull, but not black. When a correct file really does look black or washed out, it is almost always an undeclared fill value inside the stretch:

scaled = np.round(ndvi * 10_000).astype("int16")                     # as in step 2
border = np.ones_like(scaled, dtype=bool)
border[50:-50, 50:-50] = False                                       # a 500 m frame, as a clip leaves
framed = np.where(border, 32767, scaled)                             # filled with 32767, NoData not declared

share = (int(scaled.max()) - int(scaled.min())) / (int(framed.max()) - int(framed.min()))
print(f"stretch {framed.min()} … {framed.max()}: the real NDVI occupies the bottom {share:.0%} of the ramp")

stretch -3637 … 32767: the real NDVI occupies the bottom 37% of the ramp

With the fill included in the stretch, every real value is squeezed into the dark end of the ramp and the frame glows white around a murky picture. Declaring NoData removes the fill from the statistics; writing statistics and overviews into the file lets any viewer pick a sensible stretch without scanning the whole raster:

with rasterio.open("ndvi.tif", "r+") as dst:
    dst.nodata = np.nan
    dst.update_tags(1, STATISTICS_MINIMUM="-0.36", STATISTICS_MAXIMUM="1.0")
    dst.build_overviews([2, 4, 8], rasterio.enums.Resampling.average)

A single-band raster also gets a grey ramp rather than a colour map — for NDVI, choose a diverging one deliberately.

A safe write wrapper

The data used on this page

Sentinel-2 L2A from Element 84's Earth Search catalogue, cut to a window over Manchester (2.30–2.18°W, 53.44–53.52°N): the red and near-infrared bands of the clear scene S2B_30UWE_20260712 (0.06% cloud), saved exactly as the source delivers them — uint16, NoData 0, in UTM zone 30N — and the red, near-infrared and scene-classification bands of S2A_30UWE_20260706, a scene 99.9% covered by cloud.

import os

os.environ["GDAL_DISABLE_READDIR_ON_OPEN"] = "EMPTY_DIR"
os.environ["AWS_NO_SIGN_REQUEST"] = "YES"

import numpy as np
import pyproj
import rasterio
from pystac_client import Client
from rasterio.enums import Resampling
from rasterio.windows import from_bounds

BOX = (-2.30, 53.44, -2.18, 53.52)
catalog = Client.open("https://earth-search.aws.element84.com/v1")


def cut(item, asset, dst, shape=None):
    with rasterio.open(item.assets[asset].href) as src:
        to_utm = pyproj.Transformer.from_crs(4326, src.crs, always_xy=True)
        left, bottom = to_utm.transform(BOX[0], BOX[1])
        right, top = to_utm.transform(BOX[2], BOX[3])
        window = from_bounds(left, bottom, right, top, src.transform).round_offsets().round_lengths()
        kwargs = {"out_shape": shape, "resampling": Resampling.nearest} if shape else {}
        data = src.read(1, window=window, **kwargs)
        profile = src.profile | {"driver": "GTiff", "width": data.shape[1], "height": data.shape[0],
                                 "transform": src.window_transform(window) if not shape else
                                 rasterio.transform.from_bounds(*rasterio.windows.bounds(window, src.transform),
                                                                data.shape[1], data.shape[0])}
    with rasterio.open(dst, "w", **profile) as out:
        out.write(data, 1)
    return data.shape


clear = catalog.get_collection("sentinel-2-l2a").get_item("S2B_30UWE_20260712_0_L2A")
cloudy = catalog.get_collection("sentinel-2-l2a").get_item("S2A_30UWE_20260706_0_L2A")
shape = cut(clear, "red", "red.tif")
cut(clear, "nir", "nir.tif")
cut(cloudy, "red", "cloudy_red.tif")
cut(cloudy, "nir", "cloudy_nir.tif")
cut(cloudy, "scl", "cloudy_scl.tif", shape=shape)          # 20 m classes onto the 10 m grid
print(f"window {shape[1]} × {shape[0]} at 10 m; cloud cover {clear.properties['eo:cloud_cover']:.2f}% "
      f"and {cloudy.properties['eo:cloud_cover']:.1f}%")

window 787 × 898 at 10 m; cloud cover 0.06% and 99.9%

Example 1: a safe write wrapper

Every check above, in one function you can put between any computation and any file:

from pathlib import Path


def write_raster(arr, path, *, profile, nodata=None, descriptions=None, overviews=(2, 4, 8, 16)):
    """Write an array with the profile corrected to match it, and refuse empty output."""
    arr = np.asarray(arr)
    if arr.ndim == 2:
        arr = arr[np.newaxis, ...]
    count, height, width = arr.shape

    if np.issubdtype(arr.dtype, np.floating):
        finite = arr[np.isfinite(arr)]
        nodata = np.nan if nodata is None else nodata
    else:
        finite = arr.ravel() if nodata is None else arr[arr != nodata]

    if finite.size == 0:
        raise ValueError(f"refusing to write {path}: no valid cells")
    if finite.min() == finite.max():
        print(f"  warning: {path} is constant at {finite.min()}")

    profile = {**profile}
    profile.update(dtype=arr.dtype.name, count=count, height=height, width=width,
                   nodata=nodata, compress="deflate", tiled=True, blockxsize=256, blockysize=256)

    tmp = Path(path).with_suffix(".tmp.tif")
    with rasterio.open(tmp, "w", **profile) as dst:
        dst.write(arr)
        for i, name in enumerate(descriptions or [], start=1):
            dst.set_band_description(i, name)
        if overviews:
            dst.build_overviews(list(overviews), Resampling.average)
        dst.update_tags(1, STATISTICS_MINIMUM=str(float(finite.min())),
                        STATISTICS_MAXIMUM=str(float(finite.max())))
    tmp.replace(path)
    print(f"  wrote {path}  {arr.dtype} {arr.shape}  {finite.min():.3f} … {finite.max():.3f}")
    return path


with rasterio.open("red.tif") as src:
    source_profile = src.profile                        # the uint16 profile that caused the trouble
write_raster(ndvi, "ndvi_safe.tif", profile=source_profile, descriptions=["ndvi"])
try:
    write_raster(cloudy_ndvi, "ndvi_cloudy.tif", profile=source_profile)
except ValueError as err:
    print(f"  ValueError: {err}")

  wrote ndvi_safe.tif  float32 (1, 898, 787)  -0.364 … 0.999
  ValueError: refusing to write ndvi_cloudy.tif: no valid cells

Four properties make this hard to misuse. The dtype and shape come from the array, so a stale profile cannot truncate the data — handed the same uint16 profile that broke step 2, it writes float32. An all-NoData result raises instead of being written. The file is written to a temporary path and renamed on success, so a crash never leaves a plausible-looking black file behind. And the printed range is the confirmation you would otherwise go looking for.

Example 2: finding which black file broke a batch run

After a folder-scale job, check the outputs rather than trusting the exit code. This folder holds the three NDVI files this page has made — the good one, the one written with the uint16 profile, and the cloudy scene written without a check:

import shutil

Path("out/ndvi").mkdir(parents=True, exist_ok=True)
shutil.copy("ndvi_safe.tif", "out/ndvi/2026-07-12.tif")
with rasterio.open("out/ndvi/2026-07-12_uint16.tif", "w", **source_profile) as dst:
    dst.write(ndvi, 1)                                                    # step 2's mistake
with rasterio.open("out/ndvi/2026-07-06.tif", "w", **(source_profile | {"dtype": "float32", "nodata": np.nan})) as dst:
    dst.write(cloudy_ndvi.astype("float32"), 1)                           # step 3's sky


def audit_outputs(folder):
    bad = []
    paths = sorted(Path(folder).glob("*.tif"))
    for p in paths:
        with rasterio.open(p) as src:
            a = src.read(1, masked=True)
            valid = 1 - np.ma.getmaskarray(a).mean()
            vals = a.compressed()
            if vals.size == 0:
                reason = "all NoData"
            elif vals.min() == vals.max():
                reason = f"constant {vals.min()}"
            elif p.stat().st_size < 10_000:
                reason = f"suspiciously small ({p.stat().st_size} B)"
            else:
                continue
            bad.append((p.name, reason, f"{100 * valid:.1f}% valid"))
    for name, reason, valid in bad:
        print(f"  ✗ {name:<24} {reason:<14} {valid}")
    print(f"{len(bad)} bad of {len(paths)}")
    return bad


bad = audit_outputs("out/ndvi/")

  ✗ 2026-07-06.tif           all NoData     0.0% valid
  ✗ 2026-07-12_uint16.tif    all NoData     0.0% valid
2 bad of 3

Both failures read "all NoData", and they have different causes — a dtype bug and a cloudy day — which is exactly why the audit must say which files and why a person then has to look at what produced each. Collapsing everything into "failed" would lose even that. The same reporting logic belongs in the job itself — see how to log and summarise errors in a batch job.

Example 3: proving the file is fine when a viewer says otherwise

import matplotlib

matplotlib.use("Agg")
import matplotlib.pyplot as plt


def show_honestly(path, band=1, percentile=(2, 98)):
    with rasterio.open(path) as src:
        a = src.read(band, masked=True)
        crs, res = src.crs, src.res
    vals = a.compressed()
    lo, hi = np.percentile(vals, percentile)

    fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(13, 5))
    im = ax1.imshow(a, vmin=lo, vmax=hi, cmap="RdYlGn")
    ax1.set_title(f"{path}\n{crs} · {res[0]} m · stretch {lo:.2f} … {hi:.2f}")
    fig.colorbar(im, ax=ax1, shrink=0.8)

    ax2.hist(vals, bins=80, color="#0ea5e9")
    ax2.axvline(lo, color="red", ls="--")
    ax2.axvline(hi, color="red", ls="--")
    ax2.set_title(f"{vals.size:,} valid cells · {vals.min():.3f} … {vals.max():.3f}")
    plt.tight_layout()
    print(f"{path}: 2–98% stretch {lo:.2f} … {hi:.2f}; full range {vals.min():.2f} … {vals.max():.2f}")
    return fig


show_honestly("ndvi_safe.tif").savefig("ndvi_honest.png", dpi=80)

ndvi_safe.tif: 2–98% stretch 0.03 … 0.87; full range -0.36 … 1.00

The histogram is the part that settles the argument. A raster that shows a healthy spread of values is a display problem at most; one whose histogram is a single spike is a data problem. The 2–98 percentile stretch is what QGIS's "Cumulative count cut" does, and applying it manually shows whether a better stretch was all the map needed.

A GeoTIFF has no notion of looking right

A GeoTIFF has no concept of "looks right". It stores an array, a dtype, an optional NoData tag, a transform and a CRS, and anything self-consistent is a valid file. Every failure in this article is a case where the file is valid and the meaning was lost somewhere between the computation and the disk.

The dtype cases are the most instructive because they are pure type-system failures. NumPy's casting rules are silent by design: np.float64(0.87).astype("uint16") is 0, and that is documented, intended behaviour, not a bug. GDAL then faithfully stores the zeros. Nothing in the chain has enough context to know that a value of 0.87 mattered. This is why the fix is a check rather than a setting — the information that the data was float NDVI exists only in your head and in the assertion you write down.

The NoData cases are a variant of the same problem. NoData is a value that means "not a value", and its specialness lives in a header tag rather than in the number itself. Choose a sentinel that collides with real data and every real cell becomes invisible; forget to set one and the sentinel is included in statistics and colour stretches. The raster data model covers this in depth, and the practical consequence is the same everywhere: NoData must be chosen against the actual value range, not by habit.

The viewer cases deserve their own note, because "the file is black" is usually a claim about rendering, not about data. Turning an array of numbers into pixels requires a mapping from value range to colour range, and no such mapping is stored in a plain GeoTIFF. QGIS guesses, and its guess for a single-band float raster — a linear stretch over the full possible range, including undeclared NoData — produces black for almost any real dataset. Writing statistics tags and overviews gives every viewer the information to guess well. It costs a few kilobytes and removes an entire category of false alarms.

Finally, the structural point: a raster write should be verified by a raster read. The write API cannot tell you whether the result is meaningful, because it never sees the data twice. Reading the file back and printing valid and range costs milliseconds and catches every failure described here — which is why the wrapper in Example 1 does exactly that, and why the batch audit in Example 2 is worth running even after a job that reported no errors.

Live profiles, band counts and NaN nodata

• profile = src.profile is a live view. Always .copy() before updating, or you may mutate the open dataset's profile.

• count must match what you write. Unwritten bands are filled with the NoData value or zero, and an RGB viewer renders that as black.

• nodata=np.nan only works on float dtypes. On an integer raster it is silently ineffective.

• Writing outside the with block writes nothing — the dataset must be open. An exception inside the block can also leave a header-only file.

• A tiny compressed GeoTIFF is usually a header with no image data. An exception inside the with block after rasterio.open(..., "w") still leaves a file: here a 798 × 687 float32 deflate-compressed one came out at 3,864 bytes, all NoData. Uncompressed, the same failure leaves a full-size 2.2 MB file of NoData, so size alone is not a test — read it back.

• dst.write(arr) needs a 3-D array; dst.write(arr, 1) needs a 2-D one. rasterio.mask returns 3-D even for one band.

• Overviews built with the wrong resampling can make a categorical raster look black when zoomed out while the full-resolution data is fine.

• Negative floats do not reliably wrap in an unsigned dtype. np.array([-0.31]).astype("uint16") gives 0 on this machine; converting an out-of-range float to an integer is undefined in C, so the result depends on the platform. Scale and range-check rather than relying on either behaviour.

• gdalinfo -stats out.tif prints min, max, mean and standard deviation per band and settles the data-versus-display question in one command.

• Compression never changes values. compress="deflate" is lossless; if values changed, it was the dtype.

Internal links

• The raster data model explained — dtype, NoData and why both are load-bearing
• Rasterio returns the wrong values: NoData, scaling and dtype fixes — the read-side version of the same problem
• Introduction to Rasterio — profiles, reading and writing
• Raster and vector do not line up in Python — when the file is fine but in the wrong place
• Why is my GeoPandas plot blank or empty? — the vector equivalent of a black raster
• How to log and summarise errors in a batch GIS job — catching the bad output in a batch of good ones
• How to clip a raster to a polygon in Python — where 3-D out_image shapes come from
• Raster resampling explained — why overviews can look wrong on categorical data

FAQ

My raster is all zeros. What happened?

Either the array was never populated — check for a loop that computes without assigning — or float values were written to an integer dtype and truncated. Print the array's range before writing.

Why did my float NDVI become zeros?

The profile was inherited from a uint16 source. Values between −1 and 1 truncate to 0. Copy the profile and update dtype="float32" before writing.

The values look correct in Python but QGIS shows black.

A display stretch problem, not a data problem. The viewer stretches the grey ramp between the minimum and maximum it estimates, and an undeclared fill value such as 32767 or −9999 is included in that range, pushing the real data into one end of the ramp. Declare the fill as NoData, or set Min/Max to the real range by hand.

Why is my output only a few kilobytes?

For a compressed file, that usually means a header and no image data: the write was interrupted. Look for an exception inside the with block. An uncompressed file interrupted the same way is full size, so read the file back rather than trusting its size.

Should nodata be 0?

Only when 0 is genuinely impossible in your data. On elevation, reflectance or NDVI it is a real value, and setting it as NoData marks the whole raster missing. Use a sentinel outside the range, or NaN for float data.

How do I stop this happening again?

Check the array's range and validity before writing, derive the profile's dtype from the array rather than the source, and read the file back to confirm. The wrapper in Example 1 does all three.

Bands 2 and 3 are black but band 1 is fine.

The profile declares count=3 and only band 1 was written. Set count to the number of bands you actually write.