A method of fast gridding of irregular data, has been developed for spatial interpolation of large irregular spatial point data sets; for example building a 3D geographic terrain grid surface from billions of irregularly spaced xyz coordinates on the earth's surface. The method developed typically translates into many orders of magnitude gain in computational speed. For example, to produce a gridded data set (having M rows and N columns) from P irregularly located sampling points, the computational steps required can be reduced from a number of the order of O(M×N×P) to a lesser number of the order of O(M×N+P) operations. The method achieves this by ensuring that each of the P sampling points is visited only once. This is particularly significant since spatial data collection devices typically collect data points in the billions. The method described is readily extendible to any number of dimensions.