Seismic sources are commonly idealized as point-sources due to their small spatial extent relative to seismic wavelengths. The acoustic isotropic point-radiator is inadequate as a model of seismic wave generation for seismic sources that are known to exhibit directivity. Therefore, accurate modeling of seismic wavefields must include source representations generating anisotropic radiation patterns. Such seismic sources can be modeled as linear combinations of multipole point-sources. In this paper we present a method for discretizing multipole sources in a finite difference setting, an extension of the moment matching conditions developed for the Dirac delta function in other applications. We also provide the necessary analysis and numerical evidence to demonstrate the accuracy of our singular source approximations. In particular, we develop a weak convergence theory for the discretization of a family of symmetric hyperbolic systems of first-order partial differential equations, with singular source terms, solved via staggered-grid finite difference methods. Numerical experiments demonstrate a stronger result than what is presented in our convergence theory, namely, optimal convergence rates of numerical solutions are achieved point-wise in space away from the source if an appropriate source discretization is used.