Traditional parametric methods for classification of land use and land cover (LULC) types using remote sensing imagery assume a global distribution model and fail to consider local variation of categorical variables. Differently, non-parametric methods do not make any statistical assumptions but are typically sensitive to the sample sizes of training sample data that usually require a high cost to collect in the field. Geostatistical classifiers, such as indicator kriging and simulation, are local variability-based methods that exhibit great potential for image-based classification of LULC types. However, variogram models required are highly sensitive to the spatial configuration of training samples as well as sample size given a study area. Moreover, when a large number of spectral variables are considered into kriging systems, modeling the variograms and cross-variograms would be problematic. To circumvent these issues, this study extended the geostatistical methods from a 2-dimensional geographic space to a m -dimensional image feature space to derive feature-space indicator variograms (FSIVs). Moreover, a novel stochastic simulation classification algorithm, Feature-Space Indicator Simulation (FSIS), was proposed and examined for classification of LULC types in Duolun County located in Inner Mongolia and in Huang-Feng-Qiao (HFQ) forest farm, Hunan of China. In Duolun, six LULC types were involved and in HFQ a complicated forest landscape consisting of nine forest types plus water, built-up area, and agricultural/bare soil, was classified. The classification results of FSIS were compared with another feature-space geostatistical classifier – feature-space indicator kriging (FSIK), a traditional parametric method – maximum likelihood (ML), a widely used nonparametric method – support vector machine (SVM), and a recently popular method – random forest (RF). The results showed that compared with ML, SVM and RF, in both study areas FSIS statistically significantly increased the accuracy of the classifications by 10.0–29.9% for percentage correct and 19.0–47.6% for Kappa statistic. Compared with FSIK, FSIS also improved the classification accuracy but the accuracy increases were relatively smaller with the percentages correct of 3.5% and 7.6% and the Kappa values of 4.6% and 8.6% for Duolun and HFQ, respectively. Moreover, FSIS led to the spatial uncertainties of the classification estimates as the quality measure of the estimates. In addition, the results also demonstrated that FSIVs were sensitive to the within-class heterogeneity but not very much to the size of training samples. Overall, FSIS exhibited the greater potential to improve the classification accuracy of LULC and forest types using remote sensing image. [ABSTRACT FROM AUTHOR]