I think it is pretty obviously in your last line. So what are the results of the commutator $[p_r,p_\theta]$ and $[p_\theta,p_\theta]$? Note that $p_\theta$ is not commute with $\theta$. Then you should get the answer.
I think that your issues are not specifically due to classical field theory, but rather Hamiltonian mechanics in general. Formally, there is no conceptual difference between the two, you just replace the Kronecker deltas of the finite dimensional phase space with Dirac deltas in field theory.
The Hamilton's equations of motion are typically obtained by a variational principle. Your real trajectories are stationary points of the action:
$$
S = \int pdq-Hdt
$$
Using calculus of variation, you can deduce Hamilton's equations of motion the Euler-Lagrange equation gives:
$$
\dot q = \frac{\partial H}{\partial p} \\
\dot p = -\frac{\partial H}{\partial q}
$$
which you can then rewrite more conveniently using Poisson brackets. Note that this is equivalent to the Euler-Lagrange equation of the corresponding Lagrangian thanks to the property of the Legendre transform, even though the candidate trajectories in the Hamiltonian formalism are more diverse than in the Lagrangian setting. Therefore, if you capture the "right" physics in your Lagrangian (typically by symmetry considerations), then the Hamilton equations of motion will also capture the "right" physics due to this equivalence.
Back to your field theory setting, the same method applies. From the variational principle using the action $S$:
$$
S = \int \Pi\cdot dq -H dt \\
\Pi \cdot \dot q = \int \Pi(x) \dot q(x) d^nx
$$
you obtain the Hamilton's equations of motion:
$$
\begin{align}
\dot q(x) &= \frac{\delta H}{\delta \Pi(x)} \\
&= \frac{\partial \mathcal H}{\partial \Pi}(x) \\
\dot\Pi(x,) &= -\frac{\delta H}{\delta q(x)} \\
&= -\frac{\partial \mathcal H}{\partial q}(x)
\end{align}
$$
where the second equations assume that (as it is often the case) the Hamiltonian is of the form:
$$
H = \int \mathcal H(q,\Pi)(x) d^nx
$$
The correct definition of the Poisson bracket two functionals $A,B$ of $q,\Pi$ is:
$$
[A,B] = \frac{\delta A}{\delta q}\frac{\delta B}{\delta \Pi}-\frac{\delta A}{\delta \Pi}\frac{\delta B}{\delta q}
$$
The Poisson bracket is still a functional of $q,\Pi$, but is implicitly a function of $x$. Things simplify if you assume that they are of the form:
$$
A = \int [\mathcal A(\Pi,q)](x) d^nx\\
B = \int [\mathcal B(\Pi,q)](x) d^nx
$$
In which case:
$$
\begin{align}
\frac{\delta A}{\delta q(x)} &= \frac{\partial \mathcal A}{\partial q}(x) & \frac{\delta B}{\delta q(x)} &= \frac{\partial\mathcal B}{\partial q}(x) \\
\frac{\delta A}{\delta \Pi(x)} &= \frac{\partial\mathcal A}{\partial \Pi}(x) & \frac{\delta\mathcal B}{\delta \Pi(x)} &= \frac{\partial\mathcal B}{\partial \Pi}(x) \\
[A,B] &= \frac{\partial \mathcal A}{\partial q}\frac{\partial \mathcal B}{\partial \Pi}-\frac{\partial \mathcal B}{\partial q}\frac{\partial \mathcal A}{\partial \Pi}
\end{align}
$$
Even the elementary fields can be written this way:
$$
q(x) = \int q(y)\delta(y-x)d^ny \\
\Pi(x) = \int \Pi(y)\delta(y-x)d^ny
$$
Or more generally, this applies to local observables as well, which justify the equations whose accuracy you were not sure:
$$
\mathcal A(x) = \int [\mathcal A(q,\Pi)] (y)\delta(y-x)d^ny \\
\mathcal B(x) = \int [\mathcal B(q,\Pi)] (y)\delta(y-x)d^ny
$$
This is why you have:
$$
\begin{align}
[q(x),\Pi(y)] &= \delta(x-y) \\
[q(x),A] &= \frac{\partial \mathcal A}{\partial \Pi}(x) \\
[\Pi(x),A] &= -\frac{\partial \mathcal A}{\partial q}(x) \\
[\mathcal A(x),\mathcal B(y)] &= \frac{\partial \mathcal A}{\partial q}(x)\frac{\partial \mathcal B}{\partial \Pi}(x)\delta(x-y)-\frac{\partial \mathcal B}{\partial q}(x)\frac{\partial \mathcal A}{\partial \Pi}(x)\delta(x-y)
\end{align}
$$
Hope this helps.
Best Answer
I can't help you with references, but there are lots of works on classical canonical transformations in classical field theory.
It is standard for one pair of conjugate fields $X(x)$ and $V(x)$ and a functional of them, $f[X,V]$, to consider
$$ \left\{ f[X,Y] , g[X,Y]\right \}= \int\!\! dx ~~~\Bigl ( \frac{\delta f}{\delta X(x)}\frac{\delta g}{\delta V(x)}- \frac{\delta f}{\delta V(x)}\frac{\delta g}{\delta X(x)} \Bigr ), $$ generalizing the standard multivariable expression to an infinite vector of variables, the fields.
If you have several fields, you may index fields to $X^i(x)$ and $V^i(x)$, and you further multiplex the above expression to discrete vectors whose components are infinite-dimensional fields, $$ \left\{ f[X,Y] , g[X,Y]\right \}=\sum_i \int\!\! dx ~~~\Bigl ( \frac{\delta f}{\delta X^i(x)}\frac{\delta g}{\delta V^i(x)}- \frac{\delta f}{\delta V^i(x)}\frac{\delta g}{\delta X^i(x)} \Bigr ). $$
Recall $$ \left\{ X^k(z) , V^j(y)\right \}= \delta^{kj} \delta (z-y)~~. $$