Small pelagic fish aggregate within areas of suitable habitat to form patchy distributions with localized peaks in abundance. This presents challenges for geostatistical methods designed to investigate the processes underpinning the spatial distribution of stocks and simulate distributions for further analysis. In two-stage models, presence–absence is treated as separable and independent from the process explaining nonzero densities. This is appropriate where gaps in the distribution are attributable to one process and conditional abundance to another, but less so where patchiness is attributable primarily to the strong schooling tendencies of small pelagic fish within suitable habitat. We therefore developed a new modelling framework based on a truncated Gaussian random field (GRF) within a Bayesian framework. We evaluated this method using simulated test data and then applied it to acoustic survey data for Peruvian anchoveta (Engraulis ringens). We assessed the method’s performance in terms of posterior densities of spatial parameters, and the density distribution, spatial pattern, and overall spatial distribution of posterior predictions. We conclude that Bayesian posterior prediction based on a truncated GRF is effective at reproducing the patchiness of the observed spatial distribution of anchoveta.