The strong $L^2$-approximation of occupation time functionals is studied with respect to discrete observations of a $d$-dimensional càdlàg process. Upper bounds on the error are obtained under weak assumptions, generalizing previous results in the literature considerably. The approach relies on regularity for the marginals of the process and applies also to non-Markovian processes, such as fractional Brownian motion. The results are used to approximate occupation times and local times. For Brownian motion, the upper bounds are shown to be sharp up to a log-factor.