diff --git a/src/spatialdata_plot/pl/render.py b/src/spatialdata_plot/pl/render.py
index 3f276beb..8e8e32f4 100644
--- a/src/spatialdata_plot/pl/render.py
+++ b/src/spatialdata_plot/pl/render.py
@@ -1314,6 +1314,21 @@ def _render_images(
 
     n_channels = len(channels)
 
+    # Reject NaN early: silent substitution (na_color in 1ch, black in multi-channel)
+    # hides upstream data problems.
+    nan_channels: list[Any] = []
+    for ch in channels:
+        layer = img.sel(c=ch) if isinstance(ch, str) else img.isel(c=ch)
+        if np.issubdtype(layer.dtype, np.floating) and np.isnan(layer.values).any():
+            nan_channels.append(ch)
+    if nan_channels:
+        raise ValueError(
+            f"Image '{render_params.element}' contains NaN pixels in channel(s) {nan_channels}. "
+            "NaN is not supported by render_images. Replace NaN before plotting, e.g. "
+            f"`sdata.images['{render_params.element}'] = sdata.images['{render_params.element}'].fillna(0)`, "
+            "or mask the affected region."
+        )
+
     # When grayscale was applied and user didn't provide an explicit cmap,
     # default to "gray" for intuitive single-channel rendering.
     got_multiple_cmaps = isinstance(render_params.cmap_params, list)
diff --git a/tests/pl/test_render_images.py b/tests/pl/test_render_images.py
index fb723c2f..7ba01b55 100644
--- a/tests/pl/test_render_images.py
+++ b/tests/pl/test_render_images.py
@@ -531,6 +531,42 @@ def test_cmap_matches_selected_channels_not_full_image(sdata_blobs: SpatialData)
     plt.close(fig)
 
 
+# Regression for #628: NaN pixels must raise, not silently render
+# (na_color in 1ch, black in multi-channel).
+def _nan_image(n_channels: int, nan_indices: list[int]) -> SpatialData:
+    rng = np.random.default_rng(0)
+    data = rng.uniform(0, 1, (n_channels, 8, 8)).astype(np.float32)
+    for ch in nan_indices:
+        data[ch, 0:3, 0:3] = np.nan
+    img = Image2DModel.parse(data, dims=("c", "y", "x"), c_coords=list(range(n_channels)))
+    return SpatialData(images={"img": img})
+
+
+def test_nan_in_multi_channel_raises():
+    # Message must list the offending channel and include the fillna hint.
+    sdata = _nan_image(n_channels=2, nan_indices=[0])
+    with pytest.raises(ValueError, match=r"NaN.*channel\(s\) \[0\].*fillna"):
+        sdata.pl.render_images("img").pl.show()
+
+
+def test_nan_in_single_channel_raises():
+    # 1ch previously substituted na_color silently; locks the new symmetric behavior.
+    sdata = _nan_image(n_channels=1, nan_indices=[0])
+    with pytest.raises(ValueError, match="NaN"):
+        sdata.pl.render_images("img").pl.show()
+
+
+def test_integer_dtype_skips_nan_check():
+    # Integer-dtype images can't contain NaN; the check must short-circuit on dtype.
+    rng = np.random.default_rng(0)
+    data = rng.integers(0, 255, (2, 8, 8), dtype=np.uint16)
+    img = Image2DModel.parse(data, dims=("c", "y", "x"), c_coords=[0, 1])
+    sdata = SpatialData(images={"img": img})
+    fig, ax = plt.subplots()
+    sdata.pl.render_images("img").pl.show(ax=ax)
+    plt.close(fig)
+
+
 # Regression for #612: vmin/vmax kwargs are no longer accepted on any render
 # function. The check covers all four to prevent the asymmetry from re-emerging.
 @pytest.mark.parametrize("kwarg", ["vmin", "vmax"])