diff --git a/doc/source/pre/whatsnew/1.4.rst b/doc/source/pre/whatsnew/1.4.rst
index c962fb9741f8f6ac8709aa8f8bd62b8ced65f377..a2ea4f32b4a5f99232b61948a5c3dde8c1ba8866 100644
--- a/doc/source/pre/whatsnew/1.4.rst
+++ b/doc/source/pre/whatsnew/1.4.rst
@@ -35,12 +35,12 @@ ensure that the chosen gauge is consistent across the whole system
 directions). This release introduces `kwant.physics.magnetic_gauge`,
 which calculates the Peierls phases for you::
 
-  import numpy as np
-  import kwant
-
   def hopping(a, b, t, peierls):
       return -t * peierls(a, b)
 
+  def B_syst(pos):
+     return np.exp(-np.sum(pos * pos))
+
   syst = make_system(hopping)
   lead = make_lead(hopping).substituted(peierls='peierls_lead')
   syst.attach_lead(lead)
@@ -48,11 +48,7 @@ which calculates the Peierls phases for you::
 
   gauge = kwant.physics.magnetic_gauge(syst)
 
-  def B_syst(pos):
-     return np.exp(-np.sum(pos * pos))
-
   # B_syst in scattering region, 0 in lead.
-  # Ensure that the fields match at the system/lead interface!
   peierls_syst, peierls_lead = gauge(B_syst, 0)
 
   params = dict(t=1, peierls=peierls_syst, peierls_lead=peierls_lead)
diff --git a/kwant/builder.py b/kwant/builder.py
index ad482a75c838438b709296704ef3ff8ee2419022..79e1f0198c0d2f882bd9cbd736117c2e572e8d08 100644
--- a/kwant/builder.py
+++ b/kwant/builder.py
@@ -1980,7 +1980,9 @@ class FiniteSystem(_FinalizedBuilderMixin, system.FiniteSystem):
         to the integer-labeled site ``i`` of the low-level system. The sites
         are ordered first by their family and then by their tag.
     id_by_site : dict
-        The inverse of ``sites``; maps from ``i`` to ``sites[i]``.
+        The inverse of ``sites``; maps high-level `~kwant.builder.Site`
+        instances to their integer label.
+        Satisfies ``id_by_site[sites[i]] == i``.
     """
 
     def __init__(self, builder):
@@ -2097,7 +2099,9 @@ class InfiniteSystem(_FinalizedBuilderMixin, system.InfiniteSystem):
         ``sites[i]`` is the `~kwant.builder.Site` instance that corresponds
         to the integer-labeled site ``i`` of the low-level system.
     id_by_site : dict
-        The inverse of ``sites``; maps from ``i`` to ``sites[i]``.
+        The inverse of ``sites``; maps high-level `~kwant.builder.Site`
+        instances to their integer label.
+        Satisfies ``id_by_site[sites[i]] == i``.
 
     Notes
     -----
diff --git a/kwant/kpm.py b/kwant/kpm.py
index eb161de374b5a1276a7e2030b14a009cd7cd656c..a06d7c6d19af5a3a5266e5fff68103cc9f7b3b04 100644
--- a/kwant/kpm.py
+++ b/kwant/kpm.py
@@ -945,7 +945,7 @@ class LocalVectors:
         must be a list of integers with the indices where column vectors
         are nonzero.
     """
-    def __init__(self, syst, where, *args):
+    def __init__(self, syst, where=None, *args):
         self.tot_norbs, self.orbs = _normalize_orbs_where(syst, where)
         self._idx = 0
 
@@ -1006,7 +1006,11 @@ def _normalize_orbs_where(syst, where):
         tot_norbs = _get_tot_norbs(syst)
         orbs = _from_where_to_orbs(syst, where)
     else:
-        tot_norbs = csr_matrix(syst).shape[0]
+        try:
+            tot_norbs = csr_matrix(syst).shape[0]
+        except TypeError:
+            raise TypeError("'syst' is neither a matrix "
+                             "nor a Kwant system.")
         orbs = (range(tot_norbs) if where is None
                 else np.asarray(where, int))
     return tot_norbs, orbs
diff --git a/kwant/physics/tests/test_leads.py b/kwant/physics/tests/test_leads.py
index bc41f671811bd6c2cdae99d48c7b82fa8ebeeb3d..e7994941f9e7375d2b292fe8163ace4df3138553 100644
--- a/kwant/physics/tests/test_leads.py
+++ b/kwant/physics/tests/test_leads.py
@@ -281,7 +281,7 @@ def test_modes_bearded_ribbon():
 def check_equivalence(h, t, n, sym='', particle_hole=None, chiral=None,
                       time_reversal=None):
     """Compare modes stabilization algorithms for a given Hamiltonian."""
-    u, s, vh = np.linalg.svd(t)
+    u, s, vh = la.svd(t)
     u, v = u * np.sqrt(s), vh.T.conj() * np.sqrt(s)
     prop_vecs = []
     evan_vecs = []
diff --git a/kwant/qsymm.py b/kwant/qsymm.py
index 93f4bdc3de30600c93f7c69ee6a13d5e26690630..bb752e19314a18cbc7ec53e53eb163df5bf17f80 100644
--- a/kwant/qsymm.py
+++ b/kwant/qsymm.py
@@ -448,7 +448,7 @@ def find_builder_symmetries(builder, momenta=None, params=None,
         params = dict()
 
     ham = builder_to_model(builder, momenta=momenta,
-                           real_space=False, params=params)
+                           real_space=spatial_symmetries, params=params)
 
     # Use dense or sparse computation depending on Hamiltonian size
     if sparse is None:
diff --git a/kwant/tests/test_qsymm.py b/kwant/tests/test_qsymm.py
index 2a13b5e4983068ecfc53cbd8b387fe666971b799..b21803db7a89743fb7badbbd32bb08a093f65d01 100644
--- a/kwant/tests/test_qsymm.py
+++ b/kwant/tests/test_qsymm.py
@@ -427,6 +427,34 @@ def test_find_builder_discrete_symmetries():
             assert sym_dict[class_symmetry] in builder_symmetries_dense
 
 
+def test_real_space_basis():
+    lat = kwant.lattice.honeycomb(norbs=[1, 1])
+    sym = kwant.TranslationalSymmetry(lat.vec((1, 0)), lat.vec((0, 1)))
+    bulk = kwant.Builder(sym)
+    bulk[[lat.a(0, 0), lat.b(0, 0)]] = 0
+    bulk[lat.neighbors()] = 1
+
+    # Including real space symmetries
+    symmetries = find_builder_symmetries(bulk)
+    hex_group_2D = hexagonal()
+    hex_group_2D = set(PointGroupElement(np.array(s.R).astype(float),
+                                         s.conjugate, s.antisymmetry, None)
+                for s in hex_group_2D)
+    assert len(symmetries) == len(hex_group_2D)
+    assert all([s1 in symmetries and s2 in hex_group_2D
+                for s1, s2 in zip(hex_group_2D, symmetries)])
+
+    # Only onsite discrete symmetries
+    symmetries = find_builder_symmetries(bulk, spatial_symmetries=False)
+    onsites = [PointGroupElement(np.eye(2), True, False, None),  # T
+               PointGroupElement(np.eye(2), True, True, None),   # P
+               PointGroupElement(np.eye(2), False, True, None),  # C
+               PointGroupElement(np.eye(2), False, False, None)] # I
+    assert len(symmetries) == len(onsites)
+    assert all([s1 in symmetries and s2 in onsites
+                for s1, s2 in zip(onsites, symmetries)])
+
+
 def random_onsite_hop(n, rng=0):
     rng = ensure_rng(rng)
     onsite = rng.randn(n, n) + 1j * rng.randn(n, n)