(Figure presented.) is a traffic-related air pollutant. Ground (Figure presented.) monitoring stations measure (Figure presented.) concentrations at certain locations and statistical predictive methods have been developed to predict (Figure presented.) as a continuous surface. Among them, ensemble tree-based methods have shown to be powerful in capturing nonlinear relationships between (Figure presented.) measurements and geospatial predictors but it is unclear if the spatial structure of (Figure presented.) is also captured in the response-covariates relationships. We dive into the comparison between spatial and nonspatial data models accounting for prediction accuracy, model interpretation and uncertainty quantification. Moreover, we implement two new spatial and a nonspatial methods that have not been applied to air pollution mapping. We implemented our study using national ground station measurements of (Figure presented.) in Germany and the Netherlands of 2017. Our results indicate heterogeneous levels of importance of modeling the spatial process in different areas. The prediction intervals predicted with ensemble tree-based methods are more satisfactory than the geostatistical methods. The two new methods implemented each obtained better prediction accuracy compared to the original ensemble tree-based and stacking methods. The probabilistic distribution of the spatial random field estimated by the geostatistical methods could provide useful information for analyzing emission sources and the spatial process of observations.