diff --git a/docs/sphinx/source/whatsnew/v0.16.2.rst b/docs/sphinx/source/whatsnew/v0.16.2.rst
index 2c243d0921..64132a2547 100644
--- a/docs/sphinx/source/whatsnew/v0.16.2.rst
+++ b/docs/sphinx/source/whatsnew/v0.16.2.rst
@@ -19,6 +19,8 @@ Bug fixes
   (:issue:`1959`, :pull:`2912`)
 * :py:func:`pvlib.tracking.calc_surface_orientation` returns scalar when input
   is scalar (:issue:`2881`, :pull:`2882`)
+* :py:func:`pvlib.scaling.wvm` no longer returns all NaNs when ``positions``
+  contains a single site (:issue:`2895`, :pull:`2896`)
 
 Enhancements
 ~~~~~~~~~~~~
@@ -65,7 +67,5 @@ Contributors
 * Chen Yi (:ghuser:`chenyi6758-source`)
 * Rajiv Daxini (:ghuser:`rdaxini`)
 * jason-writer (:ghuser:`jason-writer`)
+* :ghuser:`bfs2021`
 * :ghuser:`xiajiahang`
-
-
-
diff --git a/pvlib/scaling.py b/pvlib/scaling.py
index 693d9f6bff..67d10f17e5 100644
--- a/pvlib/scaling.py
+++ b/pvlib/scaling.py
@@ -23,10 +23,10 @@ def wvm(clearsky_index, positions, cloud_speed, dt=None):
     clearsky_index : numeric or pandas.Series
         Clear Sky Index time series that will be smoothed.
 
-    positions : numeric
-        Array of coordinate distances as (x,y) pairs representing the
-        easting, northing of the site positions in meters [m]. Distributed
-        plants could be simulated by gridded points throughout the plant
+    positions : np.array
+        Array of shape (N, 2) of coordinates (x, y) representing the
+        easting, northing in meters [m] of the PV systems. Large PV
+        plants could be simulated by gridded points over the plant's
         footprint.
 
     cloud_speed : numeric
@@ -93,10 +93,10 @@ def _compute_vr(positions, cloud_speed, tmscales):
 
     Parameters
     ----------
-    positions : numeric
-        Array of coordinate distances as (x,y) pairs representing the
-        easting, northing of the site positions in meters [m]. Distributed
-        plants could be simulated by gridded points throughout the plant
+    positions : np.array
+        Array of shape (N, 2) of coordinates (x, y) representing the
+        easting, northing in meters [m] of the PV systems. Large PV
+        plants could be simulated by gridded points over the plant's
         footprint.
 
     cloud_speed : numeric
@@ -127,6 +127,11 @@ def _compute_vr(positions, cloud_speed, tmscales):
     # Added by Joe Ranalli (@jranalli), Penn State Hazleton, 2021
 
     pos = np.array(positions)
+    if len(pos) == 1:
+        # A single site has no pairwise distances, so the variability
+        # reduction is 1; the effective-size optimizer below would converge
+        # to n_dist=0 and divide by zero. See GH #2895.
+        return np.ones(tmscales.shape)
     dist = pdist(pos, 'euclidean')
 
     # Find effective length of position vector, 'dist' is full pairwise
diff --git a/tests/test_scaling.py b/tests/test_scaling.py
index f8638f580a..f778f0e65b 100644
--- a/tests/test_scaling.py
+++ b/tests/test_scaling.py
@@ -207,3 +207,17 @@ def test_wvm_series_xyaslist(clear_sky_index, time, positions,
 def test_wvm_invalid(clear_sky_index, positions):
     with pytest.raises(ValueError):
         scaling.wvm(clear_sky_index, positions, cloud_speed)
+
+
+def test_compute_vr_single_site(expect_tmscale):
+    vr = scaling._compute_vr(np.array([(0., 0.)]), cloud_speed,
+                             np.array(expect_tmscale))
+    assert_almost_equal(vr, np.ones(len(expect_tmscale)))
+
+
+def test_wvm_single_site(clear_sky_index):
+    # A single site has no variability reduction: the smoothed signal is the
+    # input unchanged. Previously returned all NaNs. See GH #2895.
+    single_site = np.array([(0., 0.)])
+    cs_sm, _, _ = scaling.wvm(clear_sky_index, single_site, cloud_speed, dt=dt)
+    assert_almost_equal(cs_sm[5000:5005], clear_sky_index[5000:5005])