We develop an algorithm to integrate GPS and InSAR data for a 3‐dimensional crustal deformation field at the Earth's surface. In the algorithm discrete GPS data points are interpolated to obtain a 3‐dimensional continuous velocity field, which is then combined with the InSAR line‐of‐sight (LOS) velocity data pixel by pixel using the least‐squares method. Advantages of our method over previous ones are that: 1) The GPS data points are optimally interpolated by balancing a trade‐off between spatial resolution and solution stability. 2) A new algorithm is developed to estimate realistic uncertainties for the interpolated GPS velocities, to be used as weights for GPS data in GPS‐InSAR combination. 3) Realistic uncertainties for the InSAR LOS rate data are estimated and used as weights for InSAR data in GPS‐InSAR combination. 4) The ramps and/or offsets of the InSAR data are globally estimated for all the images to minimize data misfit, particularly at regions where the data overlaps. Application of this method to real data from southern California shows its capability of successfully restoring 3‐dimensional continuous deformation field from spatially limited GPS and dimensionally limited InSAR data. The deformation field reveals water withdrawal induced subsidence and drought caused uplift at various regions in southern California.