In regression analysis, environmental data often exhibit spatial heterogeneity in regression coefficients. Spatially varying coefficient models, including geographically weighted regression and bivariate spline models, are standard tools to quantify such nonstationarity. In this paper, we propose a spatially varying coefficient model by representing the spatially-varying parameters as a mixture of local polynomials at selected locations. The local polynomial parameters have attractive interpretations, indicating various types of spatial nonstationarity. Instead of estimating the spatially-varying regression coefficients directly, we develop a penalized least squares regression procedure for the local polynomial parameter estimation, which penalizes the differences in neighboring parameters and selects the number of mixture components simultaneously. The proposed model is applied to describe the association between particulate matter concentrations (PM 2.5) and other pollutants related to the secondary aerosol formulation, such as sulfate and nitrate. The identified spatial subregions show distinct relationships between the PM 2.5 and other pollutants we have considered.