diff --git a/yt/visualization/base_plot_types.py b/yt/visualization/base_plot_types.py
index aec211d5ebb..211fb869894 100644
--- a/yt/visualization/base_plot_types.py
+++ b/yt/visualization/base_plot_types.py
@@ -219,22 +219,22 @@ def _init_image(self, data, cbnorm, cblinthresh, cmap, extent, aspect):
 
         if cbnorm == "symlog":
             # if cblinthresh is not specified, try to come up with a reasonable default
-            min_abs_val, max_abs_val = np.sort(np.abs((zmin, zmax)))
+            min_abs_val, max_abs_val = np.sort(
+                np.abs((np.nanmin(data), np.nanmax(data)))
+            )
             if cblinthresh is not None:
-                pass
+                if zmin * zmax > 0 and cblinthresh < min_abs_val:
+                    # see https://github.com/yt-project/yt/issues/3564
+                    warnings.warn(
+                        f"Cannot set a symlog norm with linear threshold {cblinthresh} "
+                        f"lower than the minimal absolute data value {min_abs_val} . "
+                        "Switching to log norm."
+                    )
+                    cbnorm = "log10"
             elif min_abs_val > 0:
                 cblinthresh = min_abs_val
-            elif zmin * zmax > 0 and cblinthresh < min_abs_val:
-                warnings.warn(
-                    f"Cannot set a symlog norm with linear threshold {cblinthresh} "
-                    f"lower than the minimal absolute data value {min_abs_val} . "
-                    "Switching to log norm."
-                )
-                cbnorm = "log10"
             else:
                 cblinthresh = max_abs_val / 1000
-            if cblinthresh is None:
-                cblinthresh = np.nanmin(np.absolute(data)[data != 0])
 
         if cbnorm == "log10":
             cbnorm_cls = matplotlib.colors.LogNorm
diff --git a/yt/visualization/tests/test_plotwindow.py b/yt/visualization/tests/test_plotwindow.py
index 550556cf564..2e9d8437837 100644
--- a/yt/visualization/tests/test_plotwindow.py
+++ b/yt/visualization/tests/test_plotwindow.py
@@ -722,6 +722,43 @@ def _neg_density(field, data):
             plot.save(f.name)
 
 
+def test_symlog_min_zero():
+    # see https://github.com/yt-project/yt/issues/3791
+    shape = (32, 16, 1)
+    a = np.linspace(0, 1, 16)
+    b = np.ones((32, 16))
+    c = np.reshape(a * b, shape)
+    data = {("gas", "density"): c}
+
+    ds = load_uniform_grid(
+        data,
+        shape,
+        bbox=np.array([[0.0, 5.0], [0, 1], [-0.1, +0.1]]),
+    )
+
+    p = SlicePlot(ds, "z", "density")
+    im_arr = p["gas", "density"].image.get_array()
+
+    # check that no data value was mapped to a NaN (log(0))
+    assert np.all(~np.isnan(im_arr))
+    # 0 should be mapped to itself since we expect a symlog norm
+    assert np.min(im_arr) == 0.0
+
+
+def test_symlog_extremely_small_vals():
+    # check that the plot can be constructed without crashing
+    # see https://github.com/yt-project/yt/issues/3858
+    shape = (64, 64, 1)
+    arr = np.full(shape, 5.0e-324)
+    arr[0, 0] = -1e12
+    arr[1, 1] = 200
+    d = {"scalar": arr}
+
+    ds = load_uniform_grid(d, shape)
+    p = SlicePlot(ds, "z", ("stream", "scalar"))
+    p["stream", "scalar"]
+
+
 def test_nan_data():
     data = np.random.random((16, 16, 16)) - 0.5
     data[:9, :9, :9] = np.nan