diff --git a/doc/source/images/closed_system.py.diff b/doc/source/images/closed_system.py.diff
index bfa0d09a472436b1fa9766d9813ea1cac97ca528..fb27cedf90497b80fafd098721b63807b03a566f 100644
--- a/doc/source/images/closed_system.py.diff
+++ b/doc/source/images/closed_system.py.diff
@@ -8,7 +8,7 @@
  
  
  def make_system(a=1, t=1.0, r=10):
-@@ -65,28 +66,40 @@
+@@ -65,29 +66,40 @@
  
          energies.append(ev)
  
@@ -39,13 +39,14 @@
      evecs = sla.eigsh(ham_mat, k=20, which='SM')[1]
  
      # Plot the probability density of the 10th eigenmode.
--    kwant.plotter.map(sys, np.abs(evecs[:, 9])**2, colorbar=False)
-+    kwant.plotter.map(sys, np.abs(evecs[:, 9])**2, colorbar=False,
-+                      file="closed_system_eigenvector.pdf",
-+                      fig_size=size, dpi=_defs.dpi)
-+    kwant.plotter.map(sys, np.abs(evecs[:, 9])**2, colorbar=False,
-+                      file="closed_system_eigenvector.png",
-+                      fig_size=size, dpi=_defs.dpi)
+-    kwant.plotter.map(sys, np.abs(evecs[:, 9])**2,
+-                      colorbar=False, oversampling=1)
++    kwant.plotter.map(
++        sys, np.abs(evecs[:, 9])**2, colorbar=False, oversampling=1,
++        file="closed_system_eigenvector.pdf", fig_size=size, dpi=_defs.dpi)
++    kwant.plotter.map(
++        sys, np.abs(evecs[:, 9])**2, colorbar=False, oversampling=1,
++        file="closed_system_eigenvector.png", fig_size=size, dpi=_defs.dpi)
  
  
  def main():
diff --git a/doc/source/tutorial/closed_system.py b/doc/source/tutorial/closed_system.py
index 8c70ff556f75dd865fe9e33be49583229e565d3e..8b46ece234e8b4535d4d898ca42439ba8dca9af8 100644
--- a/doc/source/tutorial/closed_system.py
+++ b/doc/source/tutorial/closed_system.py
@@ -85,7 +85,8 @@ def plot_wave_function(sys):
     evecs = sla.eigsh(ham_mat, k=20, which='SM')[1]
 
     # Plot the probability density of the 10th eigenmode.
-    kwant.plotter.map(sys, np.abs(evecs[:, 9])**2, colorbar=False)
+    kwant.plotter.map(sys, np.abs(evecs[:, 9])**2,
+                      colorbar=False, oversampling=1)
 #HIDDEN_END_wave
 
 
diff --git a/doc/source/tutorial/tutorial3.rst b/doc/source/tutorial/tutorial3.rst
index 814ab761abb81f2aaf47913fc0df3205a64b46b6..d0a86b6595c942660cd9eeeb845c9b2a37b7bd71 100644
--- a/doc/source/tutorial/tutorial3.rst
+++ b/doc/source/tutorial/tutorial3.rst
@@ -109,10 +109,12 @@ using `~kwant.plotter.map`:
     :start-after: #HIDDEN_BEGIN_wave
     :end-before: #HIDDEN_END_wave
 
-This yields the result:
-
 .. image:: ../images/closed_system_eigenvector.*
 
+The last two arguments to `~kwant.plotter.map` are optional.  The first prevents
+a colorbar from appearing.  The second, ``oversampling=1``, makes the image look
+better for the special case of a square lattice.
+
 .. seealso::
     The full source code can be found in
     :download:`tutorial/closed_system.py <../../../tutorial/closed_system.py>`
diff --git a/kwant/plotter.py b/kwant/plotter.py
index ae4371cc1e618cccf9eb4c556aaeed689f48decf..22d4d3097ecf8fb20331c89269b6c797fb155ca0 100644
--- a/kwant/plotter.py
+++ b/kwant/plotter.py
@@ -743,8 +743,8 @@ def mask_interpolate(coords, values, a=None, method='nearest', oversampling=3):
       exactly at the center of a pixel of the output array.
 
     - When plotting a system on a square lattice and `method` is "nearest", it
-      makes sense to set `oversampling` to ``1`` to minimize the size of the
-      output array.
+      makes sense to set `oversampling` to ``1``.  Then, each site will
+      correspond to exactly one pixel in the resulting array.
     """
     # Build the bounding box.
     cmin, cmax = coords.min(0), coords.max(0)
@@ -766,7 +766,12 @@ def mask_interpolate(coords, values, a=None, method='nearest', oversampling=3):
     grid = tuple(np.ogrid[dims])
     img = interpolate.griddata(coords, values, grid, method)
     mask = np.mgrid[dims].reshape(len(cmin), -1).T
-    mask = tree.query(mask, eps=1.)[0] > 1.5 * a
+    # The numerical values in the following line are optimized for the common
+    # case of a square lattice:
+    # * 0.99 makes sure that non-masked pixels and sites correspond 1-by-1 to
+    #   each other when oversampling == 1.
+    # * 0.4 (which is just below sqrt(2) - 1) makes tree.query() exact.
+    mask = tree.query(mask, eps=0.4)[0] > 0.99 * a
 
     return np.ma.masked_array(img, mask), cmin, cmax
 
@@ -814,10 +819,9 @@ def map(sys, value, colorbar=True, cmap=None,
 
     Notes
     -----
-    - See notes of `~kwant.plotter.show_interpolate`.
-
-    - Matplotlib's interpolation is turned off, if the keyword argument
-      `method` is not set or set to the default value "nearest".
+    - When plotting a system on a square lattice and `method` is "nearest", it
+      makes sense to set `oversampling` to ``1``.  Then, each site will
+      correspond to exactly one pixel.
     """
     sites = sys_leads_sites(sys, 0)
     coords = sys_leads_pos(sys, sites)