Abstract Seismic anisotropy provides crucial information on the stress state and geodynamic processes inside the Earth. We develop a novel adjoint-state traveltime tomography method using P-wave traveltime data to simultaneously determine velocity heterogeneity and azimuthal anisotropy of the subsurface. First, an anisotropic eikonal equation is derived to model first-arrival traveltimes in azimuthally anisotropic media. Traveltime tomography is then formulated as an optimization problem constrained by the anisotropic eikonal equation, which is subsequently solved by the adjoint-state method. Ray tracing is not required. Its high accuracy is achieved by solving the anisotropic eikonal equation and the associated adjoint equation with efficient numerical solvers. In addition, an eikonal equation-based earthquake location method for azimuthally anisotropic media is developed to solve the coupled hypocenter-velocity problem. The tomography and earthquake location methods are applied to central California near Parkfield to test their performance in practice. A total of 1,068,850 first P-wave traveltimes clearly maps the velocity heterogeneity and azimuthal anisotropy in the upper and middle crust. The average P-wave velocity model shows a striking velocity contrast across the San Andreas Fault (SAF). In the upper crust, we find structural anisotropy in the SAF zone and stress-induced anisotropy off the SAF zone. In the middle crust, the fast P-wave velocity directions are generally fault-parallel due to the decreased effect of the maximum horizontal compressive stress. In all, the real-data application suggests that the new adjoint-state traveltime tomography method can be reliably used to investigate anisotropic seismic structures.