The liquefaction potential index (LPI) is widely used for evaluating the severity of liquefaction manifestation at the ground surface (e.g., settlement, lateral spreading, sand boils, and crack) and for developing liquefaction severity maps. Over the last two decades, several methods, such as cumulative probability distribution of LPI and geostatistics-based LPI mapping, have been proposed to develop a liquefaction severity map from in situ tests (e.g., cone penetration tests, CPT), which are often sparsely performed in sites. These methods are either based on an assumption of statistical homogeneity within each geologic unit or a stationary Gaussian model. However, subsurface soils frequently show significant spatial variability and LPI data obtained at different geological units usually exhibit non-stationary characteristics. More importantly, existing methods offer little insight into the reliability level of the obtained liquefaction severity map. To address these issues, this study proposes a non-parametric and data-drivenmethod for CPT-based liquefaction severity mapping and, for the first time ever, quantification of the liquefaction severity maps’ reliability level using the probability of mis-predicting liquefaction severity from the map. Both synthetic and real-life data are used to demonstrate and validate the proposedmethod. The illustration examples indicate that the proposed method can properly deal with non-stationary LPI data from different geological units and quantify the misprediction probability of liquefaction severity at each point of the map.