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., 100 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 Polya-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 United States using our proposed model, in which we embed a detection sub-model to account for imperfect detection of bird species. We subsequently generate species richness maps of two bird communities with associated uncertainty across the continental United States.