We have developed a complete theory for the calculation of the observed Stokes parameters for radiation emitted from the surface of a rapidly rotating neutron star (NS) using the oblate Schwarzschild approximation. We accounted for the rotation of the polarization plane due to relativistic effects along the path from the stellar surface to the observer. The results were shown to agree with those obtained by performing full numerical general relativistic ray-tracing with the ARCMANCER code. We showed that the obtained polarization angle profiles may differ substantially from those derived for a spherical star. We demonstrated that assuming incorrect shape for the star can lead to biased constraints for NS parameters when fitting the polarization data. Using a simplified model, we also made a rough estimate of how accurately the geometrical parameters of an accreting NS can be determined using the X-ray polarization measurements of upcoming polarimeters like the Imaging X-ray Polarimeter Explorer or the enhanced X-ray Timing and Polarimetry mission.