Does the integral of a Wigner function over a finite region mean anything? I've recently been dipping my toes deeper into the so-called "Wigner function" formalism for quantum theory, and what I am curious about is this: ostensibly, the Wigner function is the analogue of a classical probability density function, or pdf, over phase space - with the interesting property that it can take on negative values (so-called "negative probability"). Now, I have often heard that an integral of a Wigner over a finite region, i.e.
$$P[(x, p) \in R] = \iint_R W(x, p)\ dx\ dp$$
is "meaningless" presumably because you cannot simultaneously measure $x$ and $p$, due to the Heisenberg uncertainty limit. But that's actually not correct - in a naive sense, you just can't measure them to infinite precision. After all, the Heisenberg principle doesn't say you must have all or nothing, it says there is a limit, bounded by $\hbar$, and the Wigner function likewise cannot be localized to regions smaller than about $1\ \hbar$ in area in the phase space. Moreover, in the most-upvoted anwer here:
Why is it impossible to measure position and momentum at the same time with arbitrary precision?
it is mentioned how that simultaneous measurement on $x$ and $p$ is in fact possible to arbitrary "precision" (in the Wigner formalism, presumably to arbitrarily below 1 $\hbar$ of phase space), you can still get "results", it's just that they will fail to reproduce, i.e. a second measurement immediately on the heels of the first will yield another pair of random values for $x$ and $p$, which makes perfect sense if you think of the probability distribution as becoming localized but not arbitrarily so, so that it retains non-trivial spread at all times.
Thus what I am wondering is: does the above integral, then, have any relation to such simultaneous measurements? Like, at the very least, if the area of the integration region $R$ is much, much larger than $\hbar$ (say $1000\ \hbar$ or more), and we're given a valid state that is widely spread in phase space, then does the integral converge to, say, the classical probability to find the particle with those position and momentum parameters? (Hence the suggestive notation above.) Moreover, at the small end of the scale, is the given $P$ in any way relatable to the probabilities involved in what is discussed above about simultaneous measurement on $x$ and $p$ (though note here it can be negative, suggesting further interpretation is required)?
 A: @Quantum Mechanic appears to have answered your question in his comments: of course conditioned (masked) quasiprobability distributions are meaningful, just as such probability distributions are meaningful: you just restrict the domain of xs and ps you are interested in; say, your detectors are limited in their performance. Your distribution will predict the statistical  incidence/outcome of your answer, since an isolated individual measurement is meaningless: a flash in the pan.
QM is a statistical theory predicting mere probabilities of outcomes. Heisenberg's UP, $\Delta x ~\Delta p\geq \hbar/2$, restricts the variance of such probabilistic expectation values (statistical averages of measurement outcomes) and is astoundingly embedded in the fabric of the formalism automatically—almost magically to the novice.

*

*An aside: this magic suffices to prove regions of W continuously negative can be no larger than about an ℏ—I gather you are familiar with the simple argument; so, indeed, the UP  protects them from statistically meaningful measurement.

Bounds of W for precisely the subregions you are talking about (disks) have been ingeniously worked out by Bracken and associates, in a group of unjustifiably shrugged-off papers:

*

*Bracken, A. J., Doebner, H. D., & Wood, J. G. (1999), Bounds on integrals of the Wigner function, Phys Rev Lett 83(19), 3758


*Bracken, A. J., Ellinas, D., & Wood, J. G. (2004), Non-positivity of the Wigner function and bounds on associated integrals Acta Physica Hungarica Series B, Quantum Electronics, 20(1), 121-124


*Bracken, A. J., Ellinas, D., & Wood, J. G. (2003), Group theory and quasiprobability integrals of Wigner functions Journal of Physics A: Mathematical and General 36(20), L297


*Bracken, A. J., & Wood, J. G. (2004), Non-positivity of Groenewold operators, Europhysics Letters 68(1), 1
...but they don't have that much bearing on your question, since they focus on small phase-space regions, not the "dull" large ones you are talking about. (Recall, like most professionals, they  nondimensionalize ℏ into the units of x and p; your units probably measure these in units of $\sqrt \hbar$.)
The dull large ones are basically like Liouville density PDFs, as you put it, (and you may force them to be positive definite by convolving them with a Gaussian of width larger than ℏ (turning them into Hussimis) with little alteration: the analog of low-pass filtering in signal processing.
The small ones are aggressively nonclassical. They cannot serve  to give you measurement outcome probabilities, and W is bounded above and below by $\pm {1\over \hbar}$. Unlike a classical mechanics Liouville density PDF, it cannot be a δ-function at a single point in phase space:

*

*Classical mechanics is spikey and precise, while QM is fuzzy and imprecise.

Back to the large ones now. To picture the simplest configuration for a phase-space disk region of area  ~10000ℏ, (I give it another order of magnitude to have an easy square root), a coherent state inside this area
will correspond to a Schroedinger wave packet of an effective x width of about 100$\sqrt \hbar$.
Stepping away from its profile like a bird, thus shrinking the coordinates, but normalizing the probability to still be 1, you'll see the familiar (phase-averaged) classical oscillator cookie cutter in phase space. I think I recall Leonhardt's book discusses measurements of such.
By significant contrast, a pure oscillator eigenstate confined to this region will have
$$
\langle x^2+ p^2 \rangle=\hbar (2n+1) = 10^4 \hbar,
$$
so it will correspond to the 5000th excited state; it will thus have 5000 nodes specified by the 5000th Laguerre polynomial.
Flying away and rescaling as above, you of course see the wiggles blend into something smoother- looking and flattish... a  pancake with a needle in the center. Nowhere near as nice a classical limit as the above coherent state.

Added note responding to comment

is there some minimum scale  at which, if a disc has radius >,then the region inside the disc always integrates to a nonnegative number...

Yes, as stated above.  Regions of negative value in the WF are small, no more than a few  ℏs in size, as  argued in the booklet linked. Proof by contradiction: Suppose there existed a puddle in the WF with uniformly negative value larger than a few  ℏs. It is a theorem (N Cartwright, Physica 83A (1976) 210-213) that  convolving with a phase-space Gaussian of size ℏ (so a Weierstrass transform) necessarily results in a positive-semidefinite function for WFs. But one assumed the puddle was of negative values, so, in that region, the convolution of a positive Gaussian with the negative puddle would be negative! The  conclusion is that there cannot be such large negative puddles.
A: Let me define for a density matrix $\rho$ with integral kernel $\rho(x,y)$ the Wigner measure by
$$
W_{\rho}(x,p) = \frac{1}{\hbar^d}\int_{\mathbb R^3} e^{-i\,y\cdot p/\hbar}\,\rho(x+y/2,x-y/2)\,\mathrm d y
$$
The Wigner function indeed converges to a classical distribution $f(x,p)$ in the phase space when $\hbar\to 0$, called a Wigner measure (as proved mathematically in Sur les Mesures de Wigner). As it is written in in french, let me comment on that.
First, the limiting measure is the right classical distribution on the phase space. Indeed, if $\psi$ is solution of Schrödinger equation and $ρ = |\psi\rangle\langle\psi|$, or more generally if $\rho$ is a solution of Liouville-Von Neumann equation, then $f$ is a solution of the classical Liouville equation, which is the classical PDE version of Newton's equations.
However, in general, one should be a bit careful because for a single wave function (or equivalently an operator $ρ = |\psi\rangle\langle\psi|$), since the Wigner function will converge to a measure (so it might be supported on a space of measure $0$, like being a sum of dirac deltas for example) and in a very weak sense (so the integral $\int_{R} W_{\rho}\,\mathrm d x\,\mathrm d p$ might not be very close to the $\int_{R} f\,\mathrm d x\,\mathrm d p$.
There are however some good case, such as the case of mixed states. If for example the kinetic energy and the Hilbert-Schmidt norm are finite and such that
$$
\int_{\mathbb R^6}  W_{\rho}(x,p)\,|p|^2\,\mathrm d x\,\mathrm d p = \mathrm{Tr}(-\hbar^2\Delta\rho) < C
\\
\int_{\mathbb R^6}  |W_{\rho}(x,p)|^2\,\mathrm d x\,\mathrm d p = h^{-3}\,\mathrm{Tr}(|\rho|^2) < C
$$
then indeed $\int_{R} W_{\rho}\,\mathrm d x\,\mathrm d p \to \int_{R} f\,\mathrm d x\,\mathrm d p$ (or equivalently, the difference of the two integrals is small in comparison to $\hbar$).
As indicated by Cosmas Zachos, a possibility to improve the accuracy and require less hypotheses is indeed to make a convolution with a Gaussian approximating a Dirac delta, to obtain the so called Husimi transform.
So, to summarize, the difficulties occurs when the limiting distribution is singular ... and this is the case for single wavefunctions. On the contrary, in the case of mixed states, one can a have much better correspondence. For instance, the Wigner transform cannot be very negative anymore: the $1/\hbar$ becomes $\hbar$! Indeed, it is not difficult to show that if $W_\rho$ is regular, then it is close to its Husimi transform (which is positive) and obtain a bound on the form
$$
W_{\rho}(x,p) ≥ -\tfrac{3\hbar}{2} \|\nabla^2W_{\rho}\|_{L^\infty}
$$
