Global and regional sea level variations are important indicators of climate change and are derived from accurate sea surface height measurements and precisely determined orbits of altimetry satellites. To validate and improve the quality of these orbits, comparisons with external solutions are important. Since orbit solutions of different institutions are not necessarily provided at the same time instants, interpolation is required for comparison. In this study, we investigate the appropriate interpolation method and its degree to reduce interpolation errors to sub-millimetre levels. We also assess the magnitude of errors occurring at transformations when expressing orbit differences not only in the terrestrial reference frame (Cartesian coordinates), but also in local orbital and ellipsoidal coordinates. The analyses conducted in this study provide good results for Hermite interpolation of degrees 7–11 and Newton interpolation of at least degree 9 with a three-dimensional interpolation error of 0.6 mm and a scattering of 0.2 mm on average for satellite coordinates given with an accuracy of 1 mm in the SP3 format. These interpolation settings limit transformation errors between coordinate systems to ±0.01 mm and incorrect mapping of interpolation errors into certain components in the target system to ±0.02 mm. The spectral analysis of orbit differences is affected up to 0.1 mm in magnitude with appropriate interpolation settings. Extending the number of decimal digits of the satellite position and velocity in SP3 files by one digit benefits the orbit comparisons and reduces the interpolation error by 90% from 0.6 to 0.06 mm. The results are obtained using piece-wise interpolation and a validity interval inside the interpolation interval to minimise the effects of the Runge phenomenon.