diff --git a/pvlib/irradiance.py b/pvlib/irradiance.py
index eb63371b81..e012059940 100644
--- a/pvlib/irradiance.py
+++ b/pvlib/irradiance.py
@@ -15,6 +15,7 @@
 
 from pvlib import atmosphere, solarposition, tools
 from pvlib._deprecation import deprecated
+from tools import nanmaximum
 
 # see References section of grounddiffuse function
 SURFACE_ALBEDOS = {'urban': 0.18,
@@ -1147,24 +1148,24 @@ def perez(surface_tilt, surface_azimuth, dhi, dni, dni_extra,
     F2c = np.vstack((F2c, nans))
 
     F1 = (F1c[ebin, 0] + F1c[ebin, 1] * delta + F1c[ebin, 2] * z)
-    F1 = np.maximum(F1, 0)
+    F1 = nanmaximum(F1, 0)
 
     F2 = (F2c[ebin, 0] + F2c[ebin, 1] * delta + F2c[ebin, 2] * z)
-    F2 = np.maximum(F2, 0)
+    F2 = nanmaximum(F2, 0)
 
     A = aoi_projection(surface_tilt, surface_azimuth,
                        solar_zenith, solar_azimuth)
-    A = np.maximum(A, 0)
+    A = nanmaximum(A, 0)
 
     B = tools.cosd(solar_zenith)
-    B = np.maximum(B, tools.cosd(85))
+    B = nanmaximum(B, tools.cosd(85))
 
     # Calculate Diffuse POA from sky dome
     term1 = 0.5 * (1 - F1) * (1 + tools.cosd(surface_tilt))
     term2 = F1 * A / B
     term3 = F2 * tools.sind(surface_tilt)
 
-    sky_diffuse = np.maximum(dhi * (term1 + term2 + term3), 0)
+    sky_diffuse = nanmaximum(dhi * (term1 + term2 + term3), 0)
 
     # we've preserved the input type until now, so don't ruin it!
     if isinstance(sky_diffuse, pd.Series):
diff --git a/pvlib/tools.py b/pvlib/tools.py
index fa3d368369..9d9c69cc5e 100644
--- a/pvlib/tools.py
+++ b/pvlib/tools.py
@@ -86,6 +86,25 @@ def asind(number):
     return res
 
 
+def nanmaximum(x1, x2):
+    """
+    Element-wise maximum preserving nans, without numpy warning
+
+    Parameters
+    ----------
+    array1, array2 : array_like
+        the arrays to be compared, must have the same shape or broadcast to the
+        same shape.
+
+    Returns
+    -------
+    y : array_like
+        the maximum of x1 and x2, element-wise. nan is preserved.
+    """
+
+    return np.where((np.isnan(x1) | np.isnan(x2)), np.nan, np.fmax(x1, x2))
+
+
 def localize_to_utc(time, location):
     """
     Converts or localizes a time series to UTC.