Nothing Special   »   [go: up one dir, main page]
Next Article in Journal
Isolating Anthropogenic Wetland Loss by Concurrently Tracking Inundation and Land Cover Disturbance across the Mid-Atlantic Region, U.S.
Previous Article in Journal
Using Satellite Interferometry to Infer Landslide Sliding Surface Depth and Geometry
Previous Article in Special Issue
Urban Land Use and Land Cover Classification Using Multisource Remote Sensing Images and Social Media Data
You seem to have javascript disabled. Please note that many of the page functionalities won't work as expected without javascript enabled.
 
 
Search for Articles:
Title / Keyword
Author / Affiliation / Email
Journal
Article Type
 
 
Advanced Search
 
Section
Special Issue
Volume
Issue
Number
Page
 
Journals
Remote Sensing
Volume 12
Issue 9
10.3390/rs12091463
remotesensing-logo
Submit to this Journal Review for this Journal Propose a Special Issue
► Article Menu

Article Menu

share Share announcement Help format_quote Cite question_answer Discuss in SciProfiles thumb_up
...
Endorse textsms
...
Comment
first_page
settings
Order Article Reprints
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Optimal Soybean (Glycine max L.) Land Suitability Using GIS-Based Multicriteria Analysis and Sentinel-2 Multitemporal Images

by
Dorijan Radočaj
1,*,
Mladen Jurišić
1,
Mateo Gašparović
2 and
Ivan Plaščak
1
1
Faculty of Agrobiotechnical Sciences Osijek, Josip Juraj Strossmayer University of Osijek, Vladimira Preloga 1, 31000 Osijek, Croatia
2
Faculty of Geodesy, University of Zagreb, Kačićeva 26, 10000 Zagreb, Croatia
*
Author to whom correspondence should be addressed.
Remote Sens. 2020, 12(9), 1463; https://doi.org/10.3390/rs12091463
Submission received: 24 March 2020 / Revised: 3 May 2020 / Accepted: 4 May 2020 / Published: 5 May 2020
(This article belongs to the Special Issue Applications of RS and GIS Integration in Natural Resources and Environmental Science)
Browse Figures
Graphical abstract
"> Figure 1
<p>The workflow of soybean land suitability calculation and evaluation.</p> "> Figure 2
<p>Study area: (<b>a</b>) location of Osijek-Baranja County, (<b>b</b>) location of used weather stations and soil samples, (<b>c</b>) location of centroids of test soybean parcels.</p> "> Figure 3
<p>Weather data at the sensing time of four used Sentinel-2 images from weather station Osijek provided by the Croatian Meteorological and Hydrological Service (DHMZ).</p> "> Figure 4
<p>Preprocessed criteria rasters.</p> "> Figure 5
<p>Graphs of standardization parameters per criterion for evaluated standardization methods.</p> "> Figure 6
<p>Regression graphs for best fitting functions per the suitability model.</p> "> Figure 7
<p>Final suitability maps for soybean cultivation for calculated models.</p> ">
Versions Notes

Abstract

:
Soybean is regarded as one of the most produced crops in the world, presenting a source of high-quality protein for human and animal diets. The general objective of the study was to determine the optimal soybean land suitability and conduct its mapping based on the multicriteria analysis. The multicriteria analysis was based on Geographic Information System (GIS) and Analytic Hierarchy Process (AHP) integration, using Sentinel-2 multitemporal images for suitability validation. The study area covered Osijek-Baranja County, a 4155 km2 area located in eastern Croatia. Three criteria standardization methods (fuzzy, stepwise and linear) were evaluated for soybean land suitability calculation. The delineation of soybean land suitability classes was performed by k-means unsupervised classification. An independent accuracy assessment of calculated suitability values was performed by a novel approach with peak Normalized Difference Vegetation Index (NDVI) values, derived from four Sentinel-2 multispectral satellite images. Fuzzy standardization with the combination of soil and climate criteria produced the most accurate suitability values, having the top coefficient of determination of 0.8438. A total of 14.5% of the study area (602 km2) was determined as the most suitable class for soybean cultivation based on k-means classification results, while 64.3% resulted in some degree of suitability.

Graphical Abstract">

Graphical Abstract

1. Introduction

The increasing demand for food and bioenergy in the world stimulated the development of land suitability calculation methods, as a basis for effective agricultural land management and environmental sustainability [1,2]. Conventional agriculture is characterized by a very high input of fertilizer, pesticides and herbicides, having a negative long-term impact on sustainability [3]. The selection of naturally suitable areas for the particular crop type cultivation reduces the overall application of inputs, creating an optimal environment for crop growth [4]. Agricultural land management plans based on inappropriate evaluation of natural resources limit crop yields and increase production costs [5]. Chen et al. [6] recommended the development of crop-specific evaluation indices through land suitability determination to ensure the sustainability of agricultural production. Delineation of suitability classes with homogenous characteristics is a fundamental segment of land evaluation, allowing effective implementation of land use planning in the field [7]. Mapping of such suitability classes for particular crop cultivation is essential for the transfer of knowledge to the end-users, whether to land management experts or individual farmers and farming companies [8].
Multicriteria analysis is widely recognized as the method for the selection of the most suitable (optimal) location and its alternatives in various areas, such as agriculture [9,10,11], forestry [12], land management [13] and environmental planning [14]. With the integration of spatial components from Geographic Information System (GIS), an approach of GIS-based multicriteria analysis enables a suitability modeling for any entity related to space [15]. The benefit of GIS-based multicriteria analysis in agriculture is its universal applicability regarding crop types, area size and location in the world. Its potential is conditioned only by the crop expert’s knowledge of the optimal agroecological conditions of the selected crop type and the quality of the input spatial data [16,17]. The core procedures of multicriteria analysis are based on the establishment of the relationship between relevant criteria [18]. Many methods have been developed to determine the relative relationship between criteria for the suitability calculation, one of the most notable being the Analytic Hierarchy Process (AHP) [19]. Standardization of criteria values is a necessary procedure in the calculation of suitability models [20]. Input values of all criteria are commonly transformed into a unique number interval during the standardization process for further processing. Many authors noted an impact of the selection of standardization methods in multicriteria analysis on land suitability values [21,22,23,24]. Standardization using simple linear scale transformation was usually the selected method in these studies. The hybridization of GIS-based multicriteria analysis with unsupervised classification presents a novel approach in the management of calculated suitability values, enabling the effective creation of land suitability classes [25]. The delineation of suitability classes is regarded as a necessary procedure in suitability analyses and precision agriculture, as it significantly facilitates the application of GIS-based multicriteria analysis results in the field [26]. Van Niekerk noted the superiority of computer algorithms over traditional manual mapping techniques for the delineation of suitability classes, as they allow objective and time-efficient classification [27].
Soybean is a fundamental component of agricultural land management plans worldwide, presenting a major source of protein for humans and a high-quality animal feed [28]. According to the Food and Agriculture Organization of the United Nations (FAO) publication [29], soybean accounted for about one-third of the total harvested area devoted to annual and perennial crops, while its share in global oil-crop output was 44%. According to the projection by the European Commission for the period between 2019 and 2030, the production of soybean food products will continue to grow due to demand for locally produced plant-protein food [30]. The same source stated that the soybean area would show significant land-use change, resulting in a 5% harvested area increase by 2030. The long-term projection from the United States Department of Agriculture to 2029 predicts continuous global demand for soybean oil for biodiesel production [31]. The competitive position of soybean among arable crops has steadily improved due to consistent improvements in yield and reductions in production costs [29]. However, additional investment for soybean yield improving research was urged to land policymakers and managers [32].
Medium-resolution multispectral satellite imagery presents an important data source for the observation of crop characteristics in yield improving research and mapping for agricultural land management [33,34]. Sentinel-2 is a multispectral satellite mission from the European Space Agency Copernicus program, started in 2015 [35]. The constellation of Sentinel-2 mission is based on two satellites, Sentinel-2A (S2A) and Sentinel-2B (S2B), orbiting 180° apart [36]. Sentinel-2 provides the possibility of an effective crop monitoring, having a 290 km swath width, high imaging resolution (up to 10 m) and revisit time of two to three days at mid-latitudes [36]. Crop parameters with the application of remote sensing are commonly monitored using vegetation indices [37]. The most used vegetation index is the Normalized Difference Vegetation Index (NDVI), which enables the determination of the relationship between photosynthetic and optical properties of crops [38]. NDVI is the most widely used indicator in studies regarding crop biomass and chlorophyll content [39]. Satellite-derived vegetation indices have been successfully used as a part of various crop models, such as the estimation of maximum evapotranspiration and irrigation requirements [40] and multitemporal crop monitoring [38].
The general objective of the study was the determination of optimal soybean land suitability and its mapping based on GIS-based multicriteria analysis. For this purpose, the evaluation of the three most commonly used standardization methods in recent GIS-based multicriteria analysis studies was conducted. The validation of calculated suitability values was performed using a novel peak NDVI method derived from Sentinel-2 multitemporal images as an independent accuracy assessment.

2. Materials and Methods

The proposed method of GIS-based multicriteria analysis for soybean suitability calculation and its validation consists of six major steps (Figure 1). These are the selection and preprocessing of relevant criteria, data standardization, weight determination, criteria aggregation, validation of calculated suitability models and creation of final suitability maps. Open-source GIS software was used in this research: SAGA GIS v7.4.0 (Hamburg, Germany) [41] for data preprocessing and calculations, QGIS v3.8.3 (Grüt, Switzerland) [42] for data visualization and GRASS GIS v7.8.2 (Bonn, Germany) [43] for the calculation of solar irradiation.

2.1. Study Area

The study area covers Osijek-Baranja County, a 4155 km2 area located in eastern Croatia (Figure 2). According to the Croatian Bureau of Statistics [44], Osijek-Baranja County has the largest utilized agricultural area in Croatia, covering 13.6% of total country agricultural land in 2016. Corine 2018 Land Cover data showed that agricultural areas are the dominant land cover class in the county, covering 64.5% of the county’s area. Forests are the second-largest class covering 26.7% of the county’s area, followed by artificial surfaces (3.8%), wetlands (2.7%) and water bodies (2.3%). Croatia follows a world trend of an increase in soybean cultivation. According to [45], a total of 111,316 t was produced in 2013 and 244,075 t in 2016, a 119% soybean production increase during four years. Soybean is the fourth highest cultivated crop in the study area behind maize, wheat and sunflower, covering 8.54% of the total agricultural land in the county during 2019 per Croatian Paying Agency for Agriculture, Fisheries and Rural Development data. Osijek-Baranja County has a moderately warm, rainy climate on a Köppen scale, according to the same source. Early soybean variants were commonly sown in eastern Croatia in recent years [46]. Sowing of these variants was conducted in late April, while harvest typically occurred mid-October. The Croatian agricultural and forestry advisory service recommends the application of 30 kg/ha of nitrogen (N), 60 kg/ha of phosphorus pentoxide (P2O5) and 90 kg/ha of potassium oxide (K2O) during basic fertilization in autumn and minor adjustments before sowing [47]. Irrigation of soybean fields has been extremely overlooked in the study area, with only 64 ha employing some sort of irrigation system, according to Croatian Paying Agency for Agriculture, Fisheries and Rural Development data for 2019.

2.2. Selection and Preprocessing of Relevant Criteria

Climate and soil criteria were considered to have the most significant impact on soybean cultivation, based on previous research regarding soybean land suitability [48,49,50,51,52]. Relevant climate and soil criteria selected in this study are shown in Table 1. The selected number of criteria in each criteria group is six, which meets the specifications for the application of the AHP method for the consistency of information between criteria [53]. The main data sources for the modeling of selected criteria were: Croatian Meteorological and Hydrological Service (DHMZ) data during the soybean growth period between April and October for years 2015–2019; Croatian Agency for the Environment and Nature (CAEN) soil samples collected during 2016; a basic pedologic map of Croatia; Shuttle Radar Topography Mission (SRTM) and Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) global digital elevation models.
Climate criteria consisting of Tmin, Tavg, Tmax, Precipitation and AirHumidity were interpolated from 15 DHMZ weather station data. Soil samples for modeling of soil criteria were downloaded from HAOP Web Feature Service. Samples were collected in the field by Croatian Geological Survey and Croatian Agency for Agricultural Land according to standards ISO 10390:2005 for soil pH, ISO 11277:2009 for soil texture and ISO 11466:1995 for soil C/N. From 72 regularly distributed samples in the study area, 48 were detected as agricultural area land cover and filtered for further processing. The solar irradiation criterion was defined as Global Horizontal Irradiation (GHI), according to previous research [54]. GHI was calculated using the novel method developed by Gašparović et al. [55], based on the ASTER digital elevation model, Linke turbidity factor [56] and effective cloud albedo acquired from Meteosat Second Generation satellites SARAH Edition 2 [57]. According to Gašparović et al. [55], GHI produced 3% higher accuracy than commercial solar irradiation data, compared to the data measured from ground stations. The same approach for GHI calculation was successfully applied in similar research [14]. The final GHI used in this research was an average daily mean value (kWh/m2) between April and October in 2013 to 2015. A three-year period was used to reduce the influence of weather conditions. The time frame between April and October covers the vegetation period of all major soybean varieties cultivated in the study area. Tmin represents a mean minimum air temperature for April and May, as soybean is most susceptible to frost at early vegetation stages that occur during April and May. Tmax represents mean maximum air temperatures in June, July and August and refers to the drought risk. Soil type classes were derived from a basic pedologic map of Croatia. Soil pH was selected as a criterion due to soybean demands for neutral soil, with neutral and mildly acidic soils having optimal values. Soybean does not have strict demands regarding soil C/N [58], but it has a major impact on sustainable planning of overall agricultural production in the study area. Soil C/N deficiency presents a major issue in the study area [59], so crop types with larger demands for C/N would benefit more for being cultivated on soil with higher soil organic content. Soil texture classification was conducted in twelve classes from the soil texture triangle, according to silt, sand and clay soil content [60]. Loam was considered as the most suitable soil texture, allowing the optimal root system development for soybean [61]. The process of soil texture raster generation from silt, sand and clay soil content rasters was automatized using Python v3.7.4 (Wilmington, Delaware, United States of America) [62]. The slope was determined using SRTM 1-arc second global digital elevation model. The slope is associated with the exploitation of air humidity in soybean cultivation, which is important for its development in reproductive growth stages. Terrain slope values were considered as hilly terrain restricts the adequate exploitation of agricultural machinery. Topographic wetness index (TWI) was calculated in SAGA GIS for the determination of specific catchment areas, representing the effect of inclined slopes on the soil water content [63].
Preprocessing of input criteria consisted of the conversion of polygon vector criterion to a raster format, georeferencing of weather station data and conversion to point vector format, evaluation of interpolation methods for climate and soil points, interpolation of these values using the most accurate method and clipping of rasters to study area extents. The selected projection coordinate system was the Croatian Terrestrial Reference System (HTRS96/TM, EPSG: 3765) with 250 m spatial resolution (ground sample distance), so all rasters were reprojected and resampled accordingly. Interpolation methods selected for interpolation accuracy assessment were Ordinary Kriging (OK) [64], Inverse Distance Weighting (IDW) [65] and Angular Distance Weighting (ADW) [66]. OK and IDW were selected as they resulted as the most accurate interpolation methods for the interpolation of soil parameters in studies by Yao et al. and Qiao et al. [67,68]. IDW and ADW produced the highest accuracy of tested deterministic interpolation methods for climate data in a study by Xavier et al. [69]. Descriptive statistics (Table 2) were calculated for the determination of input samples normality and stationarity with the aim of selection of optimal interpolation method [70]. The mathematical model and fitting range for OK interpolation were selected on the highest possible fitting coefficient of determination to variogram. Both IDW and ADW were performed with a weight parameter of 2. Accuracy assessment was conducted based on cross-validation with the leave-one-out procedure, with the coefficient of determination (R2), root-mean-square error (RMSE) and normalized RMSE (NRMSE) as statistical indicators. NRMSE was calculated by dividing RMSE with a respective average of measured values from input climate or soil samples.

2.3. Data Standardization

Three standardization methods were evaluated in this research: linear standardization, stepwise standardization and fuzzy standardization. The defined number intervals for standardization were closed intervals containing values from zero to one (0.00, 1.00) or one to zero (1.00, 0.00), where 0.00 designates least favorable and 1.00 designates the most favorable impact on suitability values. The integration of qualitative and quantitative criteria values was conducted in this procedure. Values from the selected number interval during standardization were designated to either thematic classes in case of qualitative criteria or input numeric values for quantitative criteria, according to the level of suitability. Linear, stepwise and fuzzy standardization were applied to all selected criteria, except for the qualitative criteria SoilType and SoilTexture. For the two qualitative criteria, stepwise standardization was selected as a standardization method, being the only method that could operate with qualitative values. Consequentially, these values were used in combination with all three applied standardization methods for quantitative criteria.
The linear standardization method, based on linear scale transformation, is the most frequently used deterministic method for standardization in GIS-based multicriteria analysis [71]. The deterministic nature of the method ensures completely objective standardization, with no crop expert’s impact on the standardization result. Linear standardization produced good results for the standardization of distances to infrastructure objects [72] and climate data for crop suitability calculation [73]. The calculation of linear standardization values was based on the score range procedure [20].
Stepwise standardization is also a commonly used standardization method in GIS-based multicriteria analysis [74]. It is based on the distinctive selection of the input values in number intervals defined by the crop expert for the standardization in new classes. An arbitrary value in the standardization interval has been assigned to each new class, allowing an extensive subjective crop expert’s impact on the result. Four to six classes were used for stepwise standardization, depending on the heterogeneity of input criteria values. Stepwise standardization was successfully applied for the standardization of soil chemical properties [75] and climate criteria [76].
Fuzzy standardization was based on the application of fuzzy membership functions. Similarly to stepwise standardization, fuzzy membership functions allowed extensive impact on the standardization result by the subjective crop expert. Three membership functions were used in this research: linear, S-shaped and J-shaped [77]. Linear membership function (2) was a base for the calculation of S-shaped (3) and J-shaped (4) membership function values. Topographic indicators, such as slope, were successfully standardized using fuzzy membership functions [78]. The same method was also applied for the standardization of soil criteria [79] and the combination of soil and climate criteria [80] in agriculture land suitability calculations. Membership function parameters were marked as a and d for the definition of support, alongside b and c for the definition of the core of the membership function [81]. The variety of selected membership functions allowed the subjective crop expert’s impact on the intensity of ascending and descending of standardized values. Intermediate values between parameter values were gradually calculated, allowing the continuous representation of the criteria effect on suitability. Equations for the calculation of the standardized values using fuzzy membership functions were derived from Novák [77].

2.4. Weight Determination

Weight determination was conducted using the AHP method in three approaches for climate criteria group, soil criteria group and combined criteria group consisting of both climate and soil criteria. Two pairwise comparison matrices were created to quantify criteria weights for climate and soil criteria, as a base of the AHP method [82]. Individual weights per criteria group were calculated using separate pairwise comparison tables. Climate and soil criteria groups were considered as equally influential in suitability result. The primary motive of using three approaches was the estimation of criteria groups’ impact on the suitability result in the validation process. This procedure allowed the criteria modifications for further application, through the elimination of redundant criteria and the addition of new criteria. AHP is regarded as one of the most advantageous weight determination methods in the multicriteria analysis by Musakwa [83], primarily due to being flexible, straightforward and comprehensive. It was successfully applied in GIS-based multicriteria analyses for crop cultivation suitability [10,84,85], irrigation suitability [76,86] and environmental analyses [14].
Each combination of two criteria per criteria group was evaluated in pairwise comparison matrices by assigning the preferred criteria value ranging from 1 (equal preference) to 9 (extreme preference) [82]. The determined values were associated with preferred criteria in the pairwise comparison matrix, while the reciprocal value was assigned to the less-preferred criteria. The consistencies of pairwise comparisons per matrix were evaluated by Consistency Ratio (CR) (1), using the Consistency Index (CI) (2) and the Random Consistency Index (RI) [87]:
C R = C I R I ,
C I = λ   −   n n   −   1 ,
where λ is the average value of consistency vector and n is the number of parameters. RI values were predefined by Saaty [87], depending on the number of compared criteria in the pairwise comparison matrix. The acceptable value of CR is 0.10, while higher values indicate that the revision and modification of the pairwise matrix should be conducted [88].

2.5. Criteria Aggregation

The weighted linear combination was the selected method for criteria aggregation in the process of soybean suitability calculation. It is the most commonly used aggregation method in the GIS-based multicriteria analyses [89]. The advantages of the selected method were simplicity and time-efficiency in the criteria aggregation, as well as the objective selection of the most suitable location and its alternatives, ranked by the suitability values [90].
The nine suitability combinations of three standardization methods and three weight determination methods were calculated, resulting in nine suitability maps. Soybean suitability index (SSI) was selected as soybean land suitability value, ranging from 0 to 1. According to Eastman [20], the equation for SSI calculation using weighted linear combination was (3):
S S I   = ∑ w i X i   ×   C
where wi denotes the weight of factor i, Xi denotes the standardized values of factor i and C denotes a Boolean layer of constraints.

2.6. Validation of Calculated Suitability Models

Calculated suitability values were validated by NDVI values derived from Sentinel-2 Level 2A images. No official database or reliable records about yield exist in Croatia, so NDVI was selected as a basis for validation. Since it was derived from Sentinel-2 images, nearly global applicability is ensured. Four images were downloaded, reprojected to HTRS96/TM and clipped to study area extents to cover the soybean seed filling growth stage (R6). This growth stage is associated with the highest total above-ground biomass during the soybean phenology cycle [91]. NDVI derived from Sentinel-2 images produced a high correlation with biomass and chlorophyll content of crops in previous research [92]. Soybean total above-ground biomass also resulted in a high correlation with yield [93]. These factors were selected as representative for validation with NDVI, as higher biomass and yield are expected on more suitable land for soybean cultivation. Not all soybean variants cultivated in the study area reach the R6 stage at the same time due to different periods of sowing, so the necessity for multitemporal images emerged. The time frame starting from 15 July to 10 August was determined for all soybean varieties in the study area to reach the R6 growth stage, based on the empirical knowledge of soybean cultivation in the study area. Four Sentinel-2 images were selected for validation from this time frame based on their availability and lowest possible cloud cover percentage (Table 3). The used images were temporally evenly distributed in the time frame of the R6 soybean growth stage and sensed in very similar meteorological conditions (Figure 3). All four NDVI images were used for the determination of peak NDVI of soybean test parcels, with slightly more than two-thirds having its peak NDVI during late July.
A total of 204 soybean test parcels, evenly distributed on the study area, were selected as the ground truth data. The total cumulative area of soybean parcels was 849.3 ha, with an average individual area of 4.16 ha. Shapefile vector data representing soybean parcels were collected from Croatian Paying Agency for Agriculture, Fisheries and Rural Development official spatial database, filtered by the soybean crop type. Four NDVI rasters derived from downloaded Sentinel-2 images were overlayed with the centroids of soybean parcels, so all NDVI values within a soybean parcel were averaged and designated to that centroid. Since spatial resolutions of suitability rasters and Sentinel-2 images were different (250 m and 10 m, respectively), centroids of reference soybean parcels with peak NDVI values were used to perform validation of suitability rasters. Centroids of soybean test parcels were forced inside polygons during creation to represent geometrically irregular parcels accurately. The comparison of these data sets was achieved with fast processing time, as point vector was overlayed with raster data. Peak NDVI values per centroid were determined as the maximum NDVI value per centroid, based on four NDVI rasters derived from multitemporal Sentinel-2 images.
The use of peak NDVI allowed the detection of the R6 growth stage for all ground truth polygons, which represented all soybean variants in the study area. Such a procedure proved to be robust regarding the atmospheric effects like clouds and haze on satellite images since validation does not depend on the individual images. Identification data of Sentinel-2 images used for the validation of suitability models are shown in Table 2.
Four regression functions were selected for the calculation of R2 between suitability values and reference peak NDVI values. Selected functions were linear (R2lin), logarithmic (R2log), exponential (R2exp), second-order polynomial (R2poly2) and third-order polynomial (R2poly3). Different regression functions were evaluated to determine the most accurate relationship model between suitability values and peak NDVI. Validation of each suitability model was based on the highest R2 calculated using these functions. RMSE and Cohen’s d index were calculated as a complementary statistical values to R2 for accuracy assessment. Cohen’s d values were interpreted according to [94]. Peak NDVI and suitability values were normalized in the (0.00, 1.00) range prior to RMSE and Cohen’s d calculation to describe suitability results while referring to the same number interval.

2.7. Creation of Final Suitability Maps

Unsupervised classification using the k-means algorithm was applied for the delineation of five suitability classes for each calculated suitability model. Created suitability classes were classified according to FAO standards in classes representing very suitable (S1), moderately suitable (S2), marginally suitable (S3), currently unsuitable (N1) and permanently unsuitable area (N2) [95]. Unsupervised classification was therefore conducted in five classes to meet the FAO specification. Ranking of class suitability values was conducted using mean suitability values of all pixels per class. The selected method approached the suitability values on a relative basis using a computer-automated classification. This procedure enabled the creation of suitability classes independently of selected criteria or study area properties. Such an approach reduced the effect of subjectivity in the pairwise comparisons of criteria in AHP, which is its major disadvantage [96]. For the application of this method in the regionalization of agricultural production, the integration of suitability values for different crop types, and consequentially from different crop experts, is necessary. By managing the suitability values using unsupervised classification, the integration was possible since the relative criteria relationship was retained. The crop experts’ judgment levels have been balanced in the process, preventing overestimation or underestimation of suitability. The applied approach also ensured the objectiveness in the creation of suitability classes, as the training data creation was not necessary, in contrast to supervised classification. Training data creation presents a common source of human error in classification due to the operator’s subjectivity in the selection of training polygons [97]. Input data in the k-means algorithm was raster representing suitability values from the most accurate suitability model. Algorithms based on unsupervised classification were successfully applied for suitability classification with a single raster input [27]. K-means was successfully applied in various combinations with GIS-based multicriteria analysis in crop suitability analyses [25,98,99] and delineation of suitability classes in precision agriculture [100]. Final suitability maps were created by the removal of constrained areas from the unsupervised classification results. CORINE 2018 Land Cover data was applied to exclude agricultural land cover, as no other land cover was available for efficient soybean cultivation. The goal of this step was to exclude permanent natural (forests, wetlands, water bodies) and built-up objects from classification in a 250 m spatial resolution, so CORINE 2018 Land Cover data was selected as the most fitting data source due to its high thematic accuracy (85%) [101]. All calculated suitability models were mutually compared using Pearson’s correlation coefficient in the correlation matrix.

3. Results

The descriptive statistics of input samples for evaluation of interpolation methods are shown in Table 3. Moderate skewness and high kurtosis deviations were observed from optimal values of zero and three, respectively. While low coefficients of variation indicated high data stationarity, low data normality was observed for both climate and soil criteria. Sand soil content was the only exception from this observation, as it produced both low normality and stationarity.
Accuracy assessment results of tested interpolation methods are shown in Table 4. Both deterministic methods outperformed OK in a case of climate data with a lower number of samples. IDW was selected as an optimal method for the interpolation of climate data, producing higher accuracy values than ADW for Tmin, Precipitation and AirHumidity. Accuracy of interpolation for Tavg and Tmax produced mixed results, as ADW produced higher R2 for Tavg and higher RMSE and NRMSE for Tmax than IDW, but lower RMSE, NRMSE and R2, respectively. All tested methods produced lower mean accuracy for soil data compared to climate data. IDW produced the highest accuracy values for the interpolation of moderately variable data: Clay, Silt and C/N. ADW produced the highest R2 for pH interpolation and second-best RMSE and NRMSE values behind IDW. OK produced slightly higher RMSE and NRMSE for interpolation of sand, but 0.0346 lower R2 than IDW. IDW was therefore also selected as an optimal interpolation method for soil data.
Soil texture classes resulting after classification were silty clay, silty clay loam, clay loam, loam and silt loam. The dominant soil texture class in the study area was silty clay loam, covering 72.1% of the total area. The next two largest classes, silty clay and silt loam, had significant silt content, covering 21.6% of the study area combined. Clay loam and loam mostly covered the transitional areas between two larger classes, resulting in 3.2% and 3.1% of the study area, respectively. All preprocessed input criteria rasters are shown in Figure 4.
Parameters for each standardization method were selected as shown in Figure 5. The presence of a few extreme input values prevented the linear standardization method to expand the full (0.00, 1.00) number interval to the majority of pixels, resulting in an inaccurate representation of these criteria. This primarily refers to slope and C/N, where the number interval covers most of the pixels in (0.85, 1.00) and (0.00, 0.25) number interval, respectively. Intermediate values for stepwise standardization were selected based on standardization classes count, alongside a minimum value of 0.00 and a maximum of 1.00. The J-shaped and linear methods were the most commonly used fuzzy membership function for standardization, both being applied for 4 out of 10 quantitative criteria. Precipitation, SolarIrradiation and AirHumidity resulted in a high variability, while Tavg, Tmin and Tmax produced moderate variability of the climate conditions in the study area. All soil criteria resulted in high variability in the study area.
Pairwise comparison matrices were created for climate (Table 5) and soil criteria (Table 6). Consistency tests for both matrices resulted in the allowed tolerance. The weight sums of interpolated criteria were 0.877 for climate and 0.431 for soil criteria group.
Calculated R2, RMSE and Cohen’s d values for the validation of suitability models are displayed in Table 7. Suitability calculated using the climate criteria group constantly produced the lowest R2, with a maximum for fuzzy standardization using exponential regression. Stepwise and fuzzy standardization models both produced the highest accuracies for soil and combined criteria groups, with combined variants being slightly more accurate. Overall, fuzzy standardization with combined criteria produced the highest R2 and RMSE, presenting the optimal model for soybean land suitability mapping. RMSE values were based on normalized values, so they resulted in (0.00, 1.00) number interval and therefore enabled objective accuracy assessment of suitability values. Cohen’s d values indicated the strongest relationship between fuzzy standardization and soil criteria with peak NDVI values. The fuzzy standardization with combined criteria also resulted in a small effect size, having a large value gap between the third-lowest d value of fuzzy standardization with climate criteria.
Visual representation of regression functions that produced the highest coefficient of determination per model is presented in Figure 6. The similarity of stepwise and fuzzy standardization with soil and combined criteria results can be observed, with fuzzy models resulting in slightly smaller dispersion. RMSE values showed similar results, as fuzzy methods produced lower RMSE than stepwise standardization for both soil (difference 0.0456) and combined criteria group (difference 0.0964). Climate criteria resulted in the lowest sensitivity to the selected standardization method, having closely distributed values for all applied standardization methods. The third-order polynomial regression function resulted in the highest R2 in the validation of all soil and combined criteria. The exponential function most accurately represented the relationship of suitability and peak NDVI values for climate criteria with fuzzy and linear standardization methods.
Shares of suitability classes calculated using k-means unsupervised classification regarding the covered area are shown in Table 8. According to the most accurate suitability values calculated using fuzzy standardization and combined criteria, 64.3% of the study area is suitable for soybean cultivation to some degree. The S1 class covers 602 km2 of the county area. The area determined as very suitable covers three major parts: the eastern part near Dalj, mainly characterized with high average minimal air temperatures; the northeastern part by the Baranja hill with optimal C/N values, pH and SoilType values; and the central-south part by the city of Đakovo with suitable SoilTexture and SoilType values. The final suitability maps after the removal of constrained areas from the unsupervised classification results are shown in Figure 7. The most accurate suitability model was exported to vector Shapefile and raster GeoTIFF format.
Suitability classes: very suitable (S1), moderately suitable (S2), marginally suitable (S3), currently unsuitable (N1) and permanently unsuitable (N2).
The correlation matrix of all combinations of calculated suitability models (Table 9) showed a constantly high correlation of suitability values calculated using soil and combined criteria for all three standardization methods. Suitability values using climate criteria produced a moderate correlation with combined criteria values. A high correlation was achieved for fuzzy and stepwise standardization with combined criteria, which produced two most accurate suitability values.

4. Discussion

Linear standardization was proven to be less effective in adjusting to standardization specifications for nearly all quantitative criteria. This occurred due to the linear standardization’s property of monotonous ascending or descending of standardized values in the specified number interval. Therefore, nonlinear varying standardization requirements of Precipitation, SolarIrradiation, Tavg, Tmax, AirHumidity, pH and TWI could not be fully adjusted. However, linear standardization performed as the second-best standardization method for climate criteria, as few extreme values were present in these criteria. Stepwise standardization allowed high subjective impact on standardized results due to the ability to join the arbitrary input value ranges to standardization classes. These classes were discretely defined, with no intermediate values between classes. The most important property of the stepwise standardization method was the ability of the standardization of qualitative criteria, which presents a necessity for the integration with quantitative criteria. Fuzzy standardization methods allowed the best adjustment of standardized values to the crop expert’s specifications for quantitative criteria, consequentially resulting in the highest accuracies of suitability models. Aside from setting the intervals of least and most favorable input values, the crop expert also had supervision over the intensity of increasing and decreasing the standardized values by selecting the optimal fuzzy model. Sigmoidal and J-shaped membership functions were also recognized as the most convenient for the standardization of similar criteria in a study by Aydi et al. [89]. Investigated properties of fuzzy standardization proved that it could be used as a universal standardization method, possibly upgrading the calculation of suitability values where linear and stepwise standardization were previously combined [74]. The FAO guidelines offered a standard for the delineation of suitability classes, which enables a comparison between similar suitability models. These standards were applied in many studies regarding suitability analyses recently [10,15,78]. The creation of suitability classes is considered as a necessary procedure in effective land management in similar studies. They were previously used in suitability analyses regarding crop cultivation and irrigation [99].
Suitability models calculated with the combined criteria group produced the highest accuracy, showing the importance of climate and soil criteria groups in crop-related suitability analysis. A low correlation between models using only climate and soil criteria was observed for all standardization methods. Climate and soil criteria groups were both necessary for reliable suitability calculation [7,102]. The possible reason for the lower accuracy of suitability models that used climate criteria compared to soil or combined criteria is a relatively small study area (4,155 km2), causing low variability for selected climate criteria. Precipitation produced the highest variability of all data from weather stations, with a CV of only 0.100. However, climate data was considered as mandatory criteria group especially in the case of the larger study area, having an impact on high suitability values for combined criteria. Many recent studies successfully integrated GIS-based multicriteria analysis with AHP as a weight determination method for crop suitability calculation [11,25,80,83,84]. Besides weight determination, three modifications of AHP in this study enabled the determination of the reliability of input data used for the modeling of criteria per criteria groups. The combined criteria group produced the highest accuracy for all standardization methods. Fuzzy and stepwise standardization methods also produced moderate accuracy using soil criteria, while climate criteria produced the lowest accuracy. Linear standardization produced similar accuracies for all AHP modifications. Fuzzy standardization resulted in higher accuracy using the combined criteria group than the other two methods, due to better adjustment of values to the crop expert’s specifications in the standardization process. All interpolation methods performed better in case of less variable climate data compared to soil data. OK usually outperformed IDW in studies comparing interpolation methods for soil and climate data but has some limitations regarding normality and stationarity of input samples [103]. IDW outperformed OK in case of low normality and high variability of sample data [104]. Kriging was unable to produce quality results in when its prerequisites of sample data normality and stationarity were not met in a different study, where IDW resulted in higher accuracy [25]. This emphasizes the importance of the selection of optimal interpolation method in similar studies as it depends on the properties of input samples [105].
No standard methodology of crop suitability values validation calculated by GIS-based multicriteria analysis was noted during the literature review. Moreover, the majority of these studies did not employ any validation procedure for crop suitability results [11,25,51,83,84,89]. Multitemporal NDVI offers a potential solution in the validation of crop suitability, due to its global accessibility from global multispectral satellite missions (Sentinel-2, Landsat 8) and being a reliable predictor of various crop properties [92]. Moreover, Sentinel-2 derived NDVI was a reliable predictor in a study of crop yield and was superior to other individual tested vegetation indices [106]. The applicability of NDVI was expanded in this study as it enabled an effective validation of suitability models, primarily due to the objective assessment of all major soybean varieties using peak NDVI values. An application of NDVI in the validation of wheat suitability combined with yield data was successfully made in a study by Dedeoğlu and Dengiz [107]. However, the application of peak NDVI values that were used for direct validation of suitability values before their classification in suitability classes represents a significant improvement for its applicability in suitability validation. Data about conducted agrotechnical operations (sowing, fertilization, spraying) for every test soybean parcel remain an ambiguity and could affect the reliability of validation using peak NDVI. Fertilization is traditionally balanced in the entire study area, but its exact effect on NDVI values for validation remain to be tested. Soybean satisfies around 60% of its nitrogen requirements by nitrogen fixation [108], which is one of the reasons for farmers in the study area to apply conventional fertilization without conducting detailed soil analysis commonly. Therefore, nitrogen fertilization was not expected to have a significant impact during validation using NDVI. The effect of pests, diseases and weather trend on crop status prevent NDVI from being fully reliable validation data, as vegetation depends on numerous individual factors [109]. A relatively small study area partially reduced the effects of weather conditions, having low variability. The averaging of NDVI values from soybean parcels to centroids allowed the lower impact of potentially affected areas by pests or disease. A possible solution for the mentioned restrictions is the integration of multiple complementary vegetation indices for validation, as the spectral resolution of Sentinel-2 allows for the calculation of most known vegetation indices [110]. The Soil-Adjusted Vegetation Index (SAVI) is considered as one of the complementary indices to NDVI, as it minimized the soil brightness effect in areas with sparser vegetation in a comparable study [111]. Of the four Sentinel-2 images used for the calculation of peak NDVI, the dataset sensed on 8 August 2019 produced the lowest amount of peak NDVI values due to the lack of precipitation at the beginning of the month. This case implies the importance of considering meteorological conditions for adjustment of the time frame for validation. The validation time frames should be slightly calibrated for each growing season, as drought could accelerate the ripening process of soybean, therefore lowering the vegetation indices values.
While the simplicity and effectiveness of the weighted linear combination for criteria aggregation were previously mentioned and successfully applied, the Ordered Weighted Averaging (OWA) method offers a potential upgrade to the applied methods. Its potential in GIS-based multicriteria analysis is primarily notable through the integration of fuzzy measures in the criteria aggregation process [112]. OWA was successfully applied for suitability analysis of crops in combination with AHP [113] and multicriteria analysis regarding soil chemical properties in combination with IDW as the interpolation method [114]. The majority of applied methods for GIS-based multicriteria analysis in this study support the automatization of the calculation process, which will be a subject of future research. Inputs consisting of criteria selection, pairwise matrix comparison and standardization specifications should be assigned at the beginning of the algorithm. The main benefits of automatization in this research were noticeable during the process of soil texture map calculation, which are time-efficiency and easier applicability of the process to future calculations [78]. Automatization is considered as a major property for making the GIS-based methods and algorithms considerably more applicable in regular practice [115]. Since the user does not have to conduct each process in the GIS environment manually, the accessibility of automated methods expands to users with limited knowledge of geospatial data.
The addition of the ecological criteria group to the multicriteria analysis presents a possibility for the improvement of the present study. Aside from the ecological criteria group, the proposed method supports the addition of multiple other criteria groups through the creation of additional pairwise comparison matrices, as long as the number of criteria per group satisfies specifications by Saaty and Ozdemir [53].

5. Conclusions

Soybean suitability values were calculated based on agronomist specifications in this study, resulting in maps of soybean suitability classes, as a base for land policy managers to improve soybean quality and reduce production costs. This was achieved by the determination of soybean land suitability using GIS-based multicriteria analysis and AHP. During the standardization process, it was observed that the fuzzy standardization method allowed the best adjustment of standardized values to soybean land suitability specifications. The separation of criteria for soybean suitability calculation to climate and soil groups allowed weight determination according to specifications of AHP. This approach also reduced the redundancy of criteria in the process of pairwise comparison matrices creation. The authors recommend the evaluation and selection of optimal interpolation methods for the preprocessing of input criteria when applicable, as was the case with the most influential criteria in this study. The validation method using Sentinel-2 derived peak NDVI values enabled the determination of the most accurate suitability model. The calculation of peak NDVI values based on four independent NDVI values allowed the objective determination of peak biomass for all major soybean varieties. This procedure was also proven to be robust regarding the cloud cover on satellite images since validation does not depend on the individual images. Fuzzy standardization resulted in the most accurate standardization method combined with all three criteria groups. The suitability model using fuzzy standardization with a combined criteria group produced the highest R2 and RMSE from nine calculated suitability values. According to suitability classes created by unsupervised classification of the most accurate suitability values, 64.3% of the study area was determined as suitable for soybean cultivation in some degree. The most suitable class for soybean cultivation covered 14.5% of the study area, an equivalent of 602 km2 in the field. Three regions of highest suitability for soybean cultivation were observed in the suitability map: the eastern part near Dalj, the northeast part by the Baranja hill and the central-south part by the city of Đakovo. The export of final suitability maps to common raster and vector formats were performed for effective knowledge sharing to agricultural land policy managers.

Author Contributions

Conceptualization: D.R. and M.G.; methodology: D.R.; validation: D.R., M.J. and M.G.; formal analysis: M.J. and I.P.; investigation: D.R.; resources: D.R. and M.G.; data curation: D.R.; writing—original draft preparation: D.R., M.G.; writing—review and editing: D.R., M.J., M.G. and I.P.; visualization: D.R.; supervision: M.J. and M.G.; funding acquisition: M.J. and I.P. All authors have read and agreed to the published version of the manuscript.

Funding

This research received no external funding.

Acknowledgments

This work was supported by the Faculty of Agrobiotechnical Sciences Osijek as a part of the scientific project ‘AgroGIT—technical and technological crop production systems, GIS and environment protection’.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Bajželj, B.; Richards, K.S.; Allwood, J.M.; Smith, P.; Dennis, J.S.; Curmi, E.; Gilligan, C.A. Importance of food-demand management for climate mitigation. Nat. Clim. Chang. 2014, 4, 924. [Google Scholar] [CrossRef] [Green Version]
  2. Pulighe, G.; Bonati, G.; Fabiani, S.; Barsali, T.; Lupia, F.; Vanino, S.; Nino, P.; Arca, P.; Roggero, P.P. Assessment of the Agronomic Feasibility of Bioenergy Crop Cultivation on Marginal and Polluted Land: A GIS-Based Suitability Study from the Sulcis Area, Italy. Energies 2016, 9, 895. [Google Scholar] [CrossRef]
  3. von Arb, C.; Bünemann, E.K.; Schmalz, H.; Portmann, M.; Adamtey, N.; Musyoka, M.W.; Frossard, E.; Fliessbach, A. Soil quality and phosphorus status after nine years of organic and conventional farming at two input levels in the Central Highlands of Kenya. Geoderma 2020, 362, 114112. [Google Scholar] [CrossRef]
  4. Davis, K.F.; Rulli, M.C.; Seveso, A.; D’Odorico, P. Increased food production and reduced water use through optimized crop distribution. Nat. Geosci. 2017, 10, 919. [Google Scholar] [CrossRef]
  5. Purnamasari, R.A.; Noguchi, R.; Ahamed, T. Land suitability assessments for yield prediction of cassava using geospatial fuzzy expert systems and remote sensing. Comput. Electron. Agric. 2019, 166, 105018. [Google Scholar] [CrossRef]
  6. Chen, J.; Yu, J.; Zhang, Y. Multivariate video analysis and Gaussian process regression model based soft sensor for online estimation and prediction of nickel pellet size distributions. Comput. Chem. Eng. 2014, 64, 13–23. [Google Scholar] [CrossRef]
  7. FAO Soils bulletin—A framework for land evaluation. Available online: http://www.fao.org/3/X5310E/x5310e00.htm (accessed on 17 April 2020).
  8. You, L.; Wood, S.; Wood-Sichra, U.; Wu, W. Generating global crop distribution maps: From census to grid. Agric. Syst. 2014, 127, 53–60. [Google Scholar] [CrossRef] [Green Version]
  9. Wang, R.; Jiang, Y.; Su, P.; Wang, J. Global Spatial Distributions of and Trends in Rice Exposure to High Temperature. Sustainability 2019, 11, 6271. [Google Scholar] [CrossRef] [Green Version]
  10. Layomi Jayasinghe, S.; Kumar, L.; Sandamali, J. Assessment of Potential Land Suitability for Tea (Camellia sinensis (L.) O. Kuntze) in Sri Lanka Using a GIS-Based Multi-Criteria Approach. Agriculture 2019, 9, 148. [Google Scholar] [CrossRef] [Green Version]
  11. Plaščak, I.; Jurišić, M.; Radočaj, D.; Barač, Ž.; Glavaš, J. Hazel plantation planning using GIS and multicriteria decision analysis. Poljoprivreda 2019, 25, 79–85. [Google Scholar] [CrossRef]
  12. Yun, H.J.; Kang, D.J.; Kim, D.-K.; Kang, Y. A GIS-Assisted Assessment and Attribute-Based Clustering of Forest Wetland Utility in South Korea. Sustainability 2019, 11, 4632. [Google Scholar] [CrossRef] [Green Version]
  13. Herzberg, R.; Pham, T.G.; Kappas, M.; Wyss, D.; Tran, C.T.M. Multi-Criteria Decision Analysis for the Land Evaluation of Potential Agricultural Land Use Types in a Hilly Area of Central Vietnam. Land 2019, 8, 90. [Google Scholar] [CrossRef] [Green Version]
  14. Gašparović, I.; Gašparović, M. Determining Optimal Solar Power Plant Locations Based on Remote Sensing and GIS Methods: A Case Study from Croatia. Remote Sens. 2019, 11, 1481. [Google Scholar] [CrossRef] [Green Version]
  15. Chen, Y.; Yu, J.; Khan, S. Spatial sensitivity analysis of multi-criteria weights in GIS-based land suitability evaluation. Environ. Modell. Softw. 2010, 25, 1582–1591. [Google Scholar] [CrossRef]
  16. Chakhar, S.; Mousseau, V. Spatial multicriteria decision making. In Encyclopedia of Geographic Information Science; Shehkar, S., Xiong, H., Eds.; Springer: New York, NY, USA, 2008; pp. 747–753. [Google Scholar] [CrossRef]
  17. Ishizaka, A.; Labib, A. Analytic hierarchy process and expert choice: Benefits and limitations. Or Insight 2009, 22, 201–220. [Google Scholar] [CrossRef] [Green Version]
  18. Yu, R.; Tzeng, G.H. A soft computing method for multi-criteria decision making with dependence and feedback. Appl. Math. Comput. 2006, 180, 63–75. [Google Scholar] [CrossRef]
  19. Aragonés-Beltrán, P.; Chaparro-González, F.; Pastor-Ferrando, J.P.; Pla-Rubio, A. An AHP (Analytic Hierarchy Process)/ANP (Analytic Network Process)-based multi-criteria decision approach for the selection of solar-thermal power plant investment projects. Energy 2014, 66, 222–238. [Google Scholar] [CrossRef]
  20. Eastman, J.R. Multi-criteria evaluation and GIS. In Geographical Informational Systems; Goodchild, M.F., Maguire, D.J., Rhind, D.W., Eds.; John Wiley and Sons: New York, NY, USA, 1999; pp. 493–502. [Google Scholar]
  21. Malczewski, J.; Chapman, T.; Flegel, C.; Walters, D.; Shrubsole, D.; Healy, M.A. GIS—Multicriteria evaluation with ordered weighted averaging (OWA): Case study of developing watershed management strategies. Environ. Plan. A 2003, 35, 1769–1784. [Google Scholar] [CrossRef]
  22. Tang, Z.; Yi, S.; Wang, C.; Xiao, Y. Incorporating probabilistic approach into local multi-criteria decision analysis for flood susceptibility assessment. Stoch. Environ. Res. Risk Assess. 2017, 32, 701–714. [Google Scholar] [CrossRef]
  23. Nguyen, T.T.; Verdoodt, A.; Van, Y.T.; Delbecque, N.; Tran, T.C.; Van Ranst, E. Design of a GIS and multi-criteria based land evaluation procedure for sustainable land-use planning at the regional level. Agric. Ecosyst. Environ. 2015, 200, 1–11. [Google Scholar] [CrossRef]
  24. Sallwey, J.; Bonilla Valverde, J.P.; Vásquez López, F.; Junghanns, R.; Stefan, C. Suitability maps for managed aquifer recharge: A review of multi-criteria decision analysis studies. Environ. Rev. 2019, 27, 138–150. [Google Scholar] [CrossRef]
  25. Jurišić, M.; Plaščak, I.; Antonić, O.; Radočaj, D. Suitability Calculation for Red Spicy Pepper Cultivation (Capsicum annum L.) Using Hybrid GIS-Based Multicriteria Analysis. Agronomy 2020, 10, 3. [Google Scholar] [CrossRef] [Green Version]
  26. Koch, B.; Khosla, R.; Frasier, W.M.; Westfall, D.G.; Inman, D. Economic feasibility of variable-rate nitrogen application utilizing site-specific suitability classes. Agron. J. 2004, 96, 1572–1580. [Google Scholar] [CrossRef] [Green Version]
  27. Van Niekerk, A. A comparison of land unit delineation techniques for land evaluation in the Western Cape, South Africa. Land Use Policy 2010, 27, 937–945. [Google Scholar] [CrossRef] [Green Version]
  28. Pagano, M.C.; Miransari, M. The importance of soybean production worldwide. In Abiotic and Biotic Stresses in Soybean Production; Miransari, M., Ed.; Academic Print: Cambridge, MA, USA, 2016; pp. 1–26. [Google Scholar]
  29. The Role of Soybean in Fighting World Hunger. Available online: http://www.fao.org/3/a-bs958e.pdf (accessed on 21 December 2019).
  30. EU Agricultural Outlook—European Commission. Available online: https://ec.europa.eu/info/sites/info/files/food-farming-fisheries/farming/documents/agricultural-outlook-2019-report_en.pdf (accessed on 17 April 2020).
  31. USDA Agricultural Projections to 2029. Available online: https://www.ers.usda.gov/webdocs/publications/95912/oce-2020-1.pdf?v=8056.6 (accessed on 17 April 2020).
  32. Masuda, T.; Goldsmith, P.D. World soybean production: Area harvested, yield, and long-term projections. Int. Food Agribus. Man. 2009, 12, 1–20. [Google Scholar]
  33. Maas, S.J.; Rajan, N. Estimating ground cover of field crops using medium-resolution multispectral satellite imagery. Agron. J. 2008, 100, 320–327. [Google Scholar] [CrossRef]
  34. Immitzer, M.; Vuolo, F.; Atzberger, C. First Experience with Sentinel-2 Data for Crop and Tree Species Classifications in Central Europe. Remote Sens. 2016, 8, 166. [Google Scholar] [CrossRef]
  35. Sentinel-2A Launch. Available online: http://marine.copernicus.eu/wp-content/uploads/2016/06/r2495_9_sentinel_2a.pdf (accessed on 21 December 2019).
  36. Sentinel-2 User Handbook. Available online: https://sentinel.esa.int/documents/247904/685211/Sentinel-2_User_Handbook (accessed on 13 April 2020).
  37. Dong, T.; Meng, J.; Shang, J.; Liu, J.; Wu, B. Evaluation of chlorophyll-related vegetation indices using simulated Sentinel-2 data for estimation of crop fraction of absorbed photosynthetically active radiation. IEEE J. Sel. Top. Appl. Earth Obs. Remote Sens. 2015, 8, 4049–4059. [Google Scholar] [CrossRef]
  38. Veloso, A.; Mermoz, S.; Bouvet, A.; Le Toan, T.; Planells, M.; Dejoux, J.F.; Ceschia, E. Understanding the temporal behavior of crops using Sentinel-1 and Sentinel-2-like data for agricultural applications. Remote Sens. Environ. 2017, 199, 415–426. [Google Scholar] [CrossRef]
  39. Walz, Y.; Wegmann, M.; Dech, S.; Vounatsou, P.; Poda, J.N.; N’Goran, E.K.; Utzinger, J.; Raso, G. Modeling and validation of environmental suitability for schistosomiasis transmission using remote sensing. PLoS Negl. Trop. Dis. 2015, 9, e0004217. [Google Scholar] [CrossRef] [Green Version]
  40. Vanino, S.; Nino, P.; De Michele, C.; Bolognesi, S.F.; D’Urso, G.; Di Bene, C.; Pennelli, B.; Vuolo, F.; Farina, R.; Napoli, R. Capability of Sentinel-2 data for estimating maximum evapotranspiration and irrigation requirements for tomato crop in Central Italy. Remote Sens. Environ. 2018, 215, 452–470. [Google Scholar] [CrossRef]
  41. Conrad, O.; Bechtel, B.; Bock, M.; Dietrich, H.; Fischer, E.; Gerlitz, L.; Wehberg, J.; Wichmann, V.; Böhner, J. System for automated geoscientific analyses (SAGA) v. 2.1.4. Geosci. Model. Dev. Discuss. 2015, 8, 1991–2015. [Google Scholar] [CrossRef] [Green Version]
  42. QGIS Development Team 2019. QGIS Geographic Information System. Open Source Geospatial Foundation Project. Available online: http://qgis.osgeo.org (accessed on 17 April 2020).
  43. GRASS Development Team 2019. Geographic Resources Analysis Support System (GRASS) Software, Version 7.8. Open Source Geospatial Foundation. Available online: http://grass.osgeo.org (accessed on 17 April 2020).
  44. Croatian Bureau of Statistics—Structure of Agricultural Holdings. Available online: https://www.dzs.hr/Hrv_Eng/Pokazatelji/Poljoprivreda%20-%20pregled%20po%20zupanijama.XLSX (accessed on 20 December 2019).
  45. Statistical Yearbook of the Republic of Croatia 2018. Available online: https://www.dzs.hr/Hrv_Eng/ljetopis/2018/sljh2018.pdf (accessed on 21 December 2019).
  46. GalićSubašić, D. Influence of Irrigation, Nitrogen Fertilization and Genotype on the Yield and Quality of Soybean (Glycine max (L.) Merr.). Ph.D. Thesis, Josip Juraj Strossmayer University of Osijek, Osijek, Croatia, 12 December 2019. [Google Scholar]
  47. Agrotechnics of Soybean Cultivation (Glycine max L.). Available online: https://www.savjetodavna.hr/wp-content/uploads/publikacije/AgrotehnikaSojeWeb102018.pdf (accessed on 17 April 2020).
  48. Gibson, K.E.; Yang, H.S.; Franz, T.; Eisenhauer, D.; Gates, J.B.; Nasta, P.; Farmaha, B.; Grassini, P. Assessing explanatory factors for variation in on-farm irrigation in US maize-soybean systems. Agric. Water Manag. 2018, 197, 34–40. [Google Scholar] [CrossRef]
  49. Bagherzadeh, A.; Ghadiri, E.; Darban, A.R.S.; Gholizadeh, A. Land suitability modeling by parametric-based neural networks and fuzzy methods for soybean production in a semi-arid region. Model. Earth Syst. Environ. 2016, 2, 104. [Google Scholar] [CrossRef] [Green Version]
  50. Munene, P.; Chabala, L.; Mweetwa, M. Land Suitability Assessment for Soybean (Glycine max (L.) Merr.) Production in Kabwe District, Central Zambia. J. Agric. Sci. 2017, 9, 1–16. [Google Scholar] [CrossRef] [Green Version]
  51. He, Y.; Liu, D.; Yao, Y.; Huang, Q.; Li, J.; Chen, Y.; Shi, S.; Wan, L.; Yu, S.; Wang, D. Spatializing growth suitability for spring soybean cultivation in northeast China. J. Appl. Meteorol. Climatol. 2013, 52, 773–783. [Google Scholar] [CrossRef]
  52. Seyedmohammadi, J.; Sarmadian, F.; Jafarzadeh, A.A.; Ghorbani, M.A.; Shahbazi, F. Application of SAW, TOPSIS and fuzzy TOPSIS models in cultivation priority planning for maize, rapeseed and soybean crops. Geoderma 2018, 310, 178–190. [Google Scholar] [CrossRef]
  53. Saaty, T.L.; Ozdemir, M.S. Why the magic number seven plus or minus two. Math. Comput. Model. 2003, 38, 233–244. [Google Scholar] [CrossRef]
  54. Huld, T. PVMAPS: Software tools and data for the estimation of solar radiation and photovoltaic module performance over large geographical areas. Sol. Energy 2017, 142, 171–181. [Google Scholar] [CrossRef]
  55. Gašparović, I.; Gašparović, M.; Medak, D. Determining and analysing solar irradiation based on freely available data: A case study from Croatia. Environ. Dev. 2018, 26, 55–67. [Google Scholar] [CrossRef]
  56. Kasten, F. The Linke turbidity factor based on improved values of the integral Rayleigh optical thickness. Sol. Energy 1996, 56, 239–244. [Google Scholar] [CrossRef]
  57. Meteosat Solar Surface Radiation and Effective Cloud Albedo Climate Data Record SARAH-2, Validation Report. Available online: https://www.cmsaf.eu/SharedDocs/Literatur/document/2016/saf_cm_dwd_val_meteosat_hel_2_1_pdf.pdf (accessed on 22 March 2020).
  58. Jantalia, C.P.; Resck, D.V.; Alves, B.J.; Zotarelli, L.; Urquiaga, S.; Boddey, R.M. Tillage effect on C stocks of a clayey Oxisol under a soybean-based crop rotation in the Brazilian Cerrado region. Soil Till Res. 2007, 95, 97–109. [Google Scholar] [CrossRef]
  59. Mesic, M.; Birkás, M.; Zgorelec, Z.; Kisic, I.; Sestak, I.; Jurisic, A.; Husnjak, S. Soil Carbon Variability in Some Hungarian and Croatian Soils. In Soil Carbon; Hartemink, A., McSweeney, K., Eds.; Springer International Publishing: Cham, Switzerland, 2014; pp. 419–426. [Google Scholar] [CrossRef]
  60. United States Department of Agriculture—Soil Survey Manual Handbook No. 18. Available online: https://www.iec.cat/mapasols/DocuInteres/PDF/Llibre50.pdf (accessed on 21 December 2019).
  61. North Dakota State University—Soybean production. Available online: https://library.ndsu.edu/ir/bitstream/handle/10365/5450/a250.pdf?sequence=1 (accessed on 22 March 2020).
  62. Python Software Foundation. Python Language Reference, Version 3.7. Available online: http://www.python.org (accessed on 17 April 2020).
  63. Böhner, J.; Selige, T. Spatial prediction of soil attributes using terrain analysis and climate regionalisation. In SAGA—Analysis and Modelling Applications; Böhner, J., McCloy, K.R., Strobl, J., Eds.; Verlag Erich Goltze: Goltze, Germany, 2006; pp. 13–27. [Google Scholar]
  64. Oliver, M.A.; Webster, R. Geostatistical prediction: Kriging. In Basic Steps in Geostatistics: The Variogram and Kriging; Oliver, M.A., Webster, R., Eds.; Springer International Publishing: Berlin, Germany, 2015; pp. 43–69. [Google Scholar] [CrossRef]
  65. Lu, G.Y.; Wong, D.W. An adaptive inverse-distance weighting spatial interpolation technique. Comput. Geosci. 2008, 34, 1044–1055. [Google Scholar] [CrossRef]
  66. Hofstra, N.; New, M. Spatial variability in correlation decay distance and influence on angular-distance weighting interpolation of daily precipitation over Europe. Int. J. Climatol. 2009, 29, 1872–1880. [Google Scholar] [CrossRef]
  67. Yao, X.; Yu, K.; Deng, Y.; Zeng, Q.; Lai, Z.; Liu, J. Spatial distribution of soil organic carbon stocks in Masson pine (Pinus massoniana) forests in subtropical China. Catena 2019, 178, 189–198. [Google Scholar] [CrossRef]
  68. Qiao, P.; Lei, M.; Yang, S.; Yang, J.; Guo, G.; Zhou, X. Comparing ordinary kriging and inverse distance weighting for soil as pollution in Beijing. Environ. Sci. Pollut. Res. 2018, 25, 15597–15608. [Google Scholar] [CrossRef]
  69. Xavier, A.C.; King, C.W.; Scanlon, B.R. Daily gridded meteorological variables in Brazil (1980–2013). Int. J. Climatol. 2016, 36, 2644–2659. [Google Scholar] [CrossRef] [Green Version]
  70. Cruz-Cárdenas, G.; López-Mata, L.; Ortiz-Solorio, C.A.; Villaseñor, J.L.; Ortiz, E.; Silva, J.T.; Estrada-Godoy, F. Interpolation of Mexican soil properties at a scale of 1: 1,000,000. Geoderma 2014, 213, 29–35. [Google Scholar] [CrossRef]
  71. Malczewski, J. GIS and Multicriteria Decision Analysis; John Wiley & Sons: New York, NY, USA, 1999; p. 408. [Google Scholar]
  72. Kaykhosravi, S.; Abogadil, K.; Khan, U.T.; Jadidi, M.A. The Low-Impact Development Demand Index: A New Approach to Identifying Locations for LID. Water 2019, 11, 2341. [Google Scholar] [CrossRef] [Green Version]
  73. Alamgir, M.; Mohsenipour, M.; Homsi, R.; Wang, X.; Shahid, S.; Shiru, M.S.; Alias, N.E.; Yuzir, A. Parametric Assessment of Seasonal Drought Risk to Crop Production in Bangladesh. Sustainability 2019, 11, 1442. [Google Scholar] [CrossRef] [Green Version]
  74. Bonilla Valverde, J.P.; Blank, C.; Roidt, M.; Schneider, L.; Stefan, C. Application of a GIS Multi-Criteria Decision Analysis for the Identification of Intrinsic Suitable Sites in Costa Rica for the Application of Managed Aquifer Recharge (MAR) through Spreading Methods. Water 2016, 8, 391. [Google Scholar] [CrossRef]
  75. Kumari, M.; Sakai, K.; Kimura, S.; Yuge, K.; Gunarathna, M. Classification of Groundwater Suitability for Irrigation in the Ulagalla Tank Cascade Landscape by GIS and the Analytic Hierarchy Process. Agronomy 2019, 9, 351. [Google Scholar] [CrossRef] [Green Version]
  76. Bezdan, A.; Blagojevic, B.; Vranesevic, M.; Benka, P.; Savic, R.; Bezdan, J. Defining Spatial Priorities for Irrigation Development Using the Soil Conservation and Water Use Efficiency Criteria. Agronomy 2019, 9, 324. [Google Scholar] [CrossRef] [Green Version]
  77. Novák, V. Fuzzy Sets and Their Applications; Taylor & Francis: London, UK, 1989; p. 190. [Google Scholar]
  78. Domazetović, F.; Šiljeg, A.; Lončar, N.; Marić, I. Development of automated multicriteria GIS analysis of gully erosion susceptibility. Appl. Geogr. 2019, 112, 102083. [Google Scholar] [CrossRef]
  79. Hamzeh, S.; Mokarram, M.; Haratian, A.; Bartholomeus, H.; Ligtenberg, A.; Bregt, A.K. Feature Selection as a Time and Cost-Saving Approach for Land Suitability Classification (Case Study of Shavur Plain, Iran). Agriculture 2016, 6, 52. [Google Scholar] [CrossRef] [Green Version]
  80. Zhang, J.; Su, Y.; Wu, J.; Liang, H. GIS based land suitability assessment for tobacco production using AHP and fuzzy set in Shandong province of China. Comput. Electron. Agric. 2015, 114, 202–211. [Google Scholar] [CrossRef]
  81. Klir, G.J. Fuzzy set theory. In Uncertainty and Information: Foundations of Generalized Information Theory; Klir, G.J., Ed.; Wiley-IEEE Press: Hoboken, NJ, USA, 2006; pp. 260–314. [Google Scholar]
  82. Saaty, T.L. A scaling method for priorities in hierarchical structures. J. Math. Psychol. 1977, 15, 234–281. [Google Scholar] [CrossRef]
  83. Musakwa, W. Identifying land suitable for agricultural land reform using GIS-MCDA in South Africa. Environ. Dev. Sustain. 2018, 20, 2281–2299. [Google Scholar] [CrossRef]
  84. Islam, M.M.; Ahamed, T.; Noguchi, R. Land Suitability and Insurance Premiums: A GIS-based Multicriteria Analysis Approach for Sustainable Rice Production. Sustainability 2018, 10, 1759. [Google Scholar] [CrossRef]
  85. Aldababseh, A.; Temimi, M.; Maghelal, P.; Branch, O.; Wulfmeyer, V. Multi-Criteria Evaluation of Irrigated Agriculture Suitability to Achieve Food Security in an Arid Environment. Sustainability 2018, 10, 803. [Google Scholar] [CrossRef] [Green Version]
  86. Nigussie, G.; Moges, M.A.; Moges, M.M.; Steenhuis, T.S. Assessment of Suitable Land for Surface Irrigation in Ungauged Catchments: Blue Nile Basin, Ethiopia. Water 2019, 11, 1465. [Google Scholar] [CrossRef] [Green Version]
  87. Saaty, T.L. Fundamentals of Decision Making and Priority Theory with the Analytic Hierarchy Process; RWS Publications: Pittsburgh, PA, USA, 2000; p. 478. [Google Scholar]
  88. Chhetri, S.; Kayastha, P. Manifestation of an analytic hierarchy process (AHP) model on fire potential zonation mapping in Kathmandu Metropolitan City, Nepal. ISPRS Int. J. Geo Inf. 2015, 4, 400–417. [Google Scholar] [CrossRef] [Green Version]
  89. Aydi, A.; Abichou, T.; Nasr, I.H.; Louati, M.; Zairi, M. Assessment of land suitability for olive mill wastewater disposal site selection by integrating fuzzy logic, AHP, and WLC in a GIS. Environ. Monit. Assess. 2016, 188, 59. [Google Scholar] [CrossRef] [PubMed]
  90. Shahabi, H.; Keihanfard, S.; Ahmad, B.B.; Amiri, M.J.T. Evaluating Boolean, AHP and WLC methods for the selection of waste landfill sites using GIS and satellite images. Environ. Earth sci. 2014, 71, 4221–4233. [Google Scholar] [CrossRef]
  91. Cui, X.; Dong, Y.; Gi, P.; Wang, H.; Xu, K.; Zhang, Z. Relationship between root vigour, photosynthesis and biomass in soybean cultivars during 87 years of genetic improvement in the northern China. Photosynthetica 2016, 54, 81–86. [Google Scholar] [CrossRef]
  92. Vizzari, M.; Santaga, F.; Benincasa, P. Sentinel 2-Based Nitrogen VRT Fertilization in Wheat: Comparison between Traditional and Simple Precision Practices. Agronomy 2019, 9, 278. [Google Scholar] [CrossRef] [Green Version]
  93. Bishop, K.A.; Betzelberger, A.M.; Long, S.P.; Ainsworth, E.A. Is there potential to adapt soybean (Glycine max Merr.) to future [CO2]? An analysis of the yield response of 18 genotypes in free-air CO2 enrichment. Plant. Cell Environ. 2015, 38, 1765–1774. [Google Scholar] [CrossRef]
  94. Sawilowsky, S.S. New effect size rules of thumb. J. Mod. Appl. Stat. Methods 2009, 8, 26. [Google Scholar] [CrossRef]
  95. The Structure of the FAO Framework Classification. Available online: http://www.fao.org/3/x5648e/x5648e0j.htm (accessed on 22 March 2020).
  96. Lolli, F.; Ishizaka, A.; Gamberini, R. New AHP-based approaches for multi-criteria inventory classification. Int. J. Prod. Econ. 2014, 156, 62–74. [Google Scholar] [CrossRef] [Green Version]
  97. Foody, G.M.; Pal, M.; Rocchini, D.; Garzon-Lopez, C.X.; Bastin, L. The Sensitivity of Mapping Methods to Reference Data Quality: Training Supervised Image Classifications with Imperfect Reference Data. ISPRS Int. J. Geo-Inf. 2016, 5, 199. [Google Scholar] [CrossRef] [Green Version]
  98. Bo, L.I.; Zhang, F.; Zhang, L.W.; Huang, J.F.; Zhi-Feng, J.I.N.; Gupta, D.K. Comprehensive suitability evaluation of tea crops using GIS and a modified land ecological suitability evaluation model. Pedosphere 2012, 22, 122–130. [Google Scholar] [CrossRef]
  99. Neji, H.B.B.; Turki, S.Y. GIS—Based multicriteria decision analysis for the delimitation of an agricultural perimeter irrigated with treated wastewater. Agric. Water Manag. 2015, 162, 78–86. [Google Scholar] [CrossRef]
  100. Gašparović, M.; Zrinjski, M.; Barković, Đ.; Radočaj, D. An automatic method for weed mapping in oat fields based on UAV imagery. Comput. Electron. Agric. 2020, 173, 105385. [Google Scholar] [CrossRef]
  101. CLC2018 Technical Guidelines. Available online: https://forum.eionet.europa.eu/nrc_land_covers/library/copernicus-2014-2020/pan-european-component/corine-land-cover-clc-2018/technical-guidelines/clc2018-technical-guidelines/download/en/2/CLC2018TechnicalGuidelines_final.pdf (accessed on 17 April 2020).
  102. Rhebergen, T.; Fairhurst, T.; Zingore, S.; Fisher, M.; Oberthür, T.; Whitbread, A. Climate, soil and land-use based land suitability evaluation for oil palm production in Ghana. Eur. J. Agron. 2016, 81, 1–14. [Google Scholar] [CrossRef] [Green Version]
  103. Gunnink, J.L.; Burrough, P.A. Interactive spatial analysis of soil attribute patterns using exploratory data analysis (EDA) and GIS. In Spatial Analytical Perspectives on GIS, 1st ed.; Ficher, M., Scholten, H., Unwin, D., Eds.; Taylor & Francis Ltd.: London, UK, 2019; pp. 87–100. [Google Scholar]
  104. Gong, G.; Mattevada, S.; O’Bryant, S.E. Comparison of the accuracy of kriging and IDW interpolations in estimating groundwater arsenic concentrations in Texas. Environ. Res. 2014, 130, 59–69. [Google Scholar] [CrossRef]
  105. Heuvelink, G.B.; Pebesma, E.J. Spatial aggregation and soil process modelling. Geoderma 1999, 89, 47–65. [Google Scholar] [CrossRef]
  106. Vrieling, A.; de Beurs, K.M.; Brown, M.E. Variability of African farming systems from phenological analysis of NDVI time series. Clim. Chang. 2011, 109, 455–477. [Google Scholar] [CrossRef] [Green Version]
  107. Dedeoğlu, M.; Dengiz, O. Generating of land suitability index for wheat with hybrid system approach using AHP and GIS. Comput. Electron. Agric. 2019, 167, 105062. [Google Scholar] [CrossRef]
  108. Schipanski, M.E.; Drinkwater, L.E.; Russelle, M.P. Understanding the variability in soybean nitrogen fixation across agroecosystems. Plant. Soil 2010, 329, 379–397. [Google Scholar] [CrossRef]
  109. Schultz, P.A.; Halpert, M.S. Global correlation of temperature, NDVI and precipitation. Adv. Space Res. 1993, 13, 277–280. [Google Scholar] [CrossRef]
  110. Frampton, W.J.; Dash, J.; Watmough, G.; Milton, E.J. Evaluating the capabilities of Sentinel-2 for quantitative estimation of biophysical variables in vegetation. ISPRS J. Photogramm. Remote Sens. 2013, 82, 83–92. [Google Scholar] [CrossRef] [Green Version]
  111. Ren, H.; Zhou, G.; Zhang, F. Using negative soil adjustment factor in soil-adjusted vegetation index (SAVI) for aboveground living biomass estimation in arid grasslands. Remote Sens. Environ. 2018, 209, 439–445. [Google Scholar] [CrossRef]
  112. Jiang, H.; Eastman, J.R. Application of fuzzy measures in multi-criteria evaluation in GIS. Int. J. Geogr. Inf. Sci. 2000, 14, 173–184. [Google Scholar] [CrossRef]
  113. Zabihi, H.; Alizadeh, M.; Kibet Langat, P.; Karami, M.; Shahabi, H.; Ahmad, A.; Nor Said, M.; Lee, S. GIS Multi-Criteria Analysis by Ordered Weighted Averaging (OWA): Toward an Integrated Citrus Management Strategy. Sustainability 2019, 11, 1009. [Google Scholar] [CrossRef] [Green Version]
  114. Mokarram, M.; Hojati, M. Using ordered weight averaging (OWA) aggregation for multi-criteria soil fertility evaluation by GIS (case study: Southeast Iran). Comput. Electron. Agric. 2017, 132, 1–13. [Google Scholar] [CrossRef]
  115. Gašparović, M.; Zrinjski, M.; Gudelj, M. Automatic cost-effective method for land cover classification (ALCC). Comput. Environ. Urban Syst. 2019, 76, 1–10. [Google Scholar] [CrossRef]
Remotesensing 12 01463 g001 550
Figure 1. The workflow of soybean land suitability calculation and evaluation.
Figure 1. The workflow of soybean land suitability calculation and evaluation.
Remotesensing 12 01463 g001
Remotesensing 12 01463 g002 550
Figure 2. Study area: (a) location of Osijek-Baranja County, (b) location of used weather stations and soil samples, (c) location of centroids of test soybean parcels.
Figure 2. Study area: (a) location of Osijek-Baranja County, (b) location of used weather stations and soil samples, (c) location of centroids of test soybean parcels.
Remotesensing 12 01463 g002
Remotesensing 12 01463 g003 550
Figure 3. Weather data at the sensing time of four used Sentinel-2 images from weather station Osijek provided by the Croatian Meteorological and Hydrological Service (DHMZ).
Figure 3. Weather data at the sensing time of four used Sentinel-2 images from weather station Osijek provided by the Croatian Meteorological and Hydrological Service (DHMZ).
Remotesensing 12 01463 g003
Remotesensing 12 01463 g004 550
Figure 4. Preprocessed criteria rasters.
Figure 4. Preprocessed criteria rasters.
Remotesensing 12 01463 g004
Remotesensing 12 01463 g005 550
Figure 5. Graphs of standardization parameters per criterion for evaluated standardization methods.
Figure 5. Graphs of standardization parameters per criterion for evaluated standardization methods.
Remotesensing 12 01463 g005
Remotesensing 12 01463 g006 550
Figure 6. Regression graphs for best fitting functions per the suitability model.
Figure 6. Regression graphs for best fitting functions per the suitability model.
Remotesensing 12 01463 g006
Remotesensing 12 01463 g007 550
Figure 7. Final suitability maps for soybean cultivation for calculated models.
Figure 7. Final suitability maps for soybean cultivation for calculated models.
Remotesensing 12 01463 g007
Table
Table 1. Selected criteria for soybean suitability calculation.
Table 1. Selected criteria for soybean suitability calculation.
Criteria GroupCriterion NameUnitDescriptionSourceNative Format
ClimateTmin°CMean minimum air temperatureDHMZtabular
Tavg°CMean average air temperatureDHMZtabular
Tmax°CMean maximum air temperatureDHMZtabular
PrecipitationmmTotal precipitation amountDHMZtabular
AirHumidity%Mean relative air humidityDHMZtabular
Solar IrradiationkWh
/m2		Mean daily global horizontal irradiationASTER raster
SoilSoilType/Soil type classesBasic pedologic map of Croatiavector (polygon)
pH/Soil pH valuesCAENvector (point)
SoilTexture/Soil texture classesCAENvector (point)
C/N/Carbon-to-nitrogen soil ratioCAENvector (point)
Slope°Terrain slopeSRTMraster
TWI/Topographic wetness indexSRTMraster
Table
Table 2. Identification of Sentinel 2 images, date and the number of peak Normalized Difference Vegetation Index (NDVI) values per image.
Table 2. Identification of Sentinel 2 images, date and the number of peak Normalized Difference Vegetation Index (NDVI) values per image.
Tile IDSatelliteSensing DateDay of YearPeak NDVI Values
T34TCRS2B19th July 201920038
T33TYLS2A24th July 2019205102
T33TYLS2A3rd August 201921549
T33TYLS2B8th August 201922015
Table
Table 3. Descriptive statistics of climate and soil samples for interpolation.
Table 3. Descriptive statistics of climate and soil samples for interpolation.
Criteria namenmeanCVSKKT
Tmin1513.80.0290.075−1.302
Tavg1518.20.0330.284−1.365
Tmax1522.90.046−0.469−1.048
Precipitation15457.50.1001.2750.818
AirHumidity1571.90.0521.0771.689
pH486.70.1590.0611.400
SoilTexture (Clay)4831.50.3320.5250.176
SoilTexture (Silt)4857.80.2150.6000.563
SoilTexture (Sand)4810.71.1901.9454.348
C/N487.60.2510.7471.537
CV: coefficient of variation, SK: skewness, KT: kurtosis.
Table
Table 4. Accuracy assessment of tested interpolation methods.
Table 4. Accuracy assessment of tested interpolation methods.
Input DataOKIDWADW
R2RMSENRMSER2RMSENRMSER2RMSENRMSE
Tmin0.69923.4790.2520.83710.7780.0560.75980.8170.059
Tavg0.74911.9650.1080.81500.7970.0430.81980.8310.045
Tmax0.60604.0910.1780.80291.1550.0500.72481.0500.045
Precipitation0.678772.4510.1580.816525.9150.0560.780626.2690.057
Air Humidity0.65898.1060.1120.73115.0490.0700.68425.9830.083
pH0.59910.6430.0960.73130.5260.0780.74070.6030.090
Clay0.55126.5700.2080.70225.7510.1820.64676.6040.209
Silt0.593010.2990.1780.66959.2470.1600.648710.4580.180
Sand0.61831.4720.1370.65291.4770.1380.64581.5310.143
C/N0.58721.8730.2460.77260.6540.0860.69370.7440.097
The most accurate statistical values per interpolation result were bolded.
Table
Table 5. Pairwise comparison matrix of climate criteria for the weight determination in AHP.
Table 5. Pairwise comparison matrix of climate criteria for the weight determination in AHP.
TminPrecipitationSolar
Irradiation		TavgTmaxAirHumidityWeight
Tmin1233470.388
Precipitation1/2123360.261
SolarIrradiation1/31/212340.156
Tavg1/31/31/21230.104
Tmax1/41/31/31/2130.082
AirHumidity1/71/61/41/31/310.042
n = 6, CI = 0.050, RI = 1.240, CR = 0.040.
Table
Table 6. Pairwise comparison matrix of soil criteria for the weight determination in Analytic Hierarchy Process (AHP).
Table 6. Pairwise comparison matrix of soil criteria for the weight determination in Analytic Hierarchy Process (AHP).
SoilTypepHSlopeSoilTextureC/NTWIWeight
SoilType1233470.383
pH1/2123460.254
Slope1/31/213450.178
SoilTexture1/31/31/31240.104
C/N1/41/41/41/2130.073
TWI1/71/61/51/41/310.039
n = 6, CI = 0.093, RI = 1.240, CR = 0.075.
Table
Table 7. Accuracy assessment values of calculated suitability models using reference peak NDVI values.
Table 7. Accuracy assessment values of calculated suitability models using reference peak NDVI values.
ValueFuzzy standardizationStepwise standardizationLinear standardization
ClimateSoilClimate + SoilClimateSoilClimate + SoilClimateSoilClimate + Soil
R2lin0.40560.68390.84160.30160.61160.69470.32900.48550.6337
R2log0.41430.66440.82730.31750.60460.69440.32490.48190.6279
R2exp0.41910.66720.84170.28440.57600.67710.34220.47550.6289
R2poly20.41610.69070.84290.33190.61170.69750.32930.48580.6337
R2poly30.41620.69230.84380.33260.62910.70950.33100.49290.6357
RMSE0.18740.14350.09620.20700.18910.19260.21560.20890.1925
d0.27170.00490.01470.29870.37200.59900.27800.34340.4738
The most accurate statistical values per suitability result are bolded. R2lin: R2 from linear regression; R2log: R2 from logarithmic regression; R2exp: R2 from exponential regression; R2poly2: R2 from second-order polynomial regression; R2poly3: R2 from third-order polynomial regression; RMSE: root mean square error, d: Cohen’s d index.
Table
Table 8. Shares of suitability classes for calculated models from unsupervised classification.
Table 8. Shares of suitability classes for calculated models from unsupervised classification.
ModelS1 (%)S2 (%)S3 (%)N1 (%)N2 (%)
Fuzzy standardizationClimate13.125.526.223.212.0
Soil23.614.730.521.79.5
Climate + Soil14.522.227.622.513.2
Stepwise standardizationClimate33.426.93.627.58.6
Soil19.123.015.918.523.5
Climate + Soil9.226.325.020.918.6
Linear standardizationClimate16.823.524.624.910.2
Soil19.220.521.615.423.3
Climate + Soil18.625.523.917.114.9
Table
Table 9. Correlation matrix between calculated suitability models.
Table 9. Correlation matrix between calculated suitability models.
FCFSFCS SCSSSCSLCLSLCS
FC1.0000.2850.5780.5780.2780.3960.7180.2420.520
FS 1.0000.9470.3430.9090.8920.3130.7790.736
FCS 1.0000.4860.8660.8910.5060.7430.799
SC 1.0000.3330.5550.4930.2920.452
SS 1.0000.9690.3260.9010.834
SCS 1.0000.4150.8700.853
LC 1.0000.2940.689
LS 1.0000.895
LCS 1.000
FC: fuzzy standardization-climate criteria, FS: fuzzy standardization-soil criteria, FCS: fuzzy standardization-combined criteria, SC: stepwise standardization-climate criteria, SS: stepwise standardization-soil criteria, SCS: stepwise standardization-combined criteria, LC: linear standardization-climate criteria, LS: linear standardization-soil criteria, LCS: linear standardization-combined criteria.

Share and Cite

MDPI and ACS Style

Radočaj, D.; Jurišić, M.; Gašparović, M.; Plaščak, I. Optimal Soybean (Glycine max L.) Land Suitability Using GIS-Based Multicriteria Analysis and Sentinel-2 Multitemporal Images. Remote Sens. 2020, 12, 1463. https://doi.org/10.3390/rs12091463

AMA Style

Radočaj D, Jurišić M, Gašparović M, Plaščak I. Optimal Soybean (Glycine max L.) Land Suitability Using GIS-Based Multicriteria Analysis and Sentinel-2 Multitemporal Images. Remote Sensing. 2020; 12(9):1463. https://doi.org/10.3390/rs12091463

Chicago/Turabian Style

Radočaj, Dorijan, Mladen Jurišić, Mateo Gašparović, and Ivan Plaščak. 2020. "Optimal Soybean (Glycine max L.) Land Suitability Using GIS-Based Multicriteria Analysis and Sentinel-2 Multitemporal Images" Remote Sensing 12, no. 9: 1463. https://doi.org/10.3390/rs12091463

APA Style

Radočaj, D., Jurišić, M., Gašparović, M., & Plaščak, I. (2020). Optimal Soybean (Glycine max L.) Land Suitability Using GIS-Based Multicriteria Analysis and Sentinel-2 Multitemporal Images. Remote Sensing, 12(9), 1463. https://doi.org/10.3390/rs12091463

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop