Detecting and measuring a non-Gaussian signature of primordial origin in the density field is a major science goal of nextgeneration galaxy surveys. The signal will permit us to determine primordial physics processes and constrain models of cosmic inflation. While traditional approaches utilise a limited set of statistical summaries of the galaxy distribution to constrain primordial non-Gaussianity, we present a field-level approach by Bayesian forward-modelling the entire three-dimensional galaxy survey. Our method naturally and fully self-consistently exploits the entirety of the large-scale structure, e.g., higherorder statistics, peculiar velocity fields, and scale-dependent galaxy bias, to extract information on the local non-Gaussianity parameter, 𝑓 nl . We demonstrate the performance of our approach through various tests with mock galaxy data emulating relevant features of the SDSS-III/BOSS-like survey, and additional tests with a Stage IV mock data set. These tests reveal that the method infers unbiased values of 𝑓 nl by accurately handling survey geometries, noise, and unknown galaxy biases. We demonstrate that our method can achieve constraints of 𝜎 𝑓 nl ≈ 8.78 for SDSS-III/BOSS-like data, an improvement of a factor ∼ 2.5 over currently published constraints. Tests with next-generation mock data show that significant further improvements are feasible with sufficiently high resolution. Furthermore, the results demonstrate that our method can consistently marginalise all nuisance parameters of the data model. The method further provides an inference of the three-dimensional primordial density field, providing opportunities to explore additional signatures of primordial physics.