We present a detailed study of the charmonium spectrum using anisotropic lattice QCD. We first derive a tree-level improved clover quark action on the anisotropic lattice for arbitrary quark mass. The heavy quark mass dependences of the improvement coefficients, i.e. the ratio of the hopping parameters $zeta=K_t/K_s$ and the clover coefficients $c_{s,t}$, are examined at the tree level. We then compute the charmonium spectrum in the quenched approximation employing $xi = a_s/a_t = 3$ anisotropic lattices. Simulations are made with the standard anisotropic gauge action and the anisotropic clover quark action at four lattice spacings in the range $a_s$=0.07-0.2 fm. The clover coefficients $c_{s,t}$ are estimated from tree-level tadpole improvement. On the other hand, for the ratio of the hopping parameters $zeta$, we adopt both the tree-level tadpole-improved value and a non-perturbative one. We calculate the spectrum of S- and P-states and their excitations. The results largely depend on the scale input even in the continuum limit, showing a quenching effect. When the lattice spacing is determined from the $1P-1S$ splitting, the deviation from the experimental value is estimated to be $sim$30% for the S-state hyperfine splitting and $sim$20% for the P-state fine structure. Our results are consistent with previous results at $xi = 2$ obtained by Chen when the lattice spacing is determined from the Sommer scale $r_0$. We also address the problem with the hyperfine splitting that different choices of the clover coefficients lead to disagreeing results in the continuum limit.