Certain patterns of symmetry fractionalization in (2+1)-dimensional [(2+1)D] topologically ordered phases of matter can be anomalous, which means that they possess an obstruction to being realized in purely (2+1)D. In this paper we demonstrate how to compute the anomaly for symmetry-enriched topological states of bosons in complete generality. We demonstrate how, given any unitary modular tensor category (UMTC) and symmetry fractionalization class for a global symmetry group G, one can define a (3+1)-dimensional [(3+1)D] topologically invariant path integral in terms of a state sum for a G-symmetry-protected topological (SPT) state. We present an exactly solvable Hamiltonian for the system and demonstrate explicitly a (2+1)D G-symmetric surface termination that hosts deconfined anyon excitations described by the given UMTC and symmetry fractionalization class. We present concrete algorithms that can be used to compute anomaly indicators in general. Our approach applies to general symmetry groups, including anyon-permuting and antiunitary symmetries. In addition to providing a general way to compute the anomaly, our result also shows, by explicit construction, that every symmetry fractionalization class for any UMTC can be realized at the surface of a (3+1)D SPT state. As a by-product, this construction also provides a way of explicitly seeing how the algebraic data that defines symmetry fractionalization in general arises in the context of exactly solvable models. In the case of unitary orientation-preserving symmetries, our results can also be viewed as providing a method to compute the H-4(G, U(1)) obstruction that arises in the theory of G-crossed braided tensor categories, for which no general method has been presented to date.