We consider monotonic, multiple regression for contiguous regions. The regression functions vary regionally and may exhibit spatial structure. We develop Bayesian nonparametric methodology that permits estimation of both continuous and discontinuous functional shapes using marked point process and reversible jump Markov chain Monte Carlo techniques. Spatial dependence is incorporated by a flexible prior distribution which is tuned using cross-validation and Bayesian optimization. We derive the mean and variance of the prior induced by the marked point process approach. Asymptotic results show consistency of the estimated functions. Posterior realizations enable variable selection, the detection of discontinuities and prediction. In simulations and in an application to a Norwegian insurance data set, our methodology shows better performance than existing approaches.