Measurements of 189 continuous and 933 campaign-mode Global Positioning System (GPS) stations with 3–16 yr data spans over the Tibetan Plateau reveal contemporary three-dimensional (3-D) crustal deformation during 1999–2016. The Empirical Orthogonal Function method was used to characterize the spatial variations in the surface deformation with distinct seasonal oscillations at the GPS sites in five regions of the Tibetan Plateau. We find that these surface variations are highly correlated with the corresponding mass load signals observed by the Gravity Recovery and Climate Experiment (GRACE) mission. The improved GPS processing strategy used to determine the 3-D velocity field includes maximum likelihood estimation, removal of common mode errors from GPS time series using Principal Component Analysis (PCA), and power law plus white noise stochastic error modeling. We determined the rates of vertical crustal movement by removing GRACE-observed non-tectonic origin load deformation, 2002–2016. The corrected vertical crustal deformation shows that the Himalaya region is uplifting at an average rate of ∼1.7 mm yr−1, and that the northeastern Tibetan Plateau is uplifting at an average rate of ∼1.3 mm yr−1. In addition, the horizontal velocity relative to the stable Eurasian plate and its corresponding dilatation throughout the Tibetan Plateau suggest that tectonic shortening and crustal thickening is occurring at −90 to −80 nanostrain yr−1 in the southern Tibetan Plateau and −30 to −20 nanostrain yr−1 in the northeastern Tibetan Plateau, which could be related to the geologic shortening and elastic strain accumulation. The interior Tibetan Plateau exhibits crustal thinning and block movement along strike-slip faults. Eastward motion of the crust north of the Xianshuihe-Xiaojiang Fault system relative to crust to its south results in shear strain and reflects eastward escape of plastic crustal material in the southeastern Tibetan Plateau.