Large, spatially correlated data often covers a vast geographical area. Consequently, it may not be reasonable to assume that the spatial process is stationary. Local spatial regions may have different mean and covariance structures. We propose a model selection methodology that accomplishes three goals: (1) cluster locations into small regions with distinct, stationary models, (2) perform Bayesian model selection within each cluster, and (3) correlate the model selection in nearby clusters. We utilize the Ising distribution to provide correlation for model selection indicators across each local Gaussian process model. We will present an application of the methodology on brook trout abundance across Pennsylvania.