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 z5Kt /Ks 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 j5as /at53 anisotropic lattices. Simulations are made with the standard anisotropic gauge action and the anisotropic clover quark action with r s51 at four lattice spacings in the range as50.07‐0.2 fm. The clover coefficientsc s,t are estimated from tree-level tadpole improvement. On the other hand, for the ratio of the hopping parameters z, 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 z 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 1 P-1S splitting, the deviation from the experimental value is estimated to be ;30% for the S-state hyperfine splitting and ;20% for the P-state fine structure. Our results are consistent with previous results at j52 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 ;95 MeV obtained by Klassen with a different choice of the clover coefficients is likely an overestimate.