The variables involved here are classical and you resolve them classically. They enter into the operator because they are parameters of the wavefunction.
So let's do this a little more broadly. For continuous systems, we want a family of solutions based on some parameters which I'll collectively identify as $\alpha \in A$; the solutions are then labeled simply as $|\alpha\rangle$.
They do not need to be eigenfunctions of a Hamiltonian or orthogonal or any such thing to use them as a "basis", they just need to satisfy one major requirement: there must be a normalization kernel, a function $\mathcal N(\alpha)$, such that they, with the kernel, "resolve the identity":$$\hat 1 = \int_{A} d\alpha ~\mathcal N(\alpha)~ |\alpha\rangle\langle\alpha|,$$ where $\hat 1$ is the identity operator.
This generalization of the usual probability rule $|\psi(\alpha)|^2 = \operatorname{Pr}(a = \alpha|\text{state} = \psi)$is necessary in certain contexts, for example when you want to talk about coherent states as a basis for your wavefunctions and they overlap so that $\langle \alpha | \beta \rangle \ne 0$ when $\alpha \ne \beta$.
Once you have the above property, you can calculate any expectation value with those coordinates. The easiest way to see this is to use the density matrix formalism, where you have a state matrix $\rho$ which for a pure state $\Psi$ is just $|\Psi\rangle\langle\Psi|$ and which in general evolves as $i \hbar \dot\rho = [\hat H, \rho].$ In this formalism, because trace is cyclically permutative, expectation values become $$\langle A \rangle_\rho = \operatorname{Tr} \left(\hat A~\rho\right)$$And by trivially inserting the above identity we can rewrite these as:$$\langle A \rangle_\rho = \int_A d\alpha~\mathcal N(\alpha)~\langle \alpha | ~\hat A ~ \rho~ |\alpha\rangle.$$When we classically do the canonical partition function we really want to sum up $e^{-\beta H}$ over all possible states with no bias towards any one state. The proper approach here is to insert $\rho = \hat 1$ as the "diagonal" elements of a normalized $\rho$ usually give probabilities for a given state: we sum all of them with equal "probability" 1.
Therefore we find ourselves trying to compute in this formalism $$Z(\beta) = \operatorname {Tr} e^{-\beta \hat H} = \int_A d\alpha~\mathcal N(\alpha)~\langle \alpha | ~e^{-\beta \hat H}~ |\alpha\rangle.$$Now tracing back through this expression: $\alpha$ is just some collection of parameters that lives in some space $A$. There is nothing quantum about $\alpha$ per se, but there are "corresponding" bras and kets labeled with $\alpha$ which are meaningful to the quantum realm. By inspection, everything "quantum" that they do has already been done by the time we take the integral.
Now, one choice which we can choose is "inspired" by the coherent states: let $\alpha = (q, p)$, $A = \mathbb R^2$, and $\lambda$ be an arbitrary constant; we can now define $$|q, p\rangle = \int_{-\infty}^\infty dx~ \left(2 \pi \lambda^2\right)^{-1/4}~ \exp\left({(x - q)^2 \over 4 \lambda^2} + i ~ \frac{p ~ x}{\hbar}\right)~|x\rangle$$ There are no operators here but it's not too hard to see that $\langle q, p | \hat x | q, p\rangle = q$ and $\langle q, p | \hat p | q, p\rangle = p.$
If I've done everything correct then $\mathcal N(q, p) = 1 / (2 \pi \hbar)$ no matter what $\lambda$ actually is, and the result that they gave follows directly from that: $$Z = \frac{1}{h} \int_{\mathbb R^2} dq~dp~\langle q, p| e^{-\beta \hat H} |q, p\rangle$$So the nuance here is that the underlying wavefunctions take position and momentum as parameters and because position and momentum aren't independent they end up, in a lot of bases, apparently "dividing the world up" into phase-space volumes $\frac {dp~dq} h $. In this case we enforce that through this "must resolve the identity" criterion which lets us insert that basis into our expectation values.
So the answers to your questions are (1) yes, as long as your basis wavefunctions $|q, p\rangle$ are sanely defined, and (2) they're normal numbers; you integrate them normally.
