Quantifying distributions of wildlife and plant species across space and time is a critical task for understanding global change effects and developing effective conservation and management actions. However, many modern ecological data sets present a variety of statistical and computational complexities when modeling ecological processes for hundreds of species across thousands of spatial locations. In this talk, I will discuss a flexible hierarchical Bayesian modeling framework designed to model geostatistical data in which both the number of locations and number of outcomes at each location is large. The proposed framework addresses these source of high dimensionality via spatial latent factors and highly scalable Nearest Neighbor Gaussian Process (NNGP) models. The approach accounts for autocorrelation within each response across spatial locations as well as among responses at any individual location. I will present two motivating case studies in which I use the proposed framework to predict forest biomass at the species-level across the US and quantify nonstationary effects of climate and land use change on a forest bird community in the eastern US.