diff --git a/doc/source/pre/whatsnew/0.2.rst b/doc/source/pre/whatsnew/0.2.rst
index 3a0e8e95a4a30a86168855fbd017be1279193cfe..2022810f63e0c00d3df442d91ff1e772f45ba0d8 100644
--- a/doc/source/pre/whatsnew/0.2.rst
+++ b/doc/source/pre/whatsnew/0.2.rst
@@ -42,7 +42,7 @@ description check `~kwant.plotter.plot` documentation.
 
 The behavior of `~kwant.plotter.plot` with low level systems has changed.
 Arguments of plot which are functions are given site numbers in place of
-`~kwant.builder.Site` objects when plotting a low level system.  This
+`~kwant.system.Site` objects when plotting a low level system.  This
 provides an easy way to make the appearance of lines and symbols depend on
 computation results.
 
diff --git a/doc/source/pre/whatsnew/1.0.rst b/doc/source/pre/whatsnew/1.0.rst
index 1eb4cf07efa8cda9b030497893b99b263afbabd7..473e414af1052f24ef271e38a2bd2b92ba34d361 100644
--- a/doc/source/pre/whatsnew/1.0.rst
+++ b/doc/source/pre/whatsnew/1.0.rst
@@ -39,7 +39,7 @@ now it suffices to write ::
     syst[lat.neighbors()] = t
 
 This is possible because `~kwant.builder.Builder` now accepts *functions* as
-keys in addition to `~kwant.builder.Site` objects and tuples of them
+keys in addition to `~kwant.system.Site` objects and tuples of them
 (hoppings).  These functions are expected to yield either sites or hoppings,
 when given a builder instance as the sole argument. The use of such keys is to
 implement sets of sites or hoppings that depend on what is already present in
diff --git a/doc/source/reference/kwant.builder.rst b/doc/source/reference/kwant.builder.rst
index 327afca9e9f44c797378b918baa0799ef30af45f..f453d40217ca545a38f2ab7744b3febd2f13786a 100644
--- a/doc/source/reference/kwant.builder.rst
+++ b/doc/source/reference/kwant.builder.rst
@@ -9,7 +9,6 @@ Types
    :toctree: generated/
 
    Builder
-   Site
    HoppingKind
    SimpleSiteFamily
    BuilderLead
@@ -23,7 +22,6 @@ Abstract base classes
 .. autosummary::
    :toctree: generated/
 
-   SiteFamily
    Symmetry
    Lead
 
diff --git a/doc/source/reference/kwant.system.rst b/doc/source/reference/kwant.system.rst
index 3b6f26274da393807c37c6c2826ce5518fe41c7d..6d6dbec6d73187f14143099cbf6ec10bda12559e 100644
--- a/doc/source/reference/kwant.system.rst
+++ b/doc/source/reference/kwant.system.rst
@@ -19,6 +19,16 @@ In practice, very often Kwant systems are finalized
 `kwant.builder.FiniteSystem` or `kwant.builder.InfiniteSystem`) and offer some
 additional attributes.
 
+Sites
+-----
+.. autosummary::
+   :toctree: generated/
+
+   Site
+   SiteFamily
+
+Systems
+--------
 .. autosummary::
    :toctree: generated/
 
diff --git a/doc/source/tutorial/first_steps.rst b/doc/source/tutorial/first_steps.rst
index 41ed21bb7fdfef0279b95f000db6d9199fa9ed63..b0760b6432e4d22003aa8bdfb9c18e5352999972 100644
--- a/doc/source/tutorial/first_steps.rst
+++ b/doc/source/tutorial/first_steps.rst
@@ -358,13 +358,13 @@ subbands that increases with energy.
      **sites** as indices. Sites themselves have a certain type, and
      belong to a **site family**. A site family is also used to convert
      something that represents a site (like a tuple) into a
-     proper `~kwant.builder.Site` object that can be used with
+     proper `~kwant.system.Site` object that can be used with
      `~kwant.builder.Builder`.
 
      In the above example, `lat` is the site family. ``lat(i, j)``
      then translates the description of a lattice site in terms of two
      integer indices (which is the natural way to do here) into
-     a proper `~kwant.builder.Site` object.
+     a proper `~kwant.system.Site` object.
 
      The concept of site families and sites allows `~kwant.builder.Builder`
      to mix arbitrary lattices and site families
diff --git a/doc/source/tutorial/operators.rst b/doc/source/tutorial/operators.rst
index ed346f178fc3582a2f8909c8bc4cc76de03cad53..1fdd27871a4f3852df6b32ea92f1b9d633bb6897 100644
--- a/doc/source/tutorial/operators.rst
+++ b/doc/source/tutorial/operators.rst
@@ -438,7 +438,7 @@ the `~kwant.operator.Current`, instead of a constant matrix:
     # evaluate the operator
     m_current = J_m(psi, params=dict(r0=25, delta=10))
 
-The function must take a `~kwant.builder.Site` as its first parameter,
+The function must take a `~kwant.system.Site` as its first parameter,
 and may optionally take other parameters (i.e. it must have the same
 signature as a Hamiltonian onsite function), and must return the square
 matrix that defines the operator we wish to calculate.
@@ -529,8 +529,8 @@ We see that the probability current is conserved across the scattering region,
 but the z-projected spin current is not due to the fact that the Hamiltonian
 does not commute with :math:`σ_z` everywhere in the scattering region.
 
-.. note:: ``where`` can also be provided as a sequence of `~kwant.builder.Site`
-          or a sequence of hoppings (i.e. pairs of `~kwant.builder.Site`),
+.. note:: ``where`` can also be provided as a sequence of `~kwant.system.Site`
+          or a sequence of hoppings (i.e. pairs of `~kwant.system.Site`),
           rather than a function.
 
 
diff --git a/doc/source/tutorial/plotting.rst b/doc/source/tutorial/plotting.rst
index bf6b51ae7e148fc6fdad8c37c5f24ed57421a079..48b6ed4796251b64aa9de10922ce911cd3aade86 100644
--- a/doc/source/tutorial/plotting.rst
+++ b/doc/source/tutorial/plotting.rst
@@ -105,7 +105,7 @@ the start end end site of hopping as arguments:
     kwant.plot(syst, site_lw=0.1, site_color=family_color, hop_lw=hopping_lw);
 
 Note that since we are using an unfinalized Builder, a ``site`` is really an
-instance of `~kwant.builder.Site`. With these adjustments we arrive at a plot
+instance of `~kwant.system.Site`. With these adjustments we arrive at a plot
 that carries the same information, but is much easier to interpret.
 
 Apart from plotting the *system* itself, `~kwant.plotter.plot` can also be used
@@ -175,7 +175,7 @@ this is interpreted as the radius of the inner circle).
 
 Finally, note that since we are dealing with a finalized system now, a site `i`
 is represented by an integer. In order to obtain the original
-`~kwant.builder.Site`, ``syst.sites[i]`` can be used.
+`~kwant.system.Site`, ``syst.sites[i]`` can be used.
 
 The way how data is presented of course influences what features of the data
 are best visible in a given plot. With `~kwant.plotter.plot` one can easily go
diff --git a/doc/source/tutorial/spin_potential_shape.rst b/doc/source/tutorial/spin_potential_shape.rst
index b87eaa7f6afced106110ca10823e48a6cb667069..23d90bf8bbf384af9d3644a34867c635c5eaf115 100644
--- a/doc/source/tutorial/spin_potential_shape.rst
+++ b/doc/source/tutorial/spin_potential_shape.rst
@@ -285,7 +285,7 @@ define the potential profile of a quantum well as:
         else:
             return 0
 
-This function takes two arguments: the first of type `~kwant.builder.Site`,
+This function takes two arguments: the first of type `~kwant.system.Site`,
 from which you can get the real-space coordinates using ``site.pos``, and the
 value of the potential as the second.  Note that in `potential` we can access
 variables `L` and `L_well` that are defined globally.
@@ -375,10 +375,10 @@ oscillatory transmission behavior through resonances in the quantum well.
 .. specialnote:: Technical details
 
   - Functions can also be used for hoppings. In this case, they take
-    two `~kwant.builder.Site`'s as arguments and then an arbitrary number
+    two `~kwant.system.Site`'s as arguments and then an arbitrary number
     of additional arguments.
 
-  - Apart from the real-space position `pos`, `~kwant.builder.Site` has also an
+  - Apart from the real-space position `pos`, `~kwant.system.Site` has also an
     attribute `tag` containing the lattice indices of the site.
 
 .. _tutorial-abring:
@@ -449,7 +449,7 @@ that returns ``True`` whenever a point is inside the shape, and
         return (r1 ** 2 < rsq < r2 ** 2)
 
 Note that this function takes a real-space position as argument (not a
-`~kwant.builder.Site`).
+`~kwant.system.Site`).
 
 We can now simply add all of the lattice points inside this shape at
 once, using the function `~kwant.lattice.Square.shape`
diff --git a/kwant/builder.py b/kwant/builder.py
index 390ca0ba3c54d30d080d67c5002b43a49b74f4bf..e28f63f4e4d0b8c941939712098576d6f7daf0a1 100644
--- a/kwant/builder.py
+++ b/kwant/builder.py
@@ -1,4 +1,4 @@
-# Copyright 2011-2016 Kwant authors.
+# Copyright 2011-2019 Kwant authors.
 #
 # This file is part of Kwant.  It is subject to the license terms in the file
 # LICENSE.rst found in the top-level directory of this distribution and at
@@ -19,6 +19,7 @@ import tinyarray as ta
 import numpy as np
 from scipy import sparse
 from . import system, graph, KwantDeprecationWarning, UserCodeError
+from .system import Site, SiteFamily
 from .linalg import lll
 from .operator import Density
 from .physics import DiscreteSymmetry, magnetic_gauge
@@ -26,182 +27,13 @@ from ._common import (ensure_isinstance, get_parameters, reraise_warnings,
                       interleave, deprecate_args, memoize)
 
 
-__all__ = ['Builder', 'Site', 'SiteFamily', 'SimpleSiteFamily', 'Symmetry',
-           'HoppingKind', 'Lead', 'BuilderLead', 'SelfEnergyLead', 'ModesLead',
-           'add_peierls_phase']
+__all__ = ['Builder', 'SimpleSiteFamily', 'Symmetry', 'HoppingKind', 'Lead',
+           'BuilderLead', 'SelfEnergyLead', 'ModesLead', 'add_peierls_phase']
 
 
-################ Sites and site families
-
-class Site(tuple):
-    """A site, member of a `SiteFamily`.
-
-    Sites are the vertices of the graph which describes the tight binding
-    system in a `Builder`.
-
-    A site is uniquely identified by its family and its tag.
-
-    Parameters
-    ----------
-    family : an instance of `SiteFamily`
-        The 'type' of the site.
-    tag : a hashable python object
-        The unique identifier of the site within the site family, typically a
-        vector of integers.
-
-    Raises
-    ------
-    ValueError
-        If `tag` is not a proper tag for `family`.
-
-    Notes
-    -----
-    For convenience, ``family(*tag)`` can be used instead of ``Site(family,
-    tag)`` to create a site.
-
-    The parameters of the constructor (see above) are stored as instance
-    variables under the same names.  Given a site ``site``, common things to
-    query are thus ``site.family``, ``site.tag``, and ``site.pos``.
-    """
-    __slots__ = ()
-
-    family = property(operator.itemgetter(0),
-                      doc="The site family to which the site belongs.")
-    tag = property(operator.itemgetter(1), doc="The tag of the site.")
-
-
-    def __new__(cls, family, tag, _i_know_what_i_do=False):
-        if _i_know_what_i_do:
-            return tuple.__new__(cls, (family, tag))
-        try:
-            tag = family.normalize_tag(tag)
-        except (TypeError, ValueError) as e:
-            msg = 'Tag {0} is not allowed for site family {1}: {2}'
-            raise type(e)(msg.format(repr(tag), repr(family), e.args[0]))
-        return tuple.__new__(cls, (family, tag))
-
-    def __repr__(self):
-        return 'Site({0}, {1})'.format(repr(self.family), repr(self.tag))
-
-    def __str__(self):
-        sf = self.family
-        return '<Site {0} of {1}>'.format(self.tag, sf.name if sf.name else sf)
-
-    def __getnewargs__(self):
-        return (self.family, self.tag, True)
-
-    @property
-    def pos(self):
-        """Real space position of the site.
-
-        This relies on ``family`` having a ``pos`` method (see `SiteFamily`).
-        """
-        return self.family.pos(self.tag)
-
+################ Site families
 
 @total_ordering
-class SiteFamily(metaclass=abc.ABCMeta):
-    """Abstract base class for site families.
-
-    Site families are the 'type' of `Site` objects.  Within a family, individual
-    sites are uniquely identified by tags.  Valid tags must be hashable Python
-    objects, further details are up to the family.
-
-    Site families must be immutable and fully defined by their initial
-    arguments.  They must inherit from this abstract base class and call its
-    __init__ function providing it with two arguments: a canonical
-    representation and a name.  The canonical representation will be returned as
-    the objects representation and must uniquely identify the site family
-    instance.  The name is a string used to distinguish otherwise identical site
-    families.  It may be empty. ``norbs`` defines the number of orbitals
-    on sites associated with this site family; it may be `None`, in which case
-    the number of orbitals is not specified.
-
-
-    All site families must define the method `normalize_tag` which brings a tag
-    to the standard format for this site family.
-
-    Site families that are intended for use with plotting should also provide a
-    method `pos(tag)`, which returns a vector with real-space coordinates of the
-    site belonging to this family with a given tag.
-
-    If the ``norbs`` of a site family are provided, and sites of this family
-    are used to populate a `~kwant.builder.Builder`, then the associated
-    Hamiltonian values must have the correct shape. That is, if a site family
-    has ``norbs = 2``, then any on-site terms for sites belonging to this
-    family should be 2x2 matrices. Similarly, any hoppings to/from sites
-    belonging to this family must have a matrix structure where there are two
-    rows/columns. This condition applies equally to Hamiltonian values that
-    are given by functions. If this condition is not satisfied, an error will
-    be raised.
-    """
-
-    def __init__(self, canonical_repr, name, norbs):
-        self.canonical_repr = canonical_repr
-        self.hash = hash(canonical_repr)
-        self.name = name
-        if norbs is None:
-            warnings.warn("Not specfying norbs is deprecated. Always specify "
-                          "norbs when creating site families.",
-                          KwantDeprecationWarning, stacklevel=3)
-        if norbs is not None:
-            if int(norbs) != norbs or norbs <= 0:
-                raise ValueError('The norbs parameter must be an integer > 0.')
-            norbs = int(norbs)
-        self.norbs = norbs
-
-    def __repr__(self):
-        return self.canonical_repr
-
-    def __str__(self):
-        if self.name:
-            msg = '<{0} site family {1}{2}>'
-        else:
-            msg = '<unnamed {0} site family{2}>'
-        orbs = ' with {0} orbitals'.format(self.norbs) if self.norbs else ''
-        return msg.format(self.__class__.__name__, self.name, orbs)
-
-    def __hash__(self):
-        return self.hash
-
-    def __eq__(self, other):
-        try:
-            return self.canonical_repr == other.canonical_repr
-        except AttributeError:
-            return False
-
-    def __ne__(self, other):
-        try:
-            return self.canonical_repr != other.canonical_repr
-        except AttributeError:
-            return True
-
-    def __lt__(self, other):
-        # If this raises an AttributeError, we were trying
-        # to compare it to something non-comparable anyway.
-        return self.canonical_repr < other.canonical_repr
-
-    @abc.abstractmethod
-    def normalize_tag(self, tag):
-        """Return a normalized version of the tag.
-
-        Raises TypeError or ValueError if the tag is not acceptable.
-        """
-        pass
-
-    def __call__(self, *tag):
-        """
-        A convenience function.
-
-        This function allows to write fam(1, 2) instead of Site(fam, (1, 2)).
-        """
-        # Catch a likely and difficult to find mistake.
-        if tag and isinstance(tag[0], tuple):
-            raise ValueError('Use site_family(1, 2) instead of '
-                             'site_family((1, 2))!')
-        return Site(self, tag)
-
-
 class SimpleSiteFamily(SiteFamily):
     """A site family used as an example and for testing.
 
@@ -412,8 +244,8 @@ class HoppingKind(tuple):
     ----------
     delta : Sequence of integers
         The sequence is interpreted as a vector with integer elements.
-    family_a : `~kwant.builder.SiteFamily`
-    family_b : `~kwant.builder.SiteFamily` or ``None`` (default)
+    family_a : `~kwant.system.SiteFamily`
+    family_b : `~kwant.system.SiteFamily` or ``None`` (default)
         The default value means: use the same family as `family_a`.
 
     Notes
@@ -570,7 +402,7 @@ class BuilderLead(Lead):
         The tight-binding system of a lead. It has to possess appropriate
         symmetry, and it may not contain hoppings between further than
         neighboring images of the fundamental domain.
-    interface : sequence of `Site` instances
+    interface : sequence of `~kwant.system.Site` instances
         Sequence of sites in the scattering region to which the lead is
         attached.
 
@@ -578,9 +410,9 @@ class BuilderLead(Lead):
     ----------
     builder : `Builder`
         The tight-binding system of a lead.
-    interface : list of `Site` instances
+    interface : list of `~kwant.system.Site` instances
         A sorted list of interface sites.
-    padding : list of `Site` instances
+    padding : list of `~kwant.system.Site` instances
         A sorted list of sites that originate from the lead, have the same
         onsite Hamiltonian, and are connected by the same hoppings as the lead
         sites.
@@ -638,7 +470,7 @@ class SelfEnergyLead(Lead):
     selfenergy_func : function
         Has the same signature as `selfenergy` (without the ``self``
         parameter) and returns the self energy matrix for the interface sites.
-    interface : sequence of `Site` instances
+    interface : sequence of `~kwant.system.Site` instances
     parameters : sequence of strings
         The parameters on which the lead depends.
     """
@@ -669,7 +501,7 @@ class ModesLead(Lead):
         and returns the modes of the lead as a tuple of
         `~kwant.physics.PropagatingModes` and `~kwant.physics.StabilizedModes`.
     interface :
-        sequence of `Site` instances
+        sequence of `~kwant.system.Site` instances
     parameters : sequence of strings
         The parameters on which the lead depends.
     """
@@ -797,10 +629,11 @@ class Builder:
     This is one of the central types in Kwant.  It is used to construct tight
     binding systems in a flexible way.
 
-    The nodes of the graph are `Site` instances.  The edges, i.e. the hoppings,
-    are pairs (2-tuples) of sites.  Each node and each edge has a value
-    associated with it.  The values associated with nodes are interpreted as
-    on-site Hamiltonians, the ones associated with edges as hopping integrals.
+    The nodes of the graph are `~kwant.system.Site` instances.  The edges,
+    i.e. the hoppings, are pairs (2-tuples) of sites.  Each node and each
+    edge has a value associated with it.  The values associated with nodes
+    are interpreted as on-site Hamiltonians, the ones associated with edges
+    as hopping integrals.
 
     To make the graph accessible in a way that is natural within the Python
     language it is exposed as a *mapping* (much like a built-in Python
@@ -1452,7 +1285,7 @@ class Builder:
             A boolean function of site returning whether the site should be
             added to the target builder or not. The shape must be compatible
             with the symmetry of the target builder.
-        start : `Site` instance or iterable thereof or iterable of numbers
+        start : `~kwant.system.Site` or iterable thereof or iterable of numbers
             The site(s) at which the the flood-fill starts.  If start is an
             iterable of numbers, the starting site will be
             ``template.closest(start)``.
@@ -1462,7 +1295,7 @@ class Builder:
 
         Returns
         -------
-        added_sites : list of `Site` objects that were added to the system.
+        added_sites : list of `~kwant.system.Site` that were added to the system.
 
         Raises
         ------
@@ -1614,7 +1447,7 @@ class Builder:
         ----------
         lead_builder : `Builder` with 1D translational symmetry
             Builder of the lead which has to be attached.
-        origin : `Site`
+        origin : `~kwant.system.Site`
             The site which should belong to a domain where the lead should
             begin. It is used to attach a lead inside the system, e.g. to an
             inner radius of a ring.
@@ -1624,7 +1457,7 @@ class Builder:
 
         Returns
         -------
-        added_sites : list of `Site` objects that were added to the system.
+        added_sites : list of `~kwant.system.Site`
 
         Raises
         ------
@@ -2155,11 +1988,11 @@ class FiniteSystem(_FinalizedBuilderMixin, system.FiniteSystem):
     Attributes
     ----------
     sites : sequence
-        ``sites[i]`` is the `~kwant.builder.Site` instance that corresponds
+        ``sites[i]`` is the `~kwant.system.Site` instance that corresponds
         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 high-level `~kwant.builder.Site`
+    id_by_site : mapping
+        The inverse of ``sites``; maps high-level `~kwant.system.Site`
         instances to their integer label.
         Satisfies ``id_by_site[sites[i]] == i``.
     """
@@ -2275,10 +2108,10 @@ class InfiniteSystem(_FinalizedBuilderMixin, system.InfiniteSystem):
     Attributes
     ----------
     sites : sequence
-        ``sites[i]`` is the `~kwant.builder.Site` instance that corresponds
+        ``sites[i]`` is the `~kwant.system.Site` instance that corresponds
         to the integer-labeled site ``i`` of the low-level system.
-    id_by_site : dict
-        The inverse of ``sites``; maps high-level `~kwant.builder.Site`
+    id_by_site : mapping
+        The inverse of ``sites``; maps high-level `~kwant.system.Site`
         instances to their integer label.
         Satisfies ``id_by_site[sites[i]] == i``.
 
diff --git a/kwant/continuum/discretizer.py b/kwant/continuum/discretizer.py
index 53f2379cd53f17588dc9e59cb0300a11e90c54c3..ce9e6fe2f65bffeed3da826b02f1618e474a1757 100644
--- a/kwant/continuum/discretizer.py
+++ b/kwant/continuum/discretizer.py
@@ -1,4 +1,4 @@
-# Copyright 2011-2017 Kwant authors.
+# Copyright 2011-2019 Kwant authors.
 #
 # This file is part of Kwant.  It is subject to the license terms in the file
 # LICENSE.rst found in the top-level directory of this distribution and at
@@ -20,7 +20,7 @@ from sympy.printing.lambdarepr import LambdaPrinter
 from sympy.printing.precedence import precedence
 from sympy.core.function import AppliedUndef
 
-from .. import builder, lattice
+from .. import builder, lattice, system
 from .. import KwantDeprecationWarning
 from .._common import reraise_warnings
 from ._common import (sympify, gcd, position_operators, momentum_operators,
@@ -63,7 +63,7 @@ class _DiscretizedBuilder(builder.Builder):
 
         for key, val in itertools.chain(self.site_value_pairs(),
                                         self.hopping_value_pairs()):
-            if isinstance(key, builder.Site):
+            if isinstance(key, system.Site):
                 result.append("# Onsite element:\n")
             else:
                 a, b = key
diff --git a/kwant/continuum/landau_levels.py b/kwant/continuum/landau_levels.py
index 13e70c4e6b7224d3dab456d5a4497587514ea8af..c653c55ef2634cd13e932f0c3f9c5fd026c10ff9 100644
--- a/kwant/continuum/landau_levels.py
+++ b/kwant/continuum/landau_levels.py
@@ -176,7 +176,7 @@ class LandauLattice(kwant.lattice.Monatomic):
     """
     A `~kwant.lattice.Monatomic` lattice with a Landau level index per site.
 
-    Site tags (see `~kwant.builder.SiteFamily`) are pairs of integers, where
+    Site tags (see `~kwant.system.SiteFamily`) are pairs of integers, where
     the first integer describes the real space position and the second the
     Landau level index.
 
diff --git a/kwant/kpm.py b/kwant/kpm.py
index 5da09b097e394bc57c07c7de97fb958a3665fc3b..44cf9758f4ddfff3f098be074ca929689c562f11 100644
--- a/kwant/kpm.py
+++ b/kwant/kpm.py
@@ -1,5 +1,5 @@
 # -*- coding: utf-8 -*-
-# Copyright 2011-2016 Kwant authors.
+# Copyright 2011-2019 Kwant authors.
 #
 # This file is part of Kwant.  It is subject to the license terms in the file
 # LICENSE.rst found in the top-level directory of this distribution and at
@@ -928,7 +928,7 @@ def RandomVectors(syst, where=None, rng=None):
     ----------
     syst : `~kwant.system.FiniteSystem` or matrix Hamiltonian
         If a system is passed, it should contain no leads.
-    where : sequence of `int` or `~kwant.builder.Site`, or callable, optional
+    where : sequence of `int` or `~kwant.system.Site`, or callable, optional
         Spatial range of the vectors produced. If ``syst`` is a
         `~kwant.system.FiniteSystem`, where behaves as in
         `~kwant.operator.Density`. If ``syst`` is a matrix, ``where``
@@ -950,7 +950,7 @@ class LocalVectors:
     ----------
     syst : `~kwant.system.FiniteSystem` or matrix Hamiltonian
         If a system is passed, it should contain no leads.
-    where : sequence of `int` or `~kwant.builder.Site`, or callable, optional
+    where : sequence of `int` or `~kwant.system.Site`, or callable, optional
         Spatial range of the vectors produced. If ``syst`` is a
         `~kwant.system.FiniteSystem`, where behaves as in
         `~kwant.operator.Density`. If ``syst`` is a matrix, ``where``
diff --git a/kwant/lattice.py b/kwant/lattice.py
index 899ce5474216421b97f35b922a51a62cd85b9602..4099890a480bc9f611148fe4e0e4f519e906e2fd 100644
--- a/kwant/lattice.py
+++ b/kwant/lattice.py
@@ -13,7 +13,7 @@ from math import sqrt
 from itertools import product
 import numpy as np
 import tinyarray as ta
-from . import builder
+from . import builder, system
 from .linalg import lll
 from ._common import ensure_isinstance
 
@@ -70,7 +70,7 @@ class Polyatomic:
     A Bravais lattice with an arbitrary number of sites in the basis.
 
     Contains `Monatomic` sublattices.  Note that an instance of ``Polyatomic`` is
-    not itself a `~kwant.builder.SiteFamily`, only its sublattices are.
+    not itself a `~kwant.system.SiteFamily`, only its sublattices are.
 
     Parameters
     ----------
@@ -171,7 +171,7 @@ class Polyatomic:
         >>> syst[lat.neighbors()] = 1
         """
         def shape_sites(symmetry=None):
-            Site = builder.Site
+            Site = system.Site
 
             if symmetry is None:
                 symmetry = builder.NoSymmetry()
@@ -306,7 +306,7 @@ class Polyatomic:
         """
         # This algorithm is not designed to be fast and can be improved,
         # however there is no real need.
-        Site = builder.Site
+        Site = system.Site
         sls = self.sublattices
         shortest_hopping = sls[0].n_closest(
             sls[0].pos(([0] * sls[0].lattice_dim)), 2)[-1]
@@ -396,12 +396,12 @@ def short_array_str(array):
     return full[1 : -1]
 
 
-class Monatomic(builder.SiteFamily, Polyatomic):
+class Monatomic(system.SiteFamily, Polyatomic):
     """
     A Bravais lattice with a single site in the basis.
 
-    Instances of this class provide the `~kwant.builder.SiteFamily` interface.
-    Site tags (see `~kwant.builder.SiteFamily`) are sequences of integers and
+    Instances of this class provide the `~kwant.system.SiteFamily` interface.
+    Site tags (see `~kwant.system.SiteFamily`) are sequences of integers and
     describe the lattice coordinates of a site.
 
     ``Monatomic`` instances are used as site families on their own or as
diff --git a/kwant/operator.pyx b/kwant/operator.pyx
index a10624d1ccf3d7d6a2a175e16d6e00e3126c3bad..5133ebb522994d288ef2c5344cb9d016366f6fa9 100644
--- a/kwant/operator.pyx
+++ b/kwant/operator.pyx
@@ -13,7 +13,6 @@ import cython
 from operator import itemgetter
 import functools as ft
 import collections
-import numbers
 
 import numpy as np
 import tinyarray as ta
@@ -25,7 +24,7 @@ from .graph.core cimport EdgeIterator
 from .graph.core import DisabledFeatureError, NodeDoesNotExistError
 from .graph.defs cimport gint
 from .graph.defs import gint_dtype
-from .system import InfiniteSystem
+from .system import InfiniteSystem, Site
 from ._common import UserCodeError, get_parameters, deprecate_args
 
 
@@ -179,8 +178,7 @@ def _normalize_site_where(syst, where):
             else:
                 where = [s for s in range(syst.graph.num_nodes) if where(s)]
     else:
-        # Cannot check for builder.Site due to circular imports
-        if not isinstance(where[0], numbers.Integral):
+        if isinstance(where[0], Site):
             try:
                 where = [syst.id_by_site[s] for s in where]
             except AttributeError:
@@ -223,8 +221,7 @@ def _normalize_hopping_where(syst, where):
         else:
             where = list(filter(lambda h: where(*h), syst.graph))
     else:
-        # Cannot check for builder.Site due to circular imports
-        if not isinstance(where[0][0], numbers.Integral):
+        if isinstance(where[0][0], Site):
             try:
                 where = list((syst.id_by_site[a], syst.id_by_site[b])
                                for a, b in where)
@@ -735,10 +732,10 @@ cdef class Density(_LocalOperator):
         maps from site families to square matrices. If a function is given it
         must take the same arguments as the onsite Hamiltonian functions of the
         system.
-    where : sequence of `int` or `~kwant.builder.Site`, or callable, optional
+    where : sequence of `int` or `~kwant.system.Site`, or callable, optional
         Where to evaluate the operator. If ``syst`` is not a finalized Builder,
         then this should be a sequence of integers. If a function is provided,
-        it should take a single `int` or `~kwant.builder.Site` (if ``syst`` is
+        it should take a single `int` or `~kwant.system.Site` (if ``syst`` is
         a finalized builder) and return True or False.  If not provided, the
         operator will be calculated over all sites in the system.
     check_hermiticity: bool
@@ -884,11 +881,11 @@ cdef class Current(_LocalOperator):
         matrices (scalars are allowed if the site family has 1 orbital per
         site). If a function is given it must take the same arguments as the
         onsite Hamiltonian functions of the system.
-    where : sequence of pairs of `int` or `~kwant.builder.Site`, or callable, optional
+    where : sequence of pairs of `int` or `~kwant.system.Site`, or callable, optional
         Where to evaluate the operator. If ``syst`` is not a finalized Builder,
         then this should be a sequence of pairs of integers. If a function is
         provided, it should take a pair of integers or a pair of
-        `~kwant.builder.Site` (if ``syst`` is a finalized builder) and return
+        `~kwant.system.Site` (if ``syst`` is a finalized builder) and return
         True or False.  If not provided, the operator will be calculated over
         all hoppings in the system.
     check_hermiticity : bool
@@ -1016,10 +1013,10 @@ cdef class Source(_LocalOperator):
         matrices (scalars are allowed if the site family has 1 orbital per
         site). If a function is given it must take the same arguments as the
         onsite Hamiltonian functions of the system.
-    where : sequence of `int` or `~kwant.builder.Site`, or callable, optional
+    where : sequence of `int` or `~kwant.system.Site`, or callable, optional
         Where to evaluate the operator. If ``syst`` is not a finalized Builder,
         then this should be a sequence of integers. If a function is provided,
-        it should take a single `int` or `~kwant.builder.Site` (if ``syst`` is
+        it should take a single `int` or `~kwant.system.Site` (if ``syst`` is
         a finalized builder) and return True or False.  If not provided, the
         operator will be calculated over all sites in the system.
     check_hermiticity : bool
diff --git a/kwant/physics/gauge.py b/kwant/physics/gauge.py
index eb4e696164ac9c7c3f9e5ce1c986df601d4d0bb2..cfb927a6fb1e3397ad4d739160f0e4031b736060 100644
--- a/kwant/physics/gauge.py
+++ b/kwant/physics/gauge.py
@@ -1,4 +1,4 @@
-# Copyright 2011-2018 Kwant authors.
+# Copyright 2011-2019 Kwant authors.
 #
 # This file is part of Kwant.  It is subject to the license terms in the file
 # LICENSE.rst found in the top-level directory of this distribution and at
@@ -303,7 +303,7 @@ def _loops_in_infinite(syst):
     extended_sites : callable : int -> Site
         Given a site index in the extended system consisting of
         two unit cells, returns the associated high-level
-        `kwant.builder.Site`.
+        `kwant.system.Site`.
     """
     assert isinstance(syst, system.InfiniteSystem)
     _check_infinite_syst(syst)
@@ -375,7 +375,7 @@ def _loops_in_composite(syst):
         -1 if the site is part of the reduced scattering region (see notes).
     extended_sites : callable : int -> Site
         Given a site index in the extended scattering region (see notes),
-        returns the associated high-level `kwant.builder.Site`.
+        returns the associated high-level `kwant.system.Site`.
 
     Notes
     -----
@@ -401,7 +401,7 @@ def _loops_in_composite(syst):
     # Get distance matrix for the extended scattering region,
     # a function that maps sites to their lead patches (-1 for sites
     # in the reduced scattering region), and a function that maps sites
-    # to high-level 'kwant.builder.Site' objects.
+    # to high-level 'kwant.system.Site' objects.
     distance_matrix, which_patch, extended_sites =\
         _extended_scattering_region(syst)
 
@@ -439,7 +439,7 @@ def _extended_scattering_region(syst):
         -1 if the site is part of the reduced scattering region.
     extended_sites : callable : int -> Site
         Given a site index in the extended scattering region, returns
-        the associated high-level `kwant.builder.Site`.
+        the associated high-level `kwant.system.Site`.
 
     Notes
     -----
@@ -1027,7 +1027,7 @@ class magnetic_gauge:
             The first callable computes the Peierls phase in the scattering
             region and the remaining callables compute the Peierls phases
             in each of the leads. Each callable takes a pair of
-            `~kwant.builder.Site` (a hopping) and returns a unit complex
+            `~kwant.system.Site` (a hopping) and returns a unit complex
             number (Peierls phase) that multiplies that hopping.
         """
         return self._peierls(syst_field, *lead_fields, tol=tol, average=False)
diff --git a/kwant/plotter.py b/kwant/plotter.py
index fabe5f61c9234b2323e0c11c860962b3b5b2bd5d..a9e5ec42793fb37eee4e5908624443cafc1db59d 100644
--- a/kwant/plotter.py
+++ b/kwant/plotter.py
@@ -1,5 +1,5 @@
 # -*- coding: utf-8 -*-
-# Copyright 2011-2018 Kwant authors.
+# Copyright 2011-2019 Kwant authors.
 #
 # This file is part of Kwant.  It is subject to the license terms in the file
 # LICENSE.rst found in the top-level directory of this distribution and at
@@ -402,7 +402,7 @@ def sys_leads_sites(sys, num_lead_cells=2):
     Returns
     -------
     sites : list of (site, lead_number, copy_number) tuples
-        A site is a `~kwant.builder.Site` instance if the system is not finalized,
+        A site is a `~kwant.system.Site` instance if the system is not finalized,
         and an integer otherwise.  For system sites `lead_number` is `None` and
         `copy_number` is `0`, for leads both are integers.
     lead_cells : list of slices
@@ -528,7 +528,7 @@ def sys_leads_hoppings(sys, num_lead_cells=2):
     Returns
     -------
     hoppings : list of (hopping, lead_number, copy_number) tuples
-        A site is a `~kwant.builder.Site` instance if the system is not finalized,
+        A site is a `~kwant.system.Site` instance if the system is not finalized,
         and an integer otherwise.  For system sites `lead_number` is `None` and
         `copy_number` is `0`, for leads both are integers.
     lead_cells : list of slices
@@ -812,7 +812,7 @@ def plot(sys, num_lead_cells=2, unit='nn',
 
     - The meaning of "site" depends on whether the system to be plotted is a
       builder or a low level system.  For builders, a site is a
-      kwant.builder.Site object.  For low level systems, a site is an integer
+      kwant.system.Site object.  For low level systems, a site is an integer
       -- the site number.
 
     - color and symbol definitions may be tuples, but not lists or arrays.
diff --git a/kwant/system.py b/kwant/system.py
index d270691508189866a0d75809917048303443c883..5498ead7987f4674f89eb2fcaf35523156331718 100644
--- a/kwant/system.py
+++ b/kwant/system.py
@@ -1,4 +1,4 @@
-# Copyright 2011-2013 Kwant authors.
+# Copyright 2011-2019 Kwant authors.
 #
 # This file is part of Kwant.  It is subject to the license terms in the file
 # LICENSE.rst found in the top-level directory of this distribution and at
@@ -8,13 +8,196 @@
 
 """Low-level interface of systems"""
 
-__all__ = ['System', 'FiniteSystem', 'InfiniteSystem']
+__all__ = [
+    'Site', 'SiteArray', 'SiteFamily',
+    'System', 'FiniteSystem', 'InfiniteSystem',
+]
 
 import abc
 import warnings
+import operator
 from copy import copy
+from collections import namedtuple
+from functools import total_ordering
+import numpy as np
 from . import _system
-from ._common  import deprecate_args
+from ._common  import deprecate_args, KwantDeprecationWarning
+
+
+
+################ Sites and Site families
+
+class Site(tuple):
+    """A site, member of a `SiteFamily`.
+
+    Sites are the vertices of the graph which describes the tight binding
+    system in a `Builder`.
+
+    A site is uniquely identified by its family and its tag.
+
+    Parameters
+    ----------
+    family : an instance of `SiteFamily`
+        The 'type' of the site.
+    tag : a hashable python object
+        The unique identifier of the site within the site family, typically a
+        vector of integers.
+
+    Raises
+    ------
+    ValueError
+        If `tag` is not a proper tag for `family`.
+
+    Notes
+    -----
+    For convenience, ``family(*tag)`` can be used instead of ``Site(family,
+    tag)`` to create a site.
+
+    The parameters of the constructor (see above) are stored as instance
+    variables under the same names.  Given a site ``site``, common things to
+    query are thus ``site.family``, ``site.tag``, and ``site.pos``.
+    """
+    __slots__ = ()
+
+    family = property(operator.itemgetter(0),
+                      doc="The site family to which the site belongs.")
+    tag = property(operator.itemgetter(1), doc="The tag of the site.")
+
+
+    def __new__(cls, family, tag, _i_know_what_i_do=False):
+        if _i_know_what_i_do:
+            return tuple.__new__(cls, (family, tag))
+        try:
+            tag = family.normalize_tag(tag)
+        except (TypeError, ValueError) as e:
+            msg = 'Tag {0} is not allowed for site family {1}: {2}'
+            raise type(e)(msg.format(repr(tag), repr(family), e.args[0]))
+        return tuple.__new__(cls, (family, tag))
+
+    def __repr__(self):
+        return 'Site({0}, {1})'.format(repr(self.family), repr(self.tag))
+
+    def __str__(self):
+        sf = self.family
+        return '<Site {0} of {1}>'.format(self.tag, sf.name if sf.name else sf)
+
+    def __getnewargs__(self):
+        return (self.family, self.tag, True)
+
+    @property
+    def pos(self):
+        """Real space position of the site.
+
+        This relies on ``family`` having a ``pos`` method (see `SiteFamily`).
+        """
+        return self.family.pos(self.tag)
+
+
+@total_ordering
+class SiteFamily(metaclass=abc.ABCMeta):
+    """Abstract base class for site families.
+
+    Site families are the 'type' of `Site` objects.  Within a family, individual
+    sites are uniquely identified by tags.  Valid tags must be hashable Python
+    objects, further details are up to the family.
+
+    Site families must be immutable and fully defined by their initial
+    arguments.  They must inherit from this abstract base class and call its
+    __init__ function providing it with two arguments: a canonical
+    representation and a name.  The canonical representation will be returned as
+    the objects representation and must uniquely identify the site family
+    instance.  The name is a string used to distinguish otherwise identical site
+    families.  It may be empty. ``norbs`` defines the number of orbitals
+    on sites associated with this site family; it may be `None`, in which case
+    the number of orbitals is not specified.
+
+
+    All site families must define the method `normalize_tag` which brings a tag
+    to the standard format for this site family.
+
+    Site families that are intended for use with plotting should also provide a
+    method `pos(tag)`, which returns a vector with real-space coordinates of the
+    site belonging to this family with a given tag.
+
+    If the ``norbs`` of a site family are provided, and sites of this family
+    are used to populate a `~kwant.builder.Builder`, then the associated
+    Hamiltonian values must have the correct shape. That is, if a site family
+    has ``norbs = 2``, then any on-site terms for sites belonging to this
+    family should be 2x2 matrices. Similarly, any hoppings to/from sites
+    belonging to this family must have a matrix structure where there are two
+    rows/columns. This condition applies equally to Hamiltonian values that
+    are given by functions. If this condition is not satisfied, an error will
+    be raised.
+    """
+
+    def __init__(self, canonical_repr, name, norbs):
+        self.canonical_repr = canonical_repr
+        self.hash = hash(canonical_repr)
+        self.name = name
+        if norbs is None:
+            warnings.warn("Not specfying norbs is deprecated. Always specify "
+                          "norbs when creating site families.",
+                          KwantDeprecationWarning, stacklevel=3)
+        if norbs is not None:
+            if int(norbs) != norbs or norbs <= 0:
+                raise ValueError('The norbs parameter must be an integer > 0.')
+            norbs = int(norbs)
+        self.norbs = norbs
+
+    def __repr__(self):
+        return self.canonical_repr
+
+    def __str__(self):
+        if self.name:
+            msg = '<{0} site family {1}{2}>'
+        else:
+            msg = '<unnamed {0} site family{2}>'
+        orbs = ' with {0} orbitals'.format(self.norbs) if self.norbs else ''
+        return msg.format(self.__class__.__name__, self.name, orbs)
+
+    def __hash__(self):
+        return self.hash
+
+    def __eq__(self, other):
+        try:
+            return self.canonical_repr == other.canonical_repr
+        except AttributeError:
+            return False
+
+    def __ne__(self, other):
+        try:
+            return self.canonical_repr != other.canonical_repr
+        except AttributeError:
+            return True
+
+    def __lt__(self, other):
+        # If this raises an AttributeError, we were trying
+        # to compare it to something non-comparable anyway.
+        return self.canonical_repr < other.canonical_repr
+
+    @abc.abstractmethod
+    def normalize_tag(self, tag):
+        """Return a normalized version of the tag.
+
+        Raises TypeError or ValueError if the tag is not acceptable.
+        """
+        pass
+
+    def __call__(self, *tag):
+        """
+        A convenience function.
+
+        This function allows to write fam(1, 2) instead of Site(fam, (1, 2)).
+        """
+        # Catch a likely and difficult to find mistake.
+        if tag and isinstance(tag[0], tuple):
+            raise ValueError('Use site_family(1, 2) instead of '
+                             'site_family((1, 2))!')
+        return Site(self, tag)
+
+
+
+################ Systems
 
 
 class System(metaclass=abc.ABCMeta):
diff --git a/kwant/tests/test_builder.py b/kwant/tests/test_builder.py
index e0ed4bfe274998b684f89d167386f1911ff9c270..f41652a85c1306fbd25e6b1d5015dd1a550d50b2 100644
--- a/kwant/tests/test_builder.py
+++ b/kwant/tests/test_builder.py
@@ -19,7 +19,7 @@ from pytest import raises, warns
 from numpy.testing import assert_almost_equal
 
 import kwant
-from kwant import builder
+from kwant import builder, system
 from kwant._common import ensure_rng
 
 
@@ -421,12 +421,12 @@ def test_finalization():
 
     # Attach lead with improper interface.
     syst.leads[-1] = builder.BuilderLead(
-        lead, 2 * tuple(builder.Site(fam, n) for n in neighbors))
+        lead, 2 * tuple(system.Site(fam, n) for n in neighbors))
     raises(ValueError, syst.finalized)
 
     # Attach lead properly.
     syst.leads[-1] = builder.BuilderLead(
-        lead, (builder.Site(fam, n) for n in neighbors))
+        lead, (system.Site(fam, n) for n in neighbors))
     fsyst = syst.finalized()
     assert len(fsyst.lead_interfaces) == 1
     assert ([fsyst.sites[i].tag for i in fsyst.lead_interfaces[0]] ==
@@ -434,15 +434,15 @@ def test_finalization():
 
     # test that we cannot finalize a system with a badly sorted interface order
     raises(ValueError, builder.InfiniteSystem, lead,
-           [builder.Site(fam, n) for n in reversed(neighbors)])
+           [system.Site(fam, n) for n in reversed(neighbors)])
     # site ordering independent of whether interface was specified
-    flead_order = builder.InfiniteSystem(lead, [builder.Site(fam, n)
+    flead_order = builder.InfiniteSystem(lead, [system.Site(fam, n)
                                                 for n in neighbors])
     assert flead.sites == flead_order.sites
 
 
     syst.leads[-1] = builder.BuilderLead(
-        lead, (builder.Site(fam, n) for n in neighbors))
+        lead, (system.Site(fam, n) for n in neighbors))
     fsyst = syst.finalized()
     assert len(fsyst.lead_interfaces) == 1
     assert ([fsyst.sites[i].tag for i in fsyst.lead_interfaces[0]] ==