We provide a technique to efficiently produce high-resolution three-dimensional surface motion maps by combining information about the motion of the Earth's surface from interferometric observations of synthetic aperture radar images and repeated Global Positioning System (GPS) geodetic measurements. Unwrapped interferograms, showing pixel-wise change in range from ground to satellite, and sparse values of three-dimensional movements are required as input. The problem of finding the full three-dimensional motion field is separated into two two-dimensional problems. Initially, the vertical component of the deformation field and its horizontal component in the look direction of the satellite are found. Later, the look direction component is resolved into north and east components. Initial values for the motion fields are assigned to each pixel of interferograms from interpolation of available GPS observations. These values are then updated and optimized by comparison with the interferograms and the GPS observations. An additional constraint is an assumption of a smoothly varying motion field. Markov random field-based regularization and simulated annealing algorithm are used for the optimization. The technique is applied to create surface motion maps for the Reykjanes Peninsula, SW Iceland.