Consider the Grassmannian manifold $G(M, N)$ of $M$-dimensional subspaces in $R^N$. I want to approximate (stochastically) an integral of the form $$ \int_{G(M, N)} f(v) \, dv, $$ where $f : G(M, N) \to R$ is some function and $dv$ is the Haar measure on the Grassmannian. I want to approximate the integral with sampling, and therefore I need a method to uniformly draw samples with respect to the measure dv.
I'm happy about hints / references on how to do that.
To expand on one of the comments to your question. Indeed, sampling each component of an $M \times N$ matrix independently from the standard normal distribution yields a random $M$-dimensional subspace of $\mathbb{R}^N$.
Using the same reference as in this (related) question, i.e., Chikuse, Y. (2003). Statistics on Special Manifolds, we have from Theorem 2.2.2. that if the elements of $Z \in \mathbb{R}^{M \times N}$ are i.i.d. from $\mathcal{N}(0, 1)$, then $P = Z (Z^\top Z)^{-1} Z^\top$ is uniformly distributed on the subset of $\mathbb{R}^{N \times N}$ containing rank-$M$ projection matrices $P = X X^\top$, which is just another representation of the Grassmann manifold.
Thus, if you instead represent points via the column space of rank-$M$ matrices from $\mathbb{R}^{N \times M}$, you can use the matrix $Z$ directly. Note that $Z^\top Z$ has almost surely full rank, so its inverse above only amounts to a change of basis which should keep you in the same equivalence class (i.e., same Grassmann point).
Similarly, if you use orthonormal $M$-frames, then Gram-Schmidt should give you yet another representation of the same subspace.