Gaussian processes are widely used for the analysis of spatial data due totheir nonparametric flexibility and ability to quantify uncertainty, andrecently developed scalable approximations have facilitated application tomassive datasets. For multivariate outcomes, linear models of coregionalizationcombine dimension reduction with spatial correlation. However, theirreal-valued latent factors and loadings are difficult to interpret because,unlike nonnegative models, they do not recover a parts-based representation. Wepresent nonnegative spatial factorization (NSF), a spatially-awareprobabilistic dimension reduction model that naturally encourages sparsity. Wecompare NSF to real-valued spatial factorizations such as MEFISTO andnonspatial dimension reduction methods using simulations and high-dimensionalspatial transcriptomics data. NSF identifies generalizable spatial patterns ofgene expression. Since not all patterns of gene expression are spatial, we alsopropose a hybrid extension of NSF that combines spatial and nonspatialcomponents, enabling quantification of spatial importance for both observationsand features. A TensorFlow implementation of NSF is available fromhttps://github.com/willtownes/nsf-paper .