diff --git a/python/lsst/ip/diffim/getTemplate.py b/python/lsst/ip/diffim/getTemplate.py
index 04e9d5ce6..1d1d12523 100644
--- a/python/lsst/ip/diffim/getTemplate.py
+++ b/python/lsst/ip/diffim/getTemplate.py
@@ -375,7 +375,8 @@ def _makeExposureCatalog(self, exposures, dataIds):
 
         return images, catalog, totalBox
 
-    def _merge(self, maskedImages, bbox, wcs):
+    @staticmethod
+    def _merge(maskedImages, bbox, wcs):
         """Merge the images that came from one tract into one larger image,
         ignoring NaN pixels and non-finite variance pixels from individual
         exposures.
@@ -398,8 +399,11 @@ def _merge(self, maskedImages, bbox, wcs):
         merged = afwImage.ExposureF(bbox, wcs)
         weights = afwImage.ImageF(bbox)
         for maskedImage in maskedImages:
-            weight = afwImage.ImageF(maskedImage.variance.array**(-0.5))
-            bad = np.isnan(maskedImage.image.array) | ~np.isfinite(maskedImage.variance.array)
+            # Catch both zero-value and NaN variance plane pixels
+            good = maskedImage.variance.array > 0
+            weight = afwImage.ImageF(maskedImage.getBBox())
+            weight.array[good] = maskedImage.variance.array[good]**(-0.5)
+            bad = np.isnan(maskedImage.image.array) | ~good
             # Note that modifying the patch MaskedImage in place is fine;
             # we're throwing it away at the end anyway.
             maskedImage.image.array[bad] = 0.0
@@ -412,16 +416,18 @@ def _merge(self, maskedImages, bbox, wcs):
             maskedImage.image *= weight
             maskedImage.variance *= weight
             merged.maskedImage[maskedImage.getBBox()] += maskedImage
-            # Clear the NaNs to ensure that areas missing from this input are
-            # masked with NO_DATA after the loop.
-            weight.array[np.isnan(weight.array)] = 0
             weights[maskedImage.getBBox()] += weight
-        # Cannot use `merged.maskedImage /= weights` because that operator
-        # divides the variance by the weight twice; in this case `weights` are
-        # the exact values we want to scale by.
-        merged.image /= weights
-        merged.variance /= weights
-        merged.mask.array |= merged.mask.getPlaneBitMask("NO_DATA") * (weights.array == 0)
+
+        inverseWeights = np.zeros_like(weights.array)
+        good = weights.array > 0
+        inverseWeights[good] = 1/weights.array[good]
+
+        # Cannot use `merged.maskedImage *= inverseWeights` because that
+        # operator divides the variance by the weight twice; in this case
+        # `inverseWeights` are the exact values we want to scale by.
+        merged.image.array *= inverseWeights
+        merged.variance.array *= inverseWeights
+        merged.mask.array |= merged.mask.getPlaneBitMask("NO_DATA") * (inverseWeights == 0)
 
         return merged
 
diff --git a/tests/test_getTemplate.py b/tests/test_getTemplate.py
index a9b00e9b1..1668c8b13 100644
--- a/tests/test_getTemplate.py
+++ b/tests/test_getTemplate.py
@@ -282,9 +282,37 @@ def testMissingPatches(self):
                           dataIds=self.dataIds,
                           physical_filter="a_test")
         no_data = (result.template.mask.array & result.template.mask.getPlaneBitMask("NO_DATA")) != 0
-        self.assertTrue(all(np.isnan(result.template.image.array[no_data])))
-        self.assertTrue(all(np.isnan(result.template.variance.array[no_data])))
-        self.assertEqual(no_data.sum(), 21548)
+        self.assertTrue(np.isfinite(result.template.image.array).all())
+        self.assertTrue(np.isfinite(result.template.variance.array).all())
+        self.assertEqual(no_data.sum(), 20990)
+
+    @lsst.utils.tests.methodParameters(
+        box=[
+            lsst.geom.Box2I(lsst.geom.Point2I(0, 0), lsst.geom.Point2I(180, 180)),
+            lsst.geom.Box2I(lsst.geom.Point2I(200, 200), lsst.geom.Point2I(600, 600)),
+        ],
+        nInput=[8, 16],
+    )
+    def testNanInputs(self, box=None, nInput=None):
+        """Test that the template has finite values when some of the input pixels have NaN as variance."""
+        for tract, patchCoadds in self.patches.items():
+            for patchCoadd in patchCoadds:
+                bbox = lsst.geom.Box2I()
+                bbox.include(lsst.geom.Point2I(patchCoadd.getBBox().getCenter()))
+                bbox.grow(3)
+                patchCoadd.variance[bbox].array *= np.nan
+
+        box = lsst.geom.Box2I(lsst.geom.Point2I(200, 200), lsst.geom.Point2I(600, 600))
+        task = lsst.ip.diffim.GetTemplateTask()
+        result = task.run(coaddExposures=self.patches,
+                          bbox=lsst.geom.Box2I(box),
+                          wcs=self.exposure.wcs,
+                          dataIds=self.dataIds,
+                          physical_filter="a_test")
+        self._checkMetadata(result.template, task.config, box, self.exposure.wcs, 16)
+        # We just check that the pixel values are all finite. We cannot check that pixel values
+        # in the template are closer to the original anymore.
+        self.assertTrue(np.isfinite(result.template.image.array).all())
 
 
 def setup_module(module):