Agricultural Land Suitability Assessment Using Satellite Remote Sensing-Derived Soil-Vegetation Indices

Satellite remote sensing technologies have a high potential in applications for evaluating land conditions and can facilitate optimized planning for agricultural sectors. However, misinformed land selection decisions limit crop yields and increase production-related costs to farmers. Therefore, the purpose of this research was to develop a land suitability assessment system using satellite remote sensing-derived soil-vegetation indicators. A multicriteria decision analysis was conducted by integrating weighted linear combinations and fuzzy multicriteria analyses in a GIS platform for suitability assessment using the following eight criteria: elevation, slope, and LST vegetation indices (SAVI, ARVI, SARVI, MSAVI, and OSAVI). The relative priorities of the indicators were identified using a fuzzy expert system. Furthermore, the results of the land suitability assessment were evaluated by ground truthed yield data. In addition, a yield estimation method was developed using indices representing influential factors. The analysis utilizing equal weights showed that 43% of the land (1832 km2) was highly suitable, 41% of the land (1747 km2) was moderately suitable, and 10% of the land (426 km2) was marginally suitable for improved yield productions. Alternatively, expert knowledge was also considered, along with references, when using the fuzzy membership function; as a result, 48% of the land (2045 km2) was identified as being highly suitable; 39% of the land (2045 km2) was identified as being moderately suitable, and 7% of the land (298 km2) was identified as being marginally suitable. Additionally, 6% (256 km2) of the land was described as not suitable by both methods. Moreover, the yield estimation using SAVI (R2 = 77.3%), ARVI (R2 = 68.9%), SARVI (R2 = 71.1%), MSAVI (R2 = 74.5%) and OSAVI (R2 = 81.2%) showed a good predictive ability. Furthermore, the combined model using these five indices reported the highest accuracy (R2 = 0.839); this model was then applied to develop yield prediction maps for the corresponding years (2017–2020). This research suggests that satellite remote sensing methods in GIS platforms are an effective and convenient way for agricultural land-use planners and land policy makers to select suitable cultivable land areas with potential for increased agricultural production.


Introduction
Proper land-use planning is essential for enhancing agricultural production and ecological conservation and for the protection of biodiversity [1]. Inappropriate land management practices lead to a higher rate of soil erosion, a diminished crop production, a hindered productivity, and a deteriorated soil quality [2]. Therefore, land management focusing on suitability should be a key issue of research and policy development mainly due to its influence on crop production. The knowledge of local land conditions has become increasingly recognized for its importance in sustainable land management [3]. Farmers of local communities assess their farmland using consistent observations and collective experiences [4]. However, for rural communities, this knowledge is usually insufficient to understand the adequacy of suitable condition, management strategies, and land-use decisions.
In addition, many conventional techniques for earth monitoring applications require specific spectral features that are defined only for multispectral data such as deep learning, exploiting both temporal and cross-sensor dependencies and deep networks achieve much better performance than traditional methods [5,6]. Furthermore, innovative farming technologies incorporate biology with computers and device exchange-based smart agriculture become achieved in a structured farm management system. The high spatial imagery from remote sensing datasets may provide an aid to systematically consider issues associated with smart farming technology. Remote sensing methods support the formation of growth profiles of plants and temporal evolution scheme of soils over their developmental phases [7]. Remote sensing indices that incorporate environmental recovery factors are useful for tracing the development of crops, their interrelatedness, and the consequences of the variables of interest for crop development. Following this concern, the application of smart agriculture and satellite remote sensing-based soil-vegetation index evaluations for agricultural land condition assessments is the key target of this research. Therefore, land suitability assessments can be performed using the multicriteria decision method. Such evaluation provides information about specific land use potentials and constraints. The Multicriteria Decision Method (MCDM) becomes more suitable when incorporating geospatial references. In recent years, computing technologies combined with GIS have enabled geospatial references using MCDM-land suitability evaluations. Furthermore, the MCDM, combined with linear combination and fuzzy set theory, has the potential to reduce subjectivity in the assessment of results. Several approaches to the MCDM that utilize equal-weighted linear combinations or fuzzy membership systems have been applied to conduct land suitability evaluations [8][9][10][11]. In addition, for sustainable land resource management, the Food and Agricultural Organization (FAO) has proposed guidelines for land evaluation [12]. According to the guidelines, land is classified into four categories: highly suitable, moderately suitable, marginally suitable, and not suitable. Additionally, the equal-weighted linear combination-fuzzy overlay technique in the GIS platform has the capabilities needed to overcome these limitations by applying the required calorie ratio (FAO recommended) to prepare land suitability assessments for smart agricultural practices. However, there is a lack of datasets in some areas of developing regions where assessments of land suitability are really challenging. In addition, recent datasets of those geographic information system have limitation, especially land uses, drainage and lack of soil sampling information for soil nutrients in a distant time period.
It is worth mentioning that quick and accurate land suitability assessments can aid in the improvement of yield prediction models. Regarding the judgment of the prediction of yield using vegetation indices (VIs), it is the most straightforward approach to establishing empirical relationships between ground-based yield measures and Vis [13][14][15]. In this regard, satellite remote sensing technologies and GIS applications for monitoring crops have the potential to establish timely assessments of changes in the growth and development of crops on regional scales [16]. The yield prediction is also helpful for making decisions on regional food security policies and production inventories to understand the availability of field crop.
Additionally, local farmer's perceptions and their assessment of land suitability can differ from scientific approaches due to the much broader contextual implications of the former and how they are often framed. This often results in differences in perceived problems and the require solutions [17]. In most cases, developing location specific descriptions by soil sampling and analysis is expensive and challenging. Following this concern, advanced and affordable smart satellite remote sensing multicriteria technologies that consider climate factors are required for land suitability and accuracy assessments.
Therefore, the purpose of this research is to develop a soil-vegetation intent land suitability assessment model based on multicriteria decision-making analysis to determine optimal land distributions according to soil-vegetation indices to ensure elevated productivity, and a yield prediction method was developed using vegetation indices. A methodology that can be applied across various countries is proposed to reduce the complexities of farmland evaluation practices.

Materials and Methods
The proposed method utilizes a GIS-based multicriteria analysis to develop a soilvegetation index-associated suitability analysis by exploiting satellite remote sensing for land suitability assessment and consists of three major steps ( Figure 1): the calculation of soil vegetation indices for land suitability mapping of diversified crops, regression analysis using ground truth yield data for validation, and the utilization of a yield prediction model to develop a yield map. ArcGIS 10.4 ® (ESRI, CA, USA) software was used for criteria aggregation, data preprocessing and calculation standardization, weight determination by an equal-weighted overlay, fuzzy membership function, fuzzy overlay, and raster calculation.

Study Area
The study area is located between 25 • 14 and 26 • 02 N latitudes and 88 • 22 and 89 • 54 E longitudes in the northern part of Bangladesh and has a total area of 8260 square kilometers. The study was conducted in the Dinajpur, Rangpur, Kurigram, and Gaibandha districts of the Rangpur Division where the inhabitants derive their livelihoods from agriculture ( Figure 2). The area consists of 36 administrative units with an overall population of 11,498,000 [18]. The population is economically active in agriculture; nevertheless, agronomic land use is highly inconsistent due to climatic factors, soil property issues, water infiltration, environmental resources, and local socioeconomic conditions. Based on weather data, the minimum and maximum mean annual temperature varies between 8.47 • C and 36.3 • C. The annual average rainfall recorded is 7650-1233 mm, with a high humidity in the range of 41-77% [19]. The elevation ranges from 5 to 30 m above sea level.

Image Acquisition
Landsat 8 is the most recently launched Landsat satellite and carries the Operational Land Imager (OLI) and Thermal Infrared Sensor (TIRS) instruments. OLI collects data in the visible, near infrared, and shortwave infrared spectral bands and in a panchromatic band. These two sensors provide seasonal coverage of the global landmass at a spatial resolution of 30 m (visible, NIR, and SWIR), 100 m (thermal), and 15 m (panchromatic). The 100 m TIRS data are registered to the OLI dataset to create radiometrically, geometrically, and terrain-corrected 12-bit data products. Images were acquired from 2017 to 2020. In this study, all satellite data were downloaded from the United States Geological Survey (USGS). The image was acquired (less than 10% cloud cover) in growing stage of the specific crop cycle (dry season irrigated rice).

Digital Image Preprocessing
All satellite images were first processed by resampling the band resolution at 30 m and then mosaicked and masked. Subsequently, an algebraic raster operation and a radiometric calibration as well as geometric and atmospheric corrections were applied to the remote sensing images using ArcGIS 10.4 ® . Image acquisition was performed for each band. After that, all selected bands were converted to a 30 m resolution using a resampling technique to ensure a similar cell size and data uniformity. The average reflectance values were calculated for the study area in each band using the raster calculator tool to compensate for the spatial variability to minimize the bias. Three different blocks collected during related time periods were mosaicked to cover the large study area.

Criteria Aggregation for Land Fertility Assessment
Three vital layers were developed by satellite data, and five vegetation indices were used as the input parameters to develop the land suitability analysis (Figure 3). The most important land and soil affiliated indices, SAVI, SARVI, ARVI, MASAVI, and OSAVI, were utilized to describe land conditions (Table 1). In addition, the slope and elevation layers were used as representative land conditions. The land surface temperature was taken as another influential criterion of agriculture land suitability assessment.  Elevation is an important factor that plays a vital role in the variability of plant cover and causes temperature changes, particularly in highland areas. Areas with higher topographic elevations are more affected by rainfall and soil erosion [20]. Soil erosion is the alarming condition of agriculture field crop. Also, it is the main problems of agricultural development, such as landslides and flood events; these disasters have been severely influenced by the soil erosion process [21]. Most of the study area is plains land, and the elevation is less than 131 m (Figure 3a). The elevation data were extracted using a data elevation model (DEM) and were downscaled to a 30 m resolution.

Slope
Slope is a vital topographic element for indicating suitable farmland in assessments. Slope indicates many landscape processes, such as soil water content, erosion potential, runoff, and surface and subsurface flow velocity. The thickness of the soil layer decreases with the increasing slope [22]. The slope gradient has an impact on the runoff and soil loss: the greater the slope gradient the higher the potential for runoff and soil loss [23]. The slope was developed by using data from the original Shuttle Radar Topography Mission (SRTM) and the digital elevation model (DEM). The DEM was downscaled to a 30 m resolution. The Universal Transverse Mercator (UTM) projection system and WGS84 datum were used as rectifying agents in ArcGIS. The slope was calculated from the maximum rate of change between each cell and its neighbors. Every cell in the output raster had a slope value. Field crops generally require flat land; only a slight slope between 0% and 8% is resistant to erosion [24]. When the slope gradient is very steep (40%), soil sediment losses remain at the same high levels after cultivation abandonment because slope gradient is the main factor controlling soil erosion [25]. In the study area, the slope range was under 10% (Figure 3b), which is a suitable condition for most farming practices [26][27][28].

Land Surface Temperature (LST)
The LST (Figure 3c) was calculated for selected land areas using temporal information from Landsat 8 OLI images when less cloud coverage was present [29]. From 2017-2020, the LST data received from the obtained images ranged between 17 • C and 33 • C. The LST calculation was based on the moving average method. A single raster was formed from multiple years of raster datasets as the multiple predictions' raster.
LST was calculated for the cropland using temporal information from Landsat 8 OLI that was selected during period of less cloud coverage. Two steps were required to calculate the LST; first, the NDVI was calculated for the given time period.
After that, the resulting NDVI value was used to analyze the proportion vegetation (PV), which can be expressed as follows: After calculating the PV, the land surface emissivity (ε λ ) could be expressed as follows [30]: Second, the thermal bands are included in band 10 and band 11 from the Landsat 8 imagery. The thermal bands were converted to digital numbers to estimate the radiance. The spectral radiance could be expressed as follows: where Lλ is the TOA spectral radiance at the sensor aperture, ML is the band-specific multiplicative rescaling factor from the metadata, Q CAL is the quantized and calibrated standard product pixel value (DN), and AL is the band-specific additive rescaling factor from the metadata. Then, the brightness temperature (BT) could be expressed as follows (Jesus and Santana, 2017): where BT is the satellite brightness temperature [Celsius]; K2 is the calibration constant 2 [Kelvin], where the band-specific thermal conversion constant is taken from the metadata; and K1 is the calibration constant 1 [Kelvin], where the band-specific thermal conversion constant is taken from the metadata. Finally, LST was calculated and expressed as follows [30]: where λ is the average wavelength of band 10; ελ is the emissivity calculated from Equation (3); and ρ is (h * c σ ), which is equal to 1.438 × 10 −2 mK, where σ is the Boltzmann constant (1.38 × 10 −23 J/K), h is Plank's constant (6.626 × 10 −34 J.s) and c is the velocity of light (3 × 108 m/s).

Soil-Adjusted Vegetation Index (SAVI)
Soil has a unique spectral signature that differs from that of other types of land cover. In the visible and near-infrared zones, reflectance increases in proportion to increases in wavelength. However, the rate of increase is affected by various factors. Soil moisture and organic matter may lower the soil reflectance. The association between red and nearinfrared reflectance remains constant for different soil type physiognomies. When the moisture content changes, the two values are related and have a linear relationship. This relationship is very specific for each type of soil. SAVI is therefore useful for monitoring soils. Furthermore, SAVI is a modification of the normalized difference vegetation index (NDVI), which corrects for the influence of soil brightness when the vegetation cover is low [31]. SAVI (Figure 3d) was extracted from Landsat 8 OLI imagery by a mask for the study area. Datasets were acquired from 2017 to 2020. These datasets were used to build a single raster using map algebra in the ArcGIS platform.
To reduce the soil background effect, modified indices were proposed using the soil adjustment factor L to account for first-order soil background variations and obtain the SAVI [32]. SAVI can be expressed as follows: where ρ N IR is the reflectance value in the near-infrared band, ρ RED is the reflectance value in the red band, and L is the soil brightness correction factor. An L value of 0.5 in the reflectance space was found to minimize soil brightness variations and eliminate the need for additional calibration for different soils. The described SAVI value was 0.798 to −0.302 for the study area ( Figure 3d).

Atmospherically Resistant Vegetation Index (ARVI)
The ARVI (Figure 3e) is obtained using a self-correction process for the atmospheric effect on the red channel, using the difference in the radiance between the blue and red channels to correct the radiance in the red channel due to the excellent atmospheric resistance of the ARVI [33].
where γ depends on the aerosol type. A good default value is γ = 1 when the aerosol model is not available. ARVI is resistant to atmospheric effects due to its self-correction process. This index uses the difference in the radiance between the blue and red bands to correct the radiance in the red band. Simulations show that ARVI has a similar dynamic range as SAVI, but on average, its sensitivity to atmospheric effects is four times less than that of NDVI [34]. The ARVI value fluctuated between 0.886 and −0.662 (Figure 3e).

Soil Adjusted and Atmospherically Resistant Vegetation Index (SARVI)
SARVI has a similar dynamic range to NDVI but is four times less sensitive to atmospheric effects than NDVI. SARVI performs much better for moderate to small sized aerosol particles (e.g., continental, urban, or smoke aerosol) than for large particles. Consequently, a single combination of blue and red channels in ARVI calculations may be used in all or most remote sensing applications [34].
where L is a correction factor similar to that in the SAVI calculation and γ is similar to that in the ARVI calculation. SARVI can minimize both canopy background and atmospheric effects [35,36]. In this research, the SARVI value was found to range from 0.679 to −0.397 (Figure 3f).

Modified Soil-Adjusted Vegetation Index (MSAVI)
Richardson and Wiegand (1977) proposed a modified secondary soil-adjusted vegetation index (MSAVI) [37], which can be expressed as follows: MSAVI does not rely on the soil line principle and has a simpler algorithm; it is mainly used in soil organic matter analysis, drought monitoring, and soil erosion analysis. In addition, it is useful for plant growth analyses, desertification studies, grassland yield estimations, and leaf area index (LAI) assessments. In the study area, the MSAVI value was observed to be between +1 and −1 (Figure 3g).

Optimized Soil-Adjusted Vegetation Index (OSAVI)
The Optimized Soil Adjusted Vegetation Index (OSAVI) is a newly developed alternative that can accommodate greater variability due to high soil background values [38]. OSAVI does not depend on the soil line and can eliminate the influence of the soil background effectively. However, the applications of OSAVI are not extensive; it is mainly used for the calculation of aboveground biomass, leaf nitrogen content, and chlorophyll content [39].
where X = 1.6. OSAVI is mainly used for the calculation of aboveground biomass, leaf nitrogen content, chlorophyll content, etc. The observed value was between 0.531 and −0.201 (Figure 3h).

Pattern Analysis
Several criteria were used for pattern analysis, which required pattern analyses from multiple years or months of data to form a predicted raster for reclassification. The single raster-based calculation was not reliable, nor did it provide enough datasets. This section shows several criteria from multiple years of data (2017 to 2020) for building a predicted raster for LST, SAVI, ARVI, SARVI, MSAVI, and OSAVI. The following section introduces the pattern analysis for multitemporal datasets into a single raster.

Moving Average
The moving average was processed after completing the digital image processing steps. The moving average was calculated in each year and can be expressed as follows: where D is the number of data points in the raster cell and n is the amount of data to average.

Multiple Predicted Raster
As a part of the point pattern analysis, a single predicted raster was made. After that, LST, SAVI, ARVI, SARVI, MSAVI, and OSAVI were computed from 2017-2020. The basic extent encompassed the overall density pattern. This is basically the ratio of the observed number of single predicted rasters of points (r) to the study region area (a) and referred to as the multiple predicted raster (MPR). The MPR was applied as a criterion for land suitability analysis. MPR can be expressed as follows: where r is the ratio of the observed number of single predicted raster points and a is the area of the study region.

Land Use/Land Cover
Land use data enable the estimation of an area's coverage with vegetated areas, settlements, forests, and water bodies. Land use data were collected from the Survey of Bangladesh (SoB), which was split into 92 blocks. After aggregation in the ArcGIS platform, the data were used to develop a more accurate land use/land cover (LULC) map for the land fertility assessment. In this study, rivers, forests, water bodies, and settlements were considered restrictions in the analysis. Subsequently, after excluding the constraints, only agricultural land was considered for land evaluation. Agricultural land was subclassified into cultivated land (80%), uncultivated land (0.5%) and vegetated land (19%) (Figure 2).

Land Fertility Assessment
The weighted overlay was used to prioritize the weights of each criterion to generate a land fertility assessment map. The reclassified raster data layers were combined with the weighted overlay in ArcGIS ® . First, the combination criteria were input as equally weighted linear combinations. Second, the land suitability analysis was carried out by a fuzzy membership function, fuzzy reclassification and fuzzy overlay to evaluate the consistency of the two outcomes.

Map Development by Weighted-Linear Combination
First, reclassification was conducted to interpret the raster data by substituting a single value as the new value or by categorizing the ranges of values into a single value. Each criteria source map was reclassified into four classifications ( Table 2). Land suitability analysis was conducted using different classification categories proposed by the FAO. As suggested by the FAO's framework for land evaluation, the first class was designated as suitable (S) or not suitable (N). The suitability classification aimed to show the suitability of each land unit for crop production. In practice, three classes, namely (S1), (S2), and (S3), are typically used to identify land that is highly suitable, moderately suitable, and marginally suitable, respectively. The analysis was built using the aforementioned criteria and reclassified into 4 classes. Finally, the classes were found based on their weighted linear combination.
where C i is each criterion (i) that has been reclassified and W n is the number of data (n) that were weighted. The fuzzification process had no sharply defined boundaries that characterized the imprecision of the classes. In this process, each value of the phenomenon central to the core of the definition of the set was set to 1. Those values that were not part of the set were set to 0. Those values that fell between these two extremes were within the transitional zone of the set, which was defined as the boundary [8,[53][54][55]. In the present study, fuzzy membership classification was used to accommodate the high uncertainty of the scoring methods in assigning the suitability classes; several fuzzy membership functions were used for normalization. For this research, fuzzy functions were determined based on references and a literature review (Table 3). Out of the seven varieties of fuzzy membership functions available, three fuzzy functions were used in this study considering ecological criteria: the fuzzy small, Gaussian, and fuzzy linear functions. These functions generate continuous fuzzy classifications under standardized criteria. The reclassification tool in ArcGIS enables the transformation of categorical data to range from 0 to 10 and then divides the resulting transformed data by 10 to derive a 0 to 1 scale. The equations for the fuzzy small (Equation (15)), fuzzy linear (Equation (16)), and Gaussian functions (Equation (17)) are as follows: The fuzzy small transformation function was used when small input values were more likely to be members of the set. The criterion slope was followed by a fuzzy small function in this research.
The layers of SAVI, ARVI, SARVI, MSAVI, and OSAVI were each followed by a fuzzy linear transformation function that related a linear function between the user-specified minimum and maximum values for reclassification.
The fuzzy Gaussian function converts the primal values into a normal distribution. If the input values decrease in membership, they move away from the midpoint. The midpoint of the fuzzy Gaussian function was set to 1 [59]. The elevation and LST layers were analyzed under the fuzzy Gaussian membership function.
In the fuzzy small and fuzzy linear membership functions, the control point included a midpoint (f2) and a spread (f1). The midpoint was a specific point that had a 0.5 value of membership in the large and small functions. Gaussian functions are determined by the user based in the references (ESRI, CA, USA). The spread was generally allocated a number between 1 and 10. The fuzzy membership curve became steeper for higher spread values. The fuzzy linear transformation function applied a linear function between the minimum and maximum values. Any value below the minimum was determined to be 0 (not a member), and any value above the maximum was 1 (a member) [61]. To analyze the relationships and interactions between all the sets for the 8 criteria in the overlay model, fuzzy overlay techniques were used. The available fuzzy set overlay techniques in ArcGIS are fuzzy And, fuzzy Or, fuzzy Product, fuzzy Sum, and fuzzy Gamma. Each of these techniques described the cell's membership related to the input sets. In this study, the fuzzy Gamma overlay assisted in developing suitability maps for three identical seasons, which were determined based on references and a literature review (ESRI, Redlands, CA, USA).

Validation Using Ground Truth Data
The detected fertile zone was verified by ground reference data. The time series datasets played a vital role in developing and validating the yield prediction models. In Bangladesh, nearly 80% of the total land is allocated solely to rice cultivation [62]. Additionally, in the northwestern part of the country, approximately 70% [63][64][65] of the land is cultivated with dry season irrigated rice (boro rice). To facilitate further analysis, the major rice crop was carefully chosen for approval. The suitable area was verified by ground truth yield data. The yield data of dry season irrigated rice were collected from the Department of Agricultural Extension (DAE), Ministry of Agriculture, Bangladesh, for the 36 subdistricts in 2017-2020 to evaluate the accuracy of the land suitability analysis (Figures 4 and 5). After preparing the data in Microsoft Excel, the correlations among the five selected indices were evaluated.

Yield Prediction and Analysis
The performances of the yield prediction models were examined by field data. After the correlations among the five selected vegetation indices were established, the yield map was developed. Simple and multiple regression analyses were carried out between the mean values of the vegetation indices and the ground referenced data of the dry season irrigated rice to determine the best-fitted models for rice production. These data were classified to evaluate the production that occurred between 2017 and 2020. The SAVI, ARVI, SARVI, MSAVI, and OSAVI values were aggregated into a time series pattern (Appendix A). The yield data were compared through regression. The five vegetation index values were collected from reference points in the 36 subdistricts. Ground truth data information was obtained from the 36 subdistricts, and the yield was reported in metric tons per hector (MT/ha).

Land Suitability Analysis
The weighted linear model was used to prioritize the weights of each criterion to generate the land suitability assessment. First, the variables were analyzed as equally weighted linear combinations. Second, the suitability assessment was carried out by a fuzzy membership function to verify the consistency of the two procedure results ( Table 4). The land fertility analysis (Figure 6a) with equal weights showed that 43% of the land (1832 km 2 ) was highly suitable, 41% of the land (1747 km 2 ) was moderately suitable, and 10% of the land (426 km 2 ) was marginally suitable. In addition, the restricted zone was defined as an unsuitable area. In this research, the unsuitable area was found to cover 6% (256 km 2 ). However, the land suitability analysis using the fuzzy membership function (Figure 6b) showed that 48% (2045 km 2 ) of the land area consisted of the most suitable area, 39% of the land (1661 km 2 ) was moderately suitable, and 7% of the land (298 km 2 ) was marginally suitable. In addition, restricted areas accounted for 6% of the land area. In fuzzy overlay analysis, 256 km 2 of the area was classified as fallow land that is not suitable for cultivation.

Yield Estimation
The predictors derived from the satellite imagery in the form of spectral bands or vegetation indices (SAVI, SAVI, ARVI, SARVI, MSAVI, and OSAVI) were the most effective spectral parameters for predicting rice yield (Figure 7). In addition, the individual index values were extracted from the ground truth information (Figures 4 and 5) for the subdistricts that were located in the highly suitable areas. A trend line approach was used to verify the index influences at different study points (Figure 8). Regression analysis was performed between the vegetation indices and the observed yield. The SAVI, SARVI, ARVI, MSAVI, and OSAVI values were obtained from satellite imagery from 2017-2020 (Table 5). The results showed good accuracy in the regression analysis using SAVI (R 2 = 77.3%), ARVI (R 2 = 68.9%), SARVI (R 2 = 71.1%), MSAVI (R 2 = 74.5%) and OSAVI (R 2 = 81.2%) (Figure 9). From the multiple regression model, it was observed that using more than one variable for the yield prediction increased the model accuracy by enhancing the R 2 value. However, the best-fitted models were obtained using the SAVI-ARVI-SARVI-MSAVI-OSAVI composite vegetation index. The yield prediction model with the composite index had a goodness of fit of R 2 = 0.839. The model was used to estimate the yield in the time series dataset (Table 6).     The developed yield map indicated that in 2017, the maximum yield was 4.59 MT/ha ( Figure 10). Furthermore, in 2018 and in 2019, it was 4.9 MT/ha and 5.08 MT/ha, respectively. For 2020, the predicted yield range appeared to be between 0.269 MT/ha and 4.537 MT/ha.

Discussion
This research provided a comprehensive strategy to create agricultural land use plans for cultivation considering suitable conditions, which were derived from satellite remote sensing data. Previously, much of the research conducted only reported land suitability for site-specific plans or single-cropping plans [66,67]. However, this research attempted to develop an overall land suitability assessment using soil-vegetation representative variables that extracted only satellite remote sensing data from the GIS platform when field soil sampling is inconvenient and expensive. Applying only remote satellite datasets to assess suitable land conditions was a source of concern that added a new feature to MCDM-land suitability analyses. In this study, the reliability of five vegetation indices was verified by a regression analysis that incorporated ground truth yield data, and the results were used for yield map preparation. Vegetation phenology analyses have potential [11,14] in estimating yield prediction with good accuracy in highly suitable areas. In addition, two topological factors (slope and elevation) and another environmental parameter, land surface temperature, were extracted from the USGS, which ensured better performance of the results in land use planning.
Either the AHP-based weighted overlay or equal-overlay technique is typically applied in most studies [68][69][70]; few studies have conducted fuzzy membership methods employed in the GIS platform by incorporation with the AHP technique [69,71,72]. Some studies have conducted farmland assessments based on soil testing [73,74]; however, considering the preparation of a suitable map for agricultural land using soil-represented remote sensing data for the linear combination of the F-MCDM approach is a new dimension of this research. A multicriteria decision-making system was applied to reduce the biases in the analysis. Variation in the land surface temperature was an important factor in this area and influenced the locations considered most suitable for crop cultivation. Moreover, atmospherically restricted vegetation indices (ARVI) and soil adjusted atmospherically restricted vegetation indices (SARVI) were used to reduce the biases associated with atmospheric effects.
Most of the suitable lands were located in the northern part, and marginally suitable lands were mostly located in the northwestern part; this result was likely due to the influence of high elevation. In addition, unsuitable zones were found mostly in the eastern part due to the presence of water bodies that are not arable for cultivation along with other adverse edaphic factors. Previous studies had the limitation of obtaining inappropriate validation results due to inadequate ground reference information. In this research, the validation of the results was accomplished by physical verification with the corresponding time series yield data of the most cultivated crop, dry season irrigated rice, which usually grows over 70% [65] of the agricultural land area. The suitable conditions were not verified by the other crop yield data, which was the main limitation of this research.

Conclusions
This research established a method to identify the most suitable agricultural land by using the potentiality of satellite remote sensing data integrated with weighted linear combinations and fuzzy multicriteria analyses in a GIS platform. The multicriteria decision analysis was performed for suitability assessment using eight criteria: elevation, slope, and LST vegetation indices (SAVI, ARVI, SARVI, MSAVI, and OSAVI). To derive a more accurate result, a land use/land cover layer was also used to mask restricted zones. The land suitability analysis with equal weights showed that 43% of the land (1832 km 2 ) was highly suitable, 41% of the land (1747 km 2 ) was moderately suitable, and 10% of the land (426 km 2 ) was marginally suitable. Conversely, expert knowledge was also considered, along with consistent assessments when using the fuzzy membership function; 48% of the land (2045 km 2 ) was highly suitable, 39% of the land (2045 km 2 ) was moderately suitable, and 7% of the land (298 km 2 ) was marginally suitable. The yield estimation using SAVI (R 2 = 77.3%), ARVI (R 2 = 68.9%), SARVI (R 2 = 71.1%), MSAVI (R 2 = 74.5%), and OSAVI (R 2 = 81.2%) showed a good accuracy. In addition, every combination of these five indices represented the best accuracy (R 2 = 0.839), which was used to develop the yield maps for the corresponding years (2017-2020). The results of the land suitability evaluation for field crops will be very useful in the decision-making process to increase production as well as for the sustainable management of agricultural lands. Thus, the influence of vegetation index evaluations, suitable condition assessments, and yield prediction models is essential for understanding future land use and production trends in the agricultural crop sector in Bangladesh, as well as other applications.

Data Availability Statement:
The data used is primarily reflected in the article. Other relevant data is available from the authors upon request.

Acknowledgments:
We would like to express their sincere gratitude to the Graduate School of Life and Environmental Sciences, the University of Tsukuba, for supporting this research in land fertility assessment and index evolution using satellite datasets. We are sincerely thankful to the United States Geological Survey (USGS) for providing Landsat 8 OLI and Shuttle Radar Topography Mission (SRTM) data free of cost. We are grateful to the Department of Agricultural Extension (DAE), Ministry of Agriculture, Bangladesh, for their cordial support by providing ground reference yield data for 36 subdistricts from 2017-2020. Furthermore, we are also thankful to the Japan International Research Center for Agricultural Sciences (JIRCAS) for technical support to pursue this research in Japan.

Conflicts of Interest:
The authors declare no conflict of interest.