This paper presents a novel application of the geostatistical multivariate method known as min–max autocorrelation factors (MAFs) for analysing fisheries survey data in a space–time context. The method was used to map essential fish habitats and evaluate the variability in time of their occupancy. Research surveys at sea on marine fish stocks have been undertaken for several decades now. The data are time series of yearly maps of fish density, making it possible to analyse the space–time variability in fish spatial distributions. Space–time models are key to addressing conservation issues requiring the characterization of variability in habitat maps over time. Here, the variability in fisheries survey data series is decomposed in space and time to address these issues, using MAFs. MAFs were originally developed for noise removal in hyperspectral multivariate data and are obtained using a specific double principal components analysis. Here, MAFs were used to extract the most continuous spatial components that are consistent in time, together with the time series of their amplitudes. MAFs formed an empirical isofactorial model of the data, which served for kriging in each year using all available information across the data series. The approach was applied on the spawning distributions of sardine in the Bay of Biscay from 2000 to 2017. A multivariate approach for dealing with space–time data was adapted here, because the evolution in time was highly variable. Maps were classified using the amplitudes of the MAFs, and groups of typical distributions were identified, which showed different occurrence probabilities in different periods.