Trying to define an operator using a system with no sites gives a confusing error message
To reproduce:
>>> import kwant; kwant.operator.Density(kwant.Builder().finalized(), [[1, 0], [0, 1]])
...
ValueError: Single `onsite` matrix provided, but there are different numbers of orbitals on different sites
Should we perhaps just refuse to finalize a system with no sites?