We discuss the flow harmonics or the elliptic and triangular flow of $J/\psi$, $\psi(2S)$, and $\chi_c(1P)$ mesons in heavy ion collisions. Starting from the investigation on transverse momentum distributions of those charmonium states, we calculate their elliptic and triangular flow when they are produced at the quark-hadron phase boundary by quark recombination. We show that the wave function distribution of charmonium states plays a significant role, especially in producing charmonium states, leading to the transverse momentum distribution of the $\psi(2S)$ meson as large as that of the $J/\psi$ meson. On the other hand, we find that the wave function effects and feed-down contributions are averaged out for elliptic and triangular flow, resulting in similar elliptic and triangular flow for all charmonium states. We further investigate the elliptic and triangular flow of charmonium states at low transverse momentum regions and discuss the quark number scaling of elliptic and triangular flow for charmonium states in heavy ion collisions.