We study the constraint on f(R) gravity that can be obtained by photometric primary probes of the Euclid mission. Our focus is the dependence of the constraint on the theoretical modelling of the nonlinear matter power spectrum. In the Hu–Sawicki f(R) gravity model, we consider four different predictions for the ratio between the power spectrum in f(R) and that in Λ cold dark matter (ΛCDM): a fitting formula, the halo model reaction approach, ReACT, and two emulators based on dark matter only N-body simulations, FORGE and e-Mantis. These predictions are added to the MontePython implementation to predict the angular power spectra for weak lensing (WL), photometric galaxy clustering, and their cross-correlation. By running Markov chain Monte Carlo, we compare constraints on parameters and investigate the bias of the recovered f(R) parameter if the data are created by a different model. For the pessimistic setting of WL, one-dimensional bias for the f(R) parameter, log10| fR0|, is found to be 0.5σ when FORGE is used to create the synthetic data with log10| fR0| = −5.301 and fitted by e-Mantis. The impact of baryonic physics on WL is studied by using a baryonification emulator, BCemu. For the optimistic setting, the f(R) parameter and two main baryonic parameters are well constrained despite the degeneracies among these parameters. However, the difference in the nonlinear dark matter prediction can be compensated for the adjustment of baryonic parameters, and the one-dimensional marginalised constraint on log10| fR0| is biased. This bias can be avoided in the pessimistic setting at the expense of weaker constraints. For the pessimistic setting, using the ΛCDM synthetic data for WL, we obtain the prior-independent upper limit of log10| fR0| < −5.6. Finally, we implement a method to include theoretical errors to avoid the bias due to inaccuracies in the nonlinear matter power spectrum prediction.