The coherent state path integral for the partition function of a theory of bosons is an integral over fields $\phi(\tau,x),\overline{\phi}(\tau,x)$; these are the classical fields corresponding to the bosonic creation annihilation operators $b(x),b^\dagger (x)$:
$\text{tr}(e^{-\beta \hat{H}}) = \int d[\phi,\overline{\phi}] e^{-S[\phi,\overline{\phi}]}$
Now perform a change of variable in the path integral $\phi(x), \overline{\phi} (x) \rightarrow n(x),\theta(x)$ where $\phi = \sqrt{n} e^{\text{i}\theta}$. This is the "number phase" representation (see Altland and Simons, chapter 6). The change of variables is a canonical mapping i.e., it preserves the classical poisson bracket. [It has some pathologies at $n=0$ which we'll ignore]
So $n,\theta$ are canonically conjugate variables. Now remember that in a theory of bosons, $\hat{b}^\dagger(x) \hat{b}(x)$, or $\overline{\phi}(x) \phi(x)$ in the path integral, is the local particle density. But note too that $\overline{\phi(x)} \phi(x) = n(x)$; hence we can identify $n(x)$ with the local particle density. What you're calling $N$ is the zeroth spatial fourier component of $n(x)$, and what you're calling $\theta_0$ is the zeroth fourier component of $\theta(x)$. That $N,\theta_0$ are canonically conjugate follows from fact $n(x),\theta(x)$ are canonically conjugate.