diff --git a/kwant/_common.py b/kwant/_common.py
index 3af56a12b2618f2cb7f7ba550e0adf13223a6ced..76e952cf6c0d2ffc437d42ca5c093a5825aac0e9 100644
--- a/kwant/_common.py
+++ b/kwant/_common.py
@@ -9,6 +9,8 @@
 import numpy as np
 import numbers
 import inspect
+import warnings
+from contextlib import contextmanager
 
 __all__ = ['KwantDeprecationWarning', 'UserCodeError']
 
@@ -87,6 +89,14 @@ def ensure_rng(rng=None):
                      "numpy.random.RandomState interface.")
 
 
+@contextmanager
+def reraise_warnings(level=3):
+    with warnings.catch_warnings(record=True) as caught_warnings:
+        yield
+    for warning in caught_warnings:
+        warnings.warn(warning.message, stacklevel=level)
+
+
 def get_parameters(func):
     """Get the names of the parameters to 'func' and whether it takes kwargs.
 
diff --git a/kwant/builder.py b/kwant/builder.py
index 096d20c8ff151db1de024f1919fc70ea610c30cd..5b053b01b9ac0436d7a92ffeba2886716d6dfb3d 100644
--- a/kwant/builder.py
+++ b/kwant/builder.py
@@ -20,7 +20,7 @@ from . import system, graph, KwantDeprecationWarning, UserCodeError
 from .linalg import lll
 from .operator import Density
 from .physics import DiscreteSymmetry
-from ._common import ensure_isinstance, get_parameters
+from ._common import ensure_isinstance, get_parameters, reraise_warnings
 
 
 __all__ = ['Builder', 'Site', 'SiteFamily', 'SimpleSiteFamily', 'Symmetry',
@@ -1517,8 +1517,9 @@ class Builder:
         if hop_range > 1:
             # Automatically increase the period, potentially warn the user.
             new_lead = Builder(sym.subgroup((hop_range,)))
-            new_lead.fill(lead_builder, lambda site: True,
-                          lead_builder.sites(), max_sites=float('inf'))
+            with reraise_warnings():
+                new_lead.fill(lead_builder, lambda site: True,
+                              lead_builder.sites(), max_sites=float('inf'))
             lead_builder = new_lead
             sym = lead_builder.symmetry
             H = lead_builder.H
@@ -1568,8 +1569,9 @@ class Builder:
         # system (this one is guaranteed to contain a complete unit cell of the
         # lead). After flood-fill we remove that domain.
         start = {sym.act((max_dom + 1,), site) for site in H}
-        all_added = self.fill(lead_builder, shape, start,
-                              max_sites=float('inf'))
+        with reraise_warnings():
+            all_added = self.fill(lead_builder, shape, start,
+                                  max_sites=float('inf'))
         all_added = [site for site in all_added if site not in start]
         del self[start]
 
diff --git a/kwant/continuum/_common.py b/kwant/continuum/_common.py
index 9b14ef129bb4ae8a12dbf8c80ba5de883e370841..5fd97c6c4cf94b2f0726d019cac7ffa595f2f15e 100644
--- a/kwant/continuum/_common.py
+++ b/kwant/continuum/_common.py
@@ -23,6 +23,8 @@ from sympy.physics.matrices import msigma as _msigma
 
 import warnings
 
+from .._common import reraise_warnings
+
 momentum_operators = sympy.symbols('k_x k_y k_z', commutative=False)
 position_operators = sympy.symbols('x y z', commutative=False)
 
@@ -73,7 +75,8 @@ def lambdify(expr, locals=None):
                [ 0.   ,  0.   ]])
 
     """
-    expr = sympify(expr, locals)
+    with reraise_warnings(level=4):
+        expr = sympify(expr, locals)
 
     args = [s.name for s in expr.atoms(sympy.Symbol)]
     args += [str(f.func) for f in expr.atoms(AppliedUndef, sympy.Function)]
@@ -151,8 +154,10 @@ def sympify(expr, locals=None):
     # if ``expr`` is already a ``sympy`` object we may terminate a code path
     if isinstance(expr, tuple(sympy_classes)):
         if locals:
-            warnings.warn('Input expression is already SymPy object: ' +
-                          '"locals" will not be used.', RuntimeWarning)
+            warnings.warn('Input expression is already SymPy object: '
+                          '"locals" will not be used.',
+                          RuntimeWarning,
+                          stacklevel=2)
         return expr
 
     # if ``expr`` is not a "sympy" then we proceed with sympifying process
diff --git a/kwant/continuum/discretizer.py b/kwant/continuum/discretizer.py
index 5f8705b6e9466b088068a3276f0984a5ff8c7e23..457842978ce57c19505ac85e3cf4022c222da670 100644
--- a/kwant/continuum/discretizer.py
+++ b/kwant/continuum/discretizer.py
@@ -19,6 +19,7 @@ from sympy.printing.precedence import precedence
 from sympy.core.function import AppliedUndef
 
 from .. import builder, lattice
+from .._common import reraise_warnings
 from ._common import (sympify, gcd, position_operators, momentum_operators,
                       monomials)
 
@@ -171,7 +172,8 @@ def discretize_symbolic(hamiltonian, coords=None, *, locals=None):
         The coordinates that have been discretized.
 
     """
-    hamiltonian = sympify(hamiltonian, locals)
+    with reraise_warnings():
+        hamiltonian = sympify(hamiltonian, locals)
 
     atoms_names = [s.name for s in hamiltonian.atoms(sympy.Symbol)]
     if any( s == 'a' for s in atoms_names):
@@ -276,8 +278,9 @@ def build_discretized(tb_hamiltonian, coords, *, grid_spacing=1, locals=None):
     if len(coords) == 0:
         raise ValueError('Discrete coordinates cannot be empty.')
 
-    for k, v in tb_hamiltonian.items():
-        tb_hamiltonian[k] = sympify(v, locals)
+    with reraise_warnings():
+        for k, v in tb_hamiltonian.items():
+            tb_hamiltonian[k] = sympify(v, locals)
 
     coords = list(coords)
     if coords != sorted(coords):
diff --git a/kwant/linalg/mumps.py b/kwant/linalg/mumps.py
index 397b4cc963bc9f9d593bce12fcd2c41904bc5912..465ef1aa359ab3c88adbf5a4a4ccf125c73094dd 100644
--- a/kwant/linalg/mumps.py
+++ b/kwant/linalg/mumps.py
@@ -299,7 +299,9 @@ class MUMPSContext:
         if reuse_analysis:
             if self.mumps_instance is None:
                 warnings.warn("Missing analysis although reuse_analysis=True. "
-                              "New analysis is performed.", RuntimeWarning)
+                              "New analysis is performed.",
+                              RuntimeWarning,
+                              stacklevel=2)
                 self.analyze(a, ordering=ordering, overwrite_a=overwrite_a)
             else:
                 dtype, row, col, data = _make_assembled_from_coo(a,
diff --git a/kwant/plotter.py b/kwant/plotter.py
index 4c17c8bbc1f06b108e3f3ee17700194e20b0b508..913936da21e51edc2735c0e75c329fac03dffaac 100644
--- a/kwant/plotter.py
+++ b/kwant/plotter.py
@@ -67,7 +67,7 @@ def matplotlib_chores():
         warnings.warn("Matplotlib 1.4.0 has a bug that makes 3D plotting "
                       "unusable (2D plotting is not affected). Please "
                       "consider using a different version of matplotlib.",
-                      RuntimeWarning)
+                      RuntimeWarning, stacklevel=2)
 
     pre_1_4_matplotlib = [int(x) for x in ver.split('.')[:2]] < [1, 4]
 
@@ -1457,7 +1457,7 @@ def mask_interpolate(coords, values, a=None, method='nearest', oversampling=3):
     if min_dist < 1e-6 * np.linalg.norm(cmax - cmin):
         warnings.warn("Some sites have nearly coinciding positions, "
                       "interpolation may be confusing.",
-                      RuntimeWarning)
+                      RuntimeWarning, stacklevel=2)
 
     if a is None:
         a = min_dist
@@ -1576,7 +1576,8 @@ def map(sys, value, colorbar=True, cmap=None, vmin=None, vmax=None, a=None,
             raise ValueError('List of values is only allowed as input '
                              'for finalized systems.')
     value = np.array(value)
-    img, min, max = mask_interpolate(coords, value, a, method, oversampling)
+    with _common.reraise_warnings():
+        img, min, max = mask_interpolate(coords, value, a, method, oversampling)
     border = 0.5 * (max - min) / (np.asarray(img.shape) - 1)
     min -= border
     max += border