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 by matching the Hamiltonian on the lattice and in the continuum. The heavy quark mass dependence of the improvement coefficients, i.e., the ratio of the hopping parameters zeta=K-t/K-s and the clover coefficients c(s,t), is examined at the tree level, and effects of the choice of the spatial Wilson parameter r(s) are discussed. 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 with r(s)=1 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. The latter employs the condition that the speed of light calculated from the meson energy-momentum relation be unity. We calculate the spectrum of S and P states and their excitations using both the pole and kinetic masses. We find that the combination of the pole mass and the tadpole-improved value of zeta to yield the smoothest approach to the continuum limit, which we then adopt for the continuum extrapolation of the spectrum. 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 similar to30% for the S-state hyperfine splitting and similar to20% 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. Making a leading order analysis based on potential models we show that a large hyperfine splitting similar to95 MeV obtained by Klassen with a different choice of the clover coefficients is likely an overestimate.