The Epidemic Type Aftershock Sequence (ETAS) model is widely employed to model the spatiotemporal distribution of earthquakes, generally using spatially invariant parameters. We propose an efficient method for the estimation of spatially varying parameters, using the expectation maximization (EM) algorithm and spatial Voronoi tessellation ensembles. We use the Bayesian information criterion (BIC) to rank inverted models given their likelihood and complexity and select the best models to finally compute an ensemble model at any location. Using a synthetic catalog, we also check that the proposed method correctly inverts the known parameters. We apply the proposed method to earthquakes included in the Advanced National Seismic System catalog that occurred within the time period 1981–2015 in a spatial polygon around California. The results indicate significant spatial variation of the ETAS parameters. We find that the efficiency of earthquakes to trigger future ones (quantified by the branching ratio) positively correlates with surface heat flow. In contrast, the rate of earthquakes triggered by far‐field tectonic loading or background seismicity rate shows no such correlation, suggesting the relevance of triggering possibly through fluid‐induced activation. Furthermore, the branching ratio and background seismicity rate are found to be uncorrelated with hypocentral depths, indicating that the seismic coupling remains invariant of hypocentral depths in the study region. Additionally, triggering seems to be mostly dominated by small earthquakes. Consequently, the static stress change studies should not only focus on the Coulomb stress changes caused by specific moderate to large earthquakes but also account for the secondary static stress changes caused by smaller earthquakes.