The LOFAR radio telescope measures air showers in the energy range $10^{17}$ to $10^{18}$ eV. For each measured shower, the depth of shower maximum Xmax is reconstructed by simulating the radio signal for an ensemble of showers using Corsika and CoREAS. Fitting their radio ‘footprints’ on the ground to the measured radio data yields an Xmax estimate to a precision of about 20 g/cm$^2$. Compared to previous works, we have improved the method in several ways. Local atmospheric data and refractive index profiles are now included into the simulations. The energy estimate and the fitting procedure are now done using the radio signals only, thus limiting systematic uncertainties due to the particle detector array (LORA). Using selection criteria from a more elaborate characterisation of the radio and particle detection, we reduce a composition bias in the Xmax reconstruction. A possible residual bias has been bounded from above. Thus, the systematic uncertainties on <Xmax> have been lowered, reducing an important limiting factor for composition studies at any level of statistics.