A double grading of catalan numbers
It appears that the zeta map which I used to answer Vince Vatter's initial question and which I describe in my answer there, see also page 50 of Jim Haglund's book, indeed solves also this problem:
As discussed in a comment above, the zeta map has the property that it sends
- the number of $0$'s in an "area sequence" $a = (a_1,\ldots,a_n)$ (item $u$ in Stanley's list, see also Property $D$ in my answer to the other question) to the number of initial north (or up) steps, and
- the number of $i$'s for which all $i$'s appear before all $i+1$'s within $a$ to the number of returns to $0$ (excluding the very last return).
On the other hand, David's bijection between rooted planar trees and Dyck paths sends
- the number of children of the root to the number of $0$'s in the corresponding area sequence, and
- the number of crucial vertices in a tree (excluding the root which is always, except for $n=1$, crucial) to the number of $i$'s for which all $i$'s appear before all $i+1$'s within the corresponding area sequence.
Combining his bijection with the zeta map yields therefore a bijection between rooted planar trees and Dyck paths that sends the bistatistic
(number of crucial vertices, number of children of the root)
to the bistatistic
(number of returns to $0$, number of initial north steps).
To extend David's initial example of the $5$ rooted planar trees of size $4$, the $5$ area sequences corresponding to them (in the ordering $$(1,2), (2,1), (3,1), (2,2), (1,3)$$ as above) are $$010, 011, 012, 001, 000.$$ Applying the zeta map yields the area sequences $$011, 001, 000, 010, 012.$$ For those, the (sequence of the) number of returns to $0$ is $1,2,3,2,1$ and the (sequence of the) number of initial north steps is $2,1,1,2,3$. We thus obtain the same bistatistical sequence.
We thus proved that $\sum c(p,q,n)x^py^q$ is indeed equal to the Tutte polynomial $T_n(x,y)$.
I wanted to add a simple combinatorial proof that the coefficient of $x^p y^q$ is dependent only on $p+q$. Given a Dyck path of length $2n$, let $q$ be the number of returns to $0$ and let the initial segment be $0 1 2 3 \cdots (p-1) p (p-1)$. I will introduce an operation which sends a path of type $(p,q)$, with $q \geq 2$, to one of type $(p+1, q-1)$.
So, take a Dyck path with $q \geq 2$. Write it as a matching parentheses sequence; we can decompose it as $$( D_1 ) ( D_2 ) ( D_3 ) \cdots ( D_q )$$ where each $D_i$ is, itself, a matching parentheses sequence. We map it to $$((D_1)D_2) (D_3) \cdots (D_q)$$ It is easy to see that this map is invertible for $p \geq 2$.