## 1 Introduction  
Precision agriculture is a new trend in modern agricultural practices as a way to improve yields and quality, thereby increasing the benefits to farmers and guaranteeing sufficient environmental protection [1; 2; 3; 4; 5]. More elaborate soil mapping that reflects the spatial variability of a target soil variable is required to solve various precision agriculture problems, such as soil quality evaluations, soil erosion risk mapping, and variable-rate fertilization [2; 6; 7; 8; 9]. The design of the soil survey is a necessary step in producing farm soil maps, including field sampling, laboratory analysis, and data processing [10; 11; 12; 13]. Field sampling, which is a critical issue in soil surveys, helps to obtain reliable soil mapping results by optimizing the sample size and recognizing the representative sampling locations [14; 15; 16]. Extensive sampling methods have been developed for soil mapping, ranging from simple random sampling to complex, advanced sampling [14; 17; 18; 19; 20; 21; 22; 23; 24; 25; 26; 27; 28; 29; 30; 31; 32]. These sampling methods can be classified into three types: geometric sampling (such as grid sampling [30; 31] and spatial coverage sampling [20]), adapted experimental sampling (such as conditional Latin hypercube sampling [26] and response surface sampling [33]), and model-based sampling (such as kriging model sampling [14]). In the existing sampling methods, the focus is mainly on two aspects that improve farm soil mapping accuracy: sample size and sample locations [14; 34]. It is known that more samples can result in more reliable soil maps, regardless of the sampling methods [35]. However, in most cases, the sample size is determined by available resources, such as field labor, soil physicochemical analysis, and funding [11; 15; 36]. When the sample size is limited by the available resources and cannot meet the requirements of the precision constraints, the determination of the representative sample locations will become critical to obtaining more accurate soil maps. Additionally, the recognition of representative samples is likely to reduce the sample size since it helps to eliminate the non-representative samples. The study by [PERSON] et al. [19] indicated that soil mapping results using representative samples had comparable accuracies to conclusions that were produced using full samples. Sampling design can provide a blueprint for selecting representative sample locations and producing reliable inputs for soil mapping [37]. Most proposed sampling designs in soil mapping create good spatial coverage on one region or good feature space coverage of the covariates for the target soil attribute [11; 16; 26; 30; 38; 39]. Meanwhile, most studies have improved the efficiency of alternative sampling designs using a given mapping method [40; 41; 42; 43; 44; 45]. For example, [PERSON] et al. [45] optimized sample locations using spatial simulated annealing for kriging with an external drift. If we use machine learning methods to produce soil maps, good spatial coverage may not be vital compared to good coverage of the environmental covariates [26]. However, accurate interpolation results using the sample data are very dependent on good geographic dispersion [43]. Nevertheless, in many situations, we may not have determined the mapping method when implementing the sampling designs. Thus, it is essential to develop a sampling design that is robust against deviations in modeling assumptions. A sampling design that considers both good spatial coverage and feature space coverage may be a reliable alternative to solve this issue. Remote sensing technology can be objective and accurate enough to provide sufficient information that is highly correlated with the target farm soil properties, including the normalized difference vegetation index (NDVI), crop yield estimations, the leaf area index (LAI), and the digital elevation model (DEM).