Spherical regression, in which both covariates and responses lie on the sphere, arises in many scientific applications and has attracted considerable methodological attention in recent years. Despite this progress, constructing flexible and expressive regression models between spherical domains remains challenging, particularly because a single global mapping is often insufficient to capture complex relationships across the entire sphere. A natural strategy is therefore to partition the spherical domain and allow distinct mappings within each region, though this introduces the additional challenge of modeling the partition structure itself. To address these issues, we propose an approach based on optimal transport to model spherical partitions, combined with parametric mappings defined locally within each region. We adopt a Bayesian framework to jointly model both the partitioning and the associated regression maps. This framework enables the identification of heterogeneous regions on the sphere while providing principled uncertainty quantification. Through real-data applications, we demonstrate that the proposed method achieves strong predictive performance, yields meaningful uncertainty estimates, and reveals interpretable clustering structure in spherical data.