Delineating spatial variations of seismic anisotropy in the crust is of great importance for the understanding of structural heterogeneities, regional stress regime and ongoing crustal dynamics. In this study, we present a 3-D anisotropic P-wave velocity model of the crust beneath northern California by using the eikonal equation-based seismic azimuthal anisotropy tomography method. The velocity heterogeneities under different geological units are well resolved. The thickness of the low-velocity sediment at the Great Valley Sequence is estimated to be about 10 km. The high-velocity anomaly underlying Great Valley probably indicates the existence of ophiolite bodies. Strong velocity contrasts are revealed across the Hayward Fault (2–9 km) and San Andreas Fault (2–12 km). In the upper crust (2–9 km), the fast velocity directions (FVDs) are generally fault-parallel in the northern Coast Range, which may be caused by geological structure; while the FVDs are mainly NE–SW in Great Valley and the northern Sierra Nevada possibly due to the regional maximum horizontal compressive stress. In contrast, seismic anisotropy in the mid-lower crust (12–22 km) may be attributed to the alignment of mica schists. The anisotropy contrast across the San Andreas Fault may imply different mechanisms of crustal deformation on the two sides of the fault. Both the strong velocity contrasts and the high angle (∼45° or above) between the FVDs and the strikes of faults suggest that the faults are mechanically weak in the San Francisco bay area (2–6 km). This study suggests that the eikonal equation-based seismic azimuthal anisotropy tomography is a valuable tool to investigate crustal heterogeneities and tectonic deformation.