Understanding spatial variation in individual species and community assemblages are fundamental goals of both applied and theoretical ecology. Bird communities globally have recently shown massive declines in abundance, leading to increased interest in quantifying variation in individual species distributions as well as biodiversity metrics of avian communities across large spatial domains. This requires modeling spatially dependent, binary detection-nondetection data for large community assemblages (e.g., 300 species) across vast spatial domains. With this as a motivating example, we develop a spatial factor Nearest Neighbor Gaussian Process model for high-dimensional binary outcomes. We leverage Pólya-Gamma latent variables to yield an efficient Markov Chain Monte Carlo sampler and discuss its implementation in the spOccupancy R package. We estimate distributions of 98 bird species across the continental U.S. using our proposed model, in which we embed a detection sub-model to account for imperfect detection of bird species. We subsequently generate maps of multiple avian species richness with associated uncertainty across the continental U.S.