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

Problem statement #

The script runs, the file appears, and it is useless:

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

Open it in QGIS and the layer is uniformly black. Open it in Python and the array is all zeros, or all -9999, or all nan. Sometimes the file is the right size on disk and still shows nothing; sometimes it is 4 KB when it should be 40 MB.

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 viewer stretches values — and that last pair means the file is often fine and only looks broken.

Quick answer #

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

import rasterio, numpy as np

with rasterio.open("ndvi.tif") as src:
    a = src.read(1, masked=True)
    print("dtype   ", src.dtypes[0])
    print("nodata  ", src.nodata)
    print("valid   ", f"{100 * (1 - np.ma.getmaskarray(a).mean()):.1f}%")
    print("range   ", a.min(), "…", a.max())
    print("unique  ", len(np.unique(a.compressed()))) 

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

Step-by-step solution #

1. Check the array before writing it #

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

import numpy as np

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()
    print(f"{name}: dtype {arr.dtype}, shape {arr.shape}, "
          f"range {finite.min()} … {finite.max()}, "
          f"{len(np.unique(finite))} distinct")
    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")
    return arr

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

Two lines of assertion turn a silent, invisible failure into a loud one at the point where the cause is still on screen.

2. dtype truncation — the commonest real cause #

ndvi = (nir - red) / (nir + red)        # float, −1.0 … 1.0
print(ndvi.min(), ndvi.max())           # -0.31 0.87

profile = src.profile                   # inherited from a uint16 source!
print(profile["dtype"])                 # 'uint16'

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

Every value between −1 and 1 becomes 0 or 65535. The file is a valid uint16 raster of almost entirely zeros.

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

profile = src.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)

Or scale deliberately if an integer dtype is required:

profile.update(dtype="int16", nodata=-32768)
dst.write(np.round(ndvi * 10_000).astype("int16"), 1)      # −3100 … 8700

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} … {hi} does not fit in {dtype} "
                         f"({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}")

assert_dtype_fits(ndvi, "uint16")     # ValueError: float data would be truncated by uint16

3. NoData that swallows the whole raster #

print(src.nodata)          # 0.0
print(a.min(), a.max())    # 0 0

Setting nodata=0 on a raster where 0 is a legitimate value marks every valid cell as missing. Viewers then render the whole thing transparent or black, and masked=True reads produce an entirely masked array.

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 (division by zero, a mask applied twice) writes a legitimately empty raster. valid 0.0% cannot distinguish the two — check what produced the array.

denom = nir + red
print(f"{(denom == 0).sum():,} cells with zero denominator")   # 3,048,192  ← the whole grid

4. The array was never populated #

out = np.zeros((height, width), dtype="float32")
for window in windows:
    result = process(src.read(window=window))    # computed…
                                                  # …and never assigned
with rasterio.open(dst_path, "w", **profile) as dst:
    dst.write(out, 1)                             # writes the zeros

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 the mistake becomes visible immediately:

out = np.full((height, width), np.nan, dtype="float32")
...
assert np.isfinite(out).any(), "no window wrote any data"

5. Band index and shape mistakes #

dst.write(arr)          # arr is 2-D → ValueError, or writes to band 1 only
dst.write(arr, 1)       # 2-D array into band 1 — correct
dst.write(arr3d)        # 3-D array, shape (count, h, w) — correct

Two related traps. 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 (0 or the NoData value), and a three-band viewer renders that as a nearly black RGB image.

profile.update(count=1)      # writing one band? say so

And mask() returns a 3-D array even for a single band:

out_image, out_transform = mask(src, shapes, crop=True)
print(out_image.shape)          # (1, 1264, 986)
dst.write(out_image)            # correct — 3-D into a count=1 dataset
dst.write(out_image[0], 1)      # also correct
dst.write(out_image, 1)         # ValueError: Source shape is inconsistent

6. The file is fine and the viewer is wrong #

This is more than half of all reports, and it is worth checking before changing any code.

with rasterio.open("ndvi.tif") as src:
    a = src.read(1, masked=True)
print(a.min(), a.max(), a.mean())      # -0.31 0.87 0.42   ← the data is fine

Why it looks black:

• QGIS defaults to a 0–255 stretch on float data. Values of 0.42 render as almost black. Right-click the layer → Properties → Symbology → set Min/Max to the actual range, or choose "Cumulative count cut".
• A single-band float raster gets a grey ramp, not a colour map. Set one deliberately.
• NoData is not declared, so -9999 is included in the stretch and compresses all real values into the bottom pixel of the range.

Declare the range in the file so every viewer gets it right:

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

Statistics tags and overviews are what let a viewer pick a sensible stretch without scanning the whole file.

Code examples #

Example 1: a safe write wrapper #

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

import numpy as np, rasterio
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), rasterio.enums.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()} … {finite.max()}")
    return path

Four properties make this hard to misuse. The dtype and shape come from the array, so a stale profile cannot truncate the data. 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:

from pathlib import Path
import numpy as np, rasterio

def audit_outputs(folder):
    bad = []
    for p in sorted(Path(folder).glob("*.tif")):
        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:<28} {reason:<22} {valid}")
    print(f"{len(bad)} bad of {len(list(Path(folder).glob('*.tif')))}")
    return bad

audit_outputs("out/ndvi/")

  ✗ scene_0042.tif             all NoData             0.0% valid
  ✗ scene_0117.tif             constant 0             100.0% valid
  ✗ scene_0203.tif             suspiciously small (4096 B) 100.0% valid
2 bad of 240

The three reasons point at three different bugs: a scene that was entirely cloud-masked, a computation that produced zeros, and a write that was truncated. Collapsing them into "failed" would lose 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.pyplot as plt
import numpy as np, rasterio

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()
    return fig

show_honestly("ndvi.tif")

The histogram is the part that settles the argument. A raster that renders black but shows a healthy spread of values between −0.3 and 0.9 is a display problem; one whose histogram is a single spike at zero is a data problem. The 2–98 percentile stretch is what QGIS's "Cumulative count cut" does, and applying it manually confirms in one plot whether the default stretch was the whole issue.

Explanation #

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.

Edge cases or notes #

• 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 4 KB GeoTIFF is a header with no image data. Almost always an exception mid-write.
• 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.
• Signed data in an unsigned dtype wraps. −0.31 as uint16 becomes 65535, which renders white, not black — a useful clue.
• 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 among 240 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. QGIS applies a 0–255 stretch by default. Set Min/Max to the real range, or write statistics tags and overviews so the viewer can pick sensibly.

Why is my output only 4 KB? #

The file has a header and no image data, which means the write was interrupted. Look for an exception inside the with block.

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.