diff --git a/kwant/builder.py b/kwant/builder.py
index d972a7c99d0f9ea7f4b281c56e5ae8c416dfe8cb..9b6a028440af07bc7a1d0b5beeb4dc671e71fdfd 100644
--- a/kwant/builder.py
+++ b/kwant/builder.py
@@ -18,7 +18,7 @@ import operator
 from itertools import izip, islice, chain
 import tinyarray as ta
 import numpy as np
-from . import system, graph
+from . import system, graph, KwantDeprecationWarning
 
 
 class BuilderKeyError(KeyError):
@@ -1385,6 +1385,9 @@ class FiniteSystem(system.FiniteSystem):
         return value
 
     def site(self, i):
+        warnings.warn("The function `site` will disappear after Kwant 1.1.  "
+                      "Use `sites` instead.", KwantDeprecationWarning,
+                      stacklevel=2)
         return self.sites[i]
 
     def pos(self, i):
@@ -1419,6 +1422,9 @@ class InfiniteSystem(system.InfiniteSystem):
         return value
 
     def site(self, i):
+        warnings.warn("The function `site` will disappear after Kwant 1.1.  "
+                      "Use `sites` instead.", KwantDeprecationWarning,
+                      stacklevel=2)
         return self.sites[i]
 
     def pos(self, i):
diff --git a/kwant/plotter.py b/kwant/plotter.py
index 87c534526a6886abeeb353557b06b9dbd733044a..5932bc515479c272b38d1627df4312a6c20a6024 100644
--- a/kwant/plotter.py
+++ b/kwant/plotter.py
@@ -761,7 +761,7 @@ def sys_leads_pos(sys, site_lead_nr):
         else:
             try:
                 sym = sys.leads[lead_nr].symmetry
-                site = sys.site(sys.lead_interfaces[lead_nr][0])
+                site = sys.sites[sys.lead_interfaces[lead_nr][0]]
             except (AttributeError, IndexError):
                 # empty leads, or leads without symmetry aren't drawn anyways
                 return (0, 0)
@@ -907,7 +907,7 @@ def sys_leads_hopping_pos(sys, hop_lead_nr):
         else:
             try:
                 sym = sys.leads[lead_nr].symmetry
-                site = sys.site(sys.lead_interfaces[lead_nr][0])
+                site = sys.sites[sys.lead_interfaces[lead_nr][0]]
             except (AttributeError, IndexError):
                 # empyt leads or leads without symmetry are not drawn anyways
                 return (0, 0)
diff --git a/kwant/tests/test_builder.py b/kwant/tests/test_builder.py
index 420291a7a2f91d90767820554aa4e8c06a91e320..8416e154cdd4b2e15d423781e819b739da2c51e4 100644
--- a/kwant/tests/test_builder.py
+++ b/kwant/tests/test_builder.py
@@ -191,7 +191,7 @@ def random_hopping_integral(rng):
 def check_onsite(fsys, sites, subset=False, check_values=True):
     freq = {}
     for node in xrange(fsys.graph.num_nodes):
-        site = fsys.site(node).tag
+        site = fsys.sites[node].tag
         freq[site] = freq.get(site, 0) + 1
         if check_values and site in sites:
             assert fsys.onsite_hamiltonians[node] is sites[site]
@@ -208,8 +208,8 @@ def check_hoppings(fsys, hops):
     assert_equal(fsys.graph.num_edges, 2 * len(hops))
     for edge_id, edge in enumerate(fsys.graph):
         tail, head = edge
-        tail = fsys.site(tail).tag
-        head = fsys.site(head).tag
+        tail = fsys.sites[tail].tag
+        head = fsys.sites[head].tag
         value = fsys.hoppings[edge_id]
         if value is builder.Other:
             assert (head, tail) in hops
@@ -320,7 +320,7 @@ def test_finalization():
         lead, (builder.Site(fam, n) for n in neighbors))
     fsys = sys.finalized()
     assert_equal(len(fsys.lead_interfaces), 1)
-    assert_equal([fsys.site(i).tag for i in fsys.lead_interfaces[0]],
+    assert_equal([fsys.sites[i].tag for i in fsys.lead_interfaces[0]],
                  neighbors)
 
     # Add a hopping to the lead which couples two next-nearest cells and check
@@ -354,14 +354,14 @@ def test_hamiltonian_evaluation():
     assert_equal(fsys.graph.num_edges, 2 * len(edges))
 
     for i in range(len(tags)):
-        site = fsys.site(i)
+        site = fsys.sites[i]
         assert site in sites
         assert_equal(fsys.hamiltonian(i, i),
                      sys[site](site))
 
     for t, h in fsys.graph:
-        tsite = fsys.site(t)
-        hsite = fsys.site(h)
+        tsite = fsys.sites[t]
+        hsite = fsys.sites[h]
         assert_equal(fsys.hamiltonian(t, h),
                      sys[tsite, hsite](tsite, hsite))
 
@@ -600,7 +600,7 @@ def test_ModesLead_and_SelfEnergyLead():
 
     # Replace lead with it's finalized copy.
     lead = fsys.leads[1]
-    interface = [lat(L-1, lead.site(i).tag[1]) for i in xrange(L)]
+    interface = [lat(L-1, lead.sites[i].tag[1]) for i in xrange(L)]
 
     # Re-attach right lead as ModesLead.
     sys.leads[1] = builder.ModesLead(lead.modes, interface)