We present a computational protocol for modeling valence photoemission spectra of liquids. We use water as an experimentally well-characterized model system, and we represent its liquid state by larger finite-sized droplets. The photoemission spectrum is evaluated for an ensemble of structures along molecular dynamics simulations. The nuclear quantum effects are accounted for by ab initio based path-integral molecular dynamics simulations that are greatly accelerated with the so-called colored noise thermostat (PI+GLE) method. The ionization energies for the valence electrons are evaluated as orbital energies of optimally tuned range-separated hybrid functionals (OT-RSH). This approach provides Koopmans-type ionization energies including relaxation energy. We show that the present protocol can quantitatively describe the valence photoemission spectrum of liquid water, i.e., the positions, shapes, and widths of the photoemission peaks. With the PI+GLE simulations, even the subtle isotope effects that have been recently observed experimentally can be modeled. The electronic properties of finite-sized droplets are shown to converge rapidly to those of liquids. We discuss the importance of proper tuning of the range-separation parameter in OT-RSH as well as possible sources of error in our simulations. The present approach seems to be a viable route to modeling photoemission spectra of liquids, especially in conjunction with efficient implementation of density functional methods on graphical processing units.