H2O has 10 electrons, two of which are oxygen 1s. This leaves you eight electrons. These eight electrons fit onto 4 spatial orbitals in a spin-restricted calculation for the singlet.
Determining localized orbitals  (for periodic boundary conditions these are called Wannier functions), the 4 orbitals localize into two covalent OH bonds, and two lone pairs ...
The following code is a python function to read the wannier90_hr.dat, from which you should figure out its data structure.
Read hopping matrix element from the wannier90 output file:wannier90_hr.dat
wan_num: number of wannier functions
wsc_num: number of wigner-sitez cells
You can read the wanniertools code. In wanniertools, to calculate surface state, they write a surfstate subroutine in surfstate.f90 file. The slab Hamiltonian is restructured from a bulk Hamiltonian in ham_qlayer2qlayer.f90 file.
enter code here
! This is a fortran code.
! H00 Hamiltonian between nearest neighbour-quintuple-layers
! the factor 2 is ...
Wannier90 might not be good at preserving the symmetry. But they probably include a few new methods to enforce symmetry in Wannier90.v.3.1.0. Maybe you can check this.
Also, WannierTools can symmetrize the hr.dat, but from my personal experience it sometimes gives you worse results than the original hr.dat.
Wannier90 should still scale cubically with size of the supercell, same as pw.x on Quantum ESPRESSO. Regarding the number of iterations needed, it depends on your system. Assuming you are generating MLWFs (Maximally Localized Wannier Functions) first, the number of iterations required to converge depends on the parameters in your system, like K-Points. Look ...