diff --git a/pvlib/bifacial/utils.py b/pvlib/bifacial/utils.py
index 8d7f5d5a71..c8d1f8cb48 100644
--- a/pvlib/bifacial/utils.py
+++ b/pvlib/bifacial/utils.py
@@ -4,7 +4,8 @@
 """
 import numpy as np
 from pvlib.tools import sind, cosd, tand
-
+import warnings
+from pvlib._deprecation import pvlibDeprecationWarning
 
 def _solar_projection_tangent(solar_zenith, solar_azimuth, surface_azimuth):
     """
@@ -172,8 +173,34 @@ def vf_ground_sky_2d(rotation, gcr, x, pitch, height, max_rows=10):
     return vf
 
 
+def _dist(p1, p2):
+    return ((p1[0] - p2[0])**2 + (p1[1] - p2[1])**2)**0.5
+
+
+def _angle(p1, p2):
+    return np.arctan2(p2[1] - p1[1], p2[0] - p1[0])
+
+
+def _obstructed_string_length(p1, p2, ob_left, ob_right):
+    # string length calculations for Hottel's crossed strings method,
+    # considering view obstructions from the left and right.
+    # all inputs are (x, y) points.
+
+    # unobstructed length
+    d = _dist(p1, p2)
+    # obstructed on the left
+    d = np.where(_angle(p1, p2) > _angle(p1, ob_left),
+                 _dist(p1, ob_left) + _dist(ob_left, p2),
+                 d)
+    # obstructed on the right
+    d = np.where(_angle(p1, p2) < _angle(p1, ob_right),
+                 _dist(p1, ob_right) + _dist(ob_right, p2),
+                 d)
+    return d
+
+
 def vf_ground_sky_2d_integ(surface_tilt, gcr, height, pitch, max_rows=10,
-                           npoints=100, vectorize=False):
+                           npoints=None, vectorize=None):
     """
     Integrated view factor to the sky from the ground underneath
     interior rows of the array.
@@ -204,23 +231,49 @@ def vf_ground_sky_2d_integ(surface_tilt, gcr, height, pitch, max_rows=10,
         Integration of view factor over the length between adjacent, interior
         rows.  Shape matches that of ``surface_tilt``. [unitless]
     """
-    # Abuse vf_ground_sky_2d by supplying surface_tilt in place
-    # of a signed rotation. This is OK because
-    # 1) z span the full distance between 2 rows, and
-    # 2) max_rows is set to be large upstream, and
-    # 3) _vf_ground_sky_2d considers [-max_rows, +max_rows]
-    # The VFs to the sky will thus be symmetric around z=0.5
-    z = np.linspace(0, 1, npoints)
-    rotation = np.atleast_1d(surface_tilt)
-    if vectorize:
-        fz_sky = vf_ground_sky_2d(rotation, gcr, z, pitch, height, max_rows)
-    else:
-        fz_sky = np.zeros((npoints, len(rotation)))
-        for k, r in enumerate(rotation):
-            vf = vf_ground_sky_2d(r, gcr, z, pitch, height, max_rows)
-            fz_sky[:, k] = vf[:, 0]  # remove spurious rotation dimension
-    # calculate the integrated view factor for all of the ground between rows
-    return np.trapz(fz_sky, z, axis=0)
+    if npoints is not None or vectorize is not None:
+        msg = (
+            "The `npoints` and `vectorize` parameters have no effect and will "
+            "be removed in a future version."  # TODO make this better
+        )
+        warnings.warn(msg, pvlibDeprecationWarning)
+
+    input_is_scalar = np.isscalar(surface_tilt)
+
+    collector_width = pitch * gcr
+    surface_tilt = np.atleast_2d(np.abs(surface_tilt))
+
+    # TODO figure out if this range is correct, or if the original code has a bug
+    k = np.arange(-max_rows+1, max_rows+1)[:, np.newaxis]
+
+    # primary crossed string points:
+    # a, b: boundaries of ground segment
+    # c, d: upper module edges
+    a = (0, 0)
+    b = (pitch, 0)
+    c = ((k+1)*pitch - 0.5 * collector_width * cosd(surface_tilt),
+         height + 0.5 * collector_width * sind(surface_tilt))
+    d = (c[0] - pitch, c[1])
+
+    # view obstruction points (lower module edges)
+    obs_left = (d[0] + collector_width * cosd(surface_tilt),
+                d[1] - collector_width * sind(surface_tilt))
+    obs_right = (obs_left[0] + pitch, obs_left[1])
+
+    # hottel string lengths, considering obstructions
+    ac = _obstructed_string_length(a, c, obs_left, obs_right)
+    ad = _obstructed_string_length(a, d, obs_left, obs_right)
+    bc = _obstructed_string_length(b, c, obs_left, obs_right)
+    bd = _obstructed_string_length(b, d, obs_left, obs_right)
+
+    # crossed string formula for VF
+    vf_per_slat = np.maximum(0.5 * (1/pitch) * ((ac + bd) - (bc + ad)), 0)
+    vf_total = np.sum(vf_per_slat, axis=0)
+
+    if input_is_scalar:
+        vf_total = vf_total.item()
+
+    return vf_total
 
 
 def _vf_poly(surface_tilt, gcr, x, delta):