text: Hi,
  The finalized system has a `sites` attribute which is a sequence of
  sites ordered in the same way as in the wavefunction. You should know
  that in the wavefunction the values on degrees of freedom which belong
  to the same site are stored sequentially.
  I see you have 4 degrees of freedom per site in your model, therefore
  you can do something like the following:
  import numpy as np
  import kwant
  sys = make_system()
  fsys = sys.finalized()
  wavefunction = ....
  ...
  # get sites in "sheet"
  def in_sheet(index, site):
  return site.pos[1] == 0 # site in x-z plane
  sheet_idx_sites = filter(in_sheet, enumerate(fsys.sites))
  n = len(sheet_sites)
  # prepare storage
  projected_wavefunction = np.empty((n, 4), dtype=np.complex)
  # get wavefunction projection
  for i, (idx, site) in enumerate(sheet_idx_sites):
  projected_wavefunction[i] = wavefunction[idx*4: (idx+1)*4]
  Similar questions have been asked several times on the mailing list,
  although never explicitly with >1 orbital per site.
  Hope that helps,
  Joe
  P.S.
  I'm actually looking to implement better support for manipulating
  quantities which are defined over
  sites/hoppings, and am talking with Christoph about having this in kwant
  2.0.