This paper presents a fast algorithm for smooth digital elevation model interpolation and approximation from scattered elevation data. The global surface is reconstructed by subdividing it into overlapping local subdomains using a perfectly balanced binary tree. In each tree leaf, a smooth local surface is reconstructed using radial basis functions. Finally a hierarchical blending is done to create the final C1-continuous surface using a family of functions called Partition of Unity. We present two terrain data sets and show that our method is robust since the number of data points in the Partition of Unity blending areas is explicitly specified. Categories and Subject Descriptors