Summary: In this paper we present a fast algorithm for two-dimensional interpolation of large scattered data sets. It is based on the partition of unity method, which makes use of compactly supported radial basis functions as local approximants. This interpolation technique is characterized by the use of an efficiently implemented nearest neighbour searching procedure, which exploits a suitable and optimal partition of the domain in strips obtained in a completely automatic way. Analysis of computational complexity and numerical results show efficiency and accuracy of the proposed algorithm, also considering an application to Earth’s topography.