In this paper, we present a rational RBF interpolation method to approximate multivariate functions with poles or other singularities on or near the domain of approximation. The method is based on scattered point layouts and is flexible with respect to the geometry of the problem’s domain. Despite the existing rational RBF-based techniques, the new method allows the use of conditionally positive definite kernels as basis functions. In particular, we use polyharmonic kernels and prove that the rational polyharmonic interpolation is scalable. The scaling property results in a stable algorithm provided that the method be implemented in a localized form. To this aim, we combine the rational polyharmonic interpolation with the partition of unity method. Sufficient number of numerical examples in one, two and three dimensions are given to show the efficiency and the accuracy of the method.