The Neutron star Interior Composition Explorer (NICER) is currently observing the x-ray pulse profiles emitted by hot spots on the surface of rotating neutron stars allowing for an inference of their radii with unprecedented precision. A critical ingredient in the pulse profile model is an analytical formula for the oblate shape of the star. These formulas require a fitting over a large ensemble of neutron star solutions, which cover a wide set of equations of state, stellar compactnesses and rotational frequencies. However, this procedure introduces a source of systematic error, as (i) the fits do not describe perfectly the surface of all stars in the ensemble and (ii) neutron stars are described by a single equation of state, whose influence on the surface shape is averaged out during the fitting procedure. Here we perform a first study of this systematic error, finding evidence that it is subdominant relative to the statistical error in the radius inference by NICER. We also find evidence that the formula currently used by NICER can be used in the inference of the radii of rapidly rotating stars, outside of the formulas domain of validity. Moreover, we employ an accurate enthalpy-based method to locate the surface of numerical solutions of rapidly rotating neutron stars and a new highly accurate formula to describe their surfaces. These results can be used in applications that require an accurate description of oblate surfaces of rapidly rotating neutron stars.