If $f$ is of bounded variation, then there are bounds given by a (multi-dimensional generalization of a) theorem of Koksma. A reference is Kuipers and Niederreiter, Uniform Distribution Of Sequences.

EDIT: Here are a couple of results from that book.

Theorem 5.5: Koksma-Hlawka Inequality. Let $f(x)$ be of bounded variation on $[0,1]^k$ in the sense of Hardy and Krause. Let $\omega$ be the finite sequence of points
$${\bf x}_1,\dots,{\bf x}_N$$

in $[0,1]^k$, and let $\omega_{j_1,m\dots,j_p}$ denote the projection of the sequence $\omega$ on the $k-p$-dimensional face of $[0,1]^k$ defined by $x^{(j_1)}=\cdots=x^{(j_p)}=1$. Then we have

$$
\left|{1\over N}\sum_{n=1}^Nf({\bf x_{\it n}})-\int_{[0,1]^k}f({\bf x})d{\bf x}\right|\le\sum_{p=1}^k\sum_{1,\dots,k;p}^*D_N^*(\omega_{p+1,\dots,k})V^{(p)}(f(\dots,1,\dots,1))
$$

where $V^{(p)}(f(\dots,1,\dots,1))$ denotes the $p$-dimensional variation of $f(x^{(1)},\dots,x^{(p)},1,\dots,1)$ on $[0,1]^p$ in the sense of Vitali and where the term of the sum corresponding to $p=k$ is understood to be $D_N^*(\omega)V^{(k)}(f)$.

Here $D$ is a discrepancy, probably very simple to calculate for the situation at hand, but I'm not up to typing it out. I'm going to bail on typing out Theorem 5.6, too; it applies when $f$ has certain continuous partial derivatives, and replaces the variation $V$ with an integral of the absolute value of said derivatives.