We developed a three-dimensional (3D) shear-wave splitting tomography method to image the spatial anisotropy distribution by back projecting shear wave splitting delay times along ray paths derived from a 3D shear velocity model, assuming the delay times are accumulated along the ray paths. The local strength of the anisotropy is indicated by a parameter of anisotropy percentage, K. Using the shear-wave splitting delay times for 575 earthquakes measured at PASO and HRSN stations, we imaged a detailed 3D anisotropy percentage model around the San Andreas Fault Observatory at Depth (SAFOD). The anisotropy percentage model shows strong heterogeneities, consistent with the strong spatial variations in both measured delay times and fast polarization directions. The San Andreas Fault (SAF) zone is highly anisotropic down to a depth of ∼4 km and then becomes less anisotropic at greater depths. Outside the fault zone, the highly anisotropic zone extends as deep as ∼7 km, consistent with the systematic depth dependence of the average time delays. To the southwest of the SAF, the Salinian granitic block shows relatively strong anisotropic anomalies that are presumably caused by aligned microcracks consistent with the direction of the regional maximum compressive horizontal stress. To the northeast of the fault zone, a strong anisotropic anomaly between depths ∼2 and ∼4 km corresponds to a serpentinite body sandwiched between Franciscan rocks.