Properly accounting for heterogeneity in evolutionary processes is crucial for the development of more realistic phylogenetic inference models. Nucleotide substitution rates and patterns are known to vary among sites in multiple sequence alignments, and progress has been made by partitioning alignments into categories corresponding to different substitution models. Determining partitions can be difficult, and flexible Bayesian nonparametric approaches have been put forth as solutions. However, these previous approaches ignore spatial information, and their use has been limited by the challenges of scaling efficiently to large data sets frequently encountered in infectious disease research. We address these issues through a generalized Polya urn process model framework that subsumes many widely-used countable mixture models, including infinite hidden Markov models that account for spatial information. To improve computational performance, we adapt data squashing techniques for mixture models to develop efficient Markov chain Monte Carlo (MCMC) transition kernels, and we exploit algorithms to parallelize computations for likelihood calculation across processing cores.