## 4 Discussion
[13] examined the effectiveness of machine learning techniques compared to traditional parametric methods for predicting the suitability of wheat and barley, finding that machine learning improved accuracy by 29% for wheat and 26% for barley. Additionally, their research highlighted the relative importance of abiotic input factors, offering a systematic approach to weighting in the AHP. Similarly, [PERSON] et al. [48] underscored the advantages of machine learning in suitability assessments, focusing on the precise identification of individual suitability components while also considering socioeconomic and environmental factors. While this study used the LAI and FAPAR, a previous study successfully enhanced traditional GIS-based multicriteria analysis by applying a novel method for evaluating cropland suitability utilizing the Normalized Difference Vegetation Index (NDVI) derived from Sentinel-2 multispectral imagery [21]. NDVI peak values were measured during the full seed growth stage (R6) of soybean, as was performed in this study, aligning with prior findings that indicated a strong correlation between the NDVI and soybean grain yield at this growth phase [28]. This method presents a promising alternative for validation processes traditionally employed in agricultural assessments and can be adapted for different crop types exhibiting significant correlations between satellite-derived vegetation indices and actual yield metrics during various growth stages. Additionally, the LAI and FAPAR were also utilized successfully in a classification-based approach from the PROBA-V satellite data to establish reference datasets for machine learning models predicting cropland suitability, given their strong associations with biomass and crop yield [49]. There are three main drawbacks to the proposed cropland suitability prediction method, which should be addressed in future studies:  
* The coverage of training and test data from 2015 to 2021 relative to the total agricultural land in continental Croatia resulted in range of 0.17-0.38% (with a minimum threshold of 0.25%). In addition to having lower coverage for 2016, this procedure might not be viable for crop types that are less represented than soybean in continental Croatia. Additionally, a correct procedure is based on 300 m spatial resolution, which requires the presence of agricultural parcels larger than 4.5 ha, which was seldom present and varies on a year-to-year basis due to the necessity of complying to crop rotation systems. These observations strongly suggest that 300 m spatial resolution, while applicable for cropland suitability prediction, is not optimal and viable for less represented crops. Therefore, the training and test data should be modeled in up to 30 m spatial resolution, as this would ensure stable training and test data coverage, without any modification of present cropland suitability prediction method other than modeling environmental criteria in the same spatial resolution. * The minimum coverage of ground truth soybean parcels of 50% per pixel of 300 m spatial resolution from the reference grid potentially caused spectral mixing with forests or other vegetation with high peak LAI and FAPAR values, as higher cropland suitability for soybean cultivation frequently occurred in close proximity to forests. This is an additional reason for introducing an up to 30 m spatial resolution reference grid for training and test data, as well as requiring 100% coverage of ground-truth soybean parcels within a pixel of a reference grid. The selection of optimal spatial resolution could be performed according to external accuracy assessment using independent ground-truth data, such as crop yield data. * The translation from K-means classes to FAO land suitability standards in this study was performed using a relative-based approach, ranking five resulting classes based on their relative mean sums of LAI and FAPAR values.