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 17a0e0570fd9a23a4a79c4693f1c89b4ea199f10..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)]
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/plotter.py b/kwant/plotter.py
index d164dfcfcf3af4cf211c976420f3ecd9c251b273..913936da21e51edc2735c0e75c329fac03dffaac 100644
--- a/kwant/plotter.py
+++ b/kwant/plotter.py
@@ -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