Ecological datasets tabulating hundreds or thousands of species are now becoming available. Identifying and understanding the underlying clustering or communities of species which comprise these datasets is difficult with existing joint species distribution models as they usually model dependence on covariates or spatial effects at the species level rather than at the community level. This study describes a novel Bayesian statistical model incorporating multiple model components from existing work designed to represent the communities of species and the dependence of each community on covariates. This model is also extended to produce per-community maps of spatial residuals designed to help identify useful missing covariates. We assess this model via a simulation study and it is applied to analyse a subset of the Florabank1 plant species database with approximately 2,000 sampling units across 900 species. We present a streamlined workflow for interpreting clusters of co-occurring species along with their salient environmental and spatial characteristics.