We present a Bayesian method to identify multiple (chemodynamic) stellar populations in dwarf spheroidal galaxies (dSphs) using velocity, metallicity, and positional stellar data without the assumption of spherical symmetry. We apply this method to a new Keck/Deep Imaging Multi-Object Spectrograph (DEIMOS) spectroscopic survey of the Ursa Minor (UMi) dSph. We identify 892 likely members, making this the largest UMi sample with line-of-sight velocity and metallicity measurements. Our Bayesian method detects two distinct chemodynamic populations with high significance (in logarithmic Bayes factor, lnB~33). The metal-rich ([Fe/H]=-2.05+/-0.03) population is kinematically colder (radial velocity dispersion of {sigma}_v_=4.9_-1.0_^+0.8^km/s) and more centrally concentrated than the metal-poor ([Fe/H]=-2.29_-0.06_^+0.05^) and kinematically hotter population ({sigma}_v_=11.5_-0.8_^+0.9^km/s). Furthermore, we apply the same analysis to an independent Multiple Mirror Telescope (MMT)/Hectochelle data set and confirm the existence of two chemodynamic populations in UMi. In both data sets, the metal-rich population is significantly flattened ({epsilon}=0.75+/-0.03) and the metal-poor population is closer to spherical ({epsilon}=0.33_-0.09_^+0.12^). Despite the presence of two populations, we are able to robustly estimate the slope of the dynamical mass profile. We found hints for prolate rotation of order ~2km/s in the MMT data set, but further observations are required to verify this. The flattened metal-rich population invalidates assumptions built into simple dynamical mass estimators, so we computed new astrophysical dark matter annihilation (J) and decay profiles based on the rounder, hotter metal-poor population and inferred log_10_(J(0.5{deg})/GeV^2^cm^-5^)~19.1 for the Keck data set. Our results paint a more complex picture of the evolution of UMi than previously discussed.