The free energetics of water density fluctuations near a surface, and the rare low-density fluctuations in particular, serve as reliable indicators of surface hydrophobicity; the easier it is to displace the interfacial waters, the more hydrophobic the underlying surface. However, characterizing the free energetics of such rare fluctuations requires computationally expensive, non-Boltzmann sampling methods like umbrella sampling. This inherent computational expense associated with umbrella sampling makes it challenging to investigate the role of polarizability or electronic structure effects in influencing interfacial fluctuations. Importantly, it also limits the size of the volume, which can be used to probe interfacial fluctuations. The latter can be particularly important in characterizing the hydrophobicity of large surfaces with molecular-level heterogeneities, such as those presented by proteins. To overcome these challenges, here we present a method for the sparse sampling of water density fluctuations, which is roughly two orders of magnitude more efficient than umbrella sampling. We employ thermodynamic integration to estimate the free energy differences between biased ensembles, thereby circumventing the umbrella sampling requirement of overlap between adjacent biased distributions. Further, a judicious choice of the biasing potential allows such free energy differences to be estimated using short simulations, so that the free energetics of water density fluctuations are obtained using only a few, short simulations. Leveraging the efficiency of the method, we characterize water density fluctuations in the entire hydration shell of the protein, ubiquitin; a large volume containing an average of more than six hundred waters.