This article proposes a geostatistical interpolation method, namely, block indicator kriging, as a means to model uncertainty-quantified presettlement vegetation surfaces. We demonstrate its potential in landscape ecology for improving the quality of the resulting surfaces using indicator-coded
presettlement land survey record data and the areal proportion of each species. The geostatistical interpolation method presented in this study explicitly models support differences between the source data and the prediction surface, while taking into account spatial dependence present in
multiscale presettlement land survey record data. In this case study, we demonstrate that block indicator kriging with areal proportion data substantially increases the prediction accuracy and coherence using witness tree species data recorded in the northeast of Minnesota. The relative merit
of the proposed geostatistical interpolation method to other conventional geostatistical approaches is illustrated through a comparative analysis, where continuous vegetation surfaces obtained from other approaches are compared with those obtained from the block indicator kriging with areal
proportion data in terms of the prediction accuracy and consistency.