This paper introduces a novel data-driven methodology for estimating the full random field (random function) of a regionalized variable, building upon and extending the concepts of ensemble spatial interpolation and adaptive ensemble spatial analysis. The core contribution is the development of a rigorous theoretical framework that proves how the method implicitly learns the underlying spatial dependence structure. By inferring local predictive distributions through adaptively defined, overlapping spatial partitions, the approach ensures that spatial coherence is governed by a learned statistical copula. Formal proof is provided demonstrating that this dependency structure converges to the actual structure of the regionalized variable. Therefore, the resulting spatial distributional estimator is consistent, captures all the uncertainty of the random field, and naturally allows for the generation of multiple equally probable realizations. This non-parametric strategy offers a flexible alternative, inherently preserving spatial patterns and capturing fine-scale variability without requiring prior model specification. Experiments conducted in a detailed case study applied to geostatistical simulation, using both synthetic and real datasets, confirm the effectiveness and computational efficiency of the method, demonstrating its ability to recover local statistics and spatial structure with greater robustness compared to conventional techniques.