A search for the Standard Model Higgs boson produced via vector boson fusion in association with a photon and decaying to a bottom quark-antiquark pair has been conducted using \unit[12.6]{fb$^{-1}$} of proton-proton collision data collected by the ATLAS detector at a center-of-mass energy $\sqrt{s}=\unit[13]{TeV}$ at the LHC in 2015 and 2016. The search benefits from a large reduction of the non-resonant multijet background relative to similar searches that do not explicitly require a photon and from efficient triggering aided by the presence of the photon. Multivariate techniques trained using simulated Monte Carlo samples are used to enhance the sensitivity of the analysis by constructing regions of phase space with higher expected signal fractions relative to the background. Data-driven techniques are used in those regions to provide a reliable estimate of the di-$b$-jet invariant mass spectrum of the non-resonant background. The final statistical analysis of the data relies on a profile likelihood fit to the di-$b$-jet invariant mass distribution, searching for a signal bump in an otherwise smoothly falling distribution. There is no observed excess above the background-only expectation, and the observed (expected) 95\% confidence level upper limit on the production cross section times branching ratio for a Higgs boson mass of \unit[125]{GeV} is 4.0 (6.0$^{+2.3}_{-1.7}$) times the Standard Model expectation. The observed signal strength is $-3.9^{+2.8}_{-2.7}$ times the Standard Model value. This analysis is also sensitive to $\Zboson+\gamma$ production: the observed (expected) upper limit is 2.0 (1.8$^{+0.7}_{-0.5}$) times the Standard Model expectation with an observed signal strength of $0.3\pm0.8$.