When radial basis functions (RBFs) are made increasingly flat, the interpolation error typically decreases steadily until some point when Runge-type oscillations either halt or reverse this trend. Because the most obvious method to calculate an RBF interpolant becomes a numerically unstable algorithm for a stable problem in the case of near-flat basis functions, there will typically also be a separate point at which disastrous ill-conditioning enters. We introduce here a new method, RBF-QR, which entirely eliminates such ill-conditioning, and we apply it in the special case when the data points are distributed over the surface of a sphere. This algorithm works even for thousands of node points, and it allows the RBF shape parameter to be optimized without the limitations imposed by stability concerns. Since interpolation in the flat RBF limit on a sphere is found to coincide with spherical harmonics interpolation, new insights are gained as to why the RBF approach (with nonflat basis functions) often is the more accurate of the two methods.