The objective of this study was to incorporate geostatistics, remote sensing and geographic information system methods due to improving the quantitative land suitability assessment in Arsanjan plain, southern Iran. The primary data was collected from 85 soil samples from tree depths (030, 3060 and 6090 cm) and the secondary information from remotely sensed data “LISSIII receiver from IRSP6 satellite”. In order to identify the spatial dependence of soil important parameters, we used ordinary kriging and simple kriging with varying local means (SKVLM) methods. The results indicated that best method with the lowest mean square error for mapping pH and electrical conductivity (ECe) (030 cm) obtained from SKVLM method that spectral values of band 1 of LISSIII receiver was used as secondary variable. While, other soil properties indicated moderate to strong spatial dependence in the study area and interpolated in unstamped point by ordinary kriging method with the reliable accuracy. The land suitability evaluation method (parametric) has applied on the density points (150 ¥ 150 m2) that obtained by kriging or SKVLM methods, instead of applying on the limited representative profiles conventionally. Overlaying the information layers of dada was used by GIS for preparing the final land suitability evaluation. Therefore, changes in land characteristics as locally could be identified even in the same soil units. In addition, it is recognized that many of the land characteristics vary over a short distance within soil uniform mapping units. In general, this new method can easily present the squares and limitation factors of different land suitability classes with considerable accuracy in arbitrary land indices.