Hybrid deep learning with a random forest system for sustainable agricultural land cover classification using DEM in Najran, Saudi Arabia

3.1 Study area

Najran is a culturally and economically significant city located in southern Saudi Arabia, and it is ideal for its wealthy records, lifestyles, and beautiful landscapes. The town mixes contemporary tendencies and historical records, making it a unique region. Najran, a city in southwestern Saudi Arabia, is approximately 17.49°N latitude and 44.13°E longitude [39]. Known for its rapid development, the city boasts a rich cultural heritage and diverse landscapes, including valleys, mountains, and desert terrain. Najran is well known for its old websites and historical al-Uukardud ruins, reducing the date again by a few years. Metropolis also has traditional Madbrica palaces, which include the luxurious al-Aa-Palas, which show the architectural splendor of the region. The panorama of Najran is a mixture of mountains, valleys, and deserts. The region has a medium climate with vanity weather and warm summer. In Najran, juicy palm and agricultural areas produced luxurious dates and fruits. Overall, Najran is a city that fully combines records, lifestyle, and natural beauty, making it a remarkable and captivating part of Saudi Arabia [40].

3.2 Satellite and DEM data preprocessing

Radiometric correction: The raw digital numbers of the IRS AWiFS imagery were converted to top-of-atmosphere (TOA) reflectance using the radiometric rescaling coefficients provided in the product metadata. The conversion was performed using equation (1):

where ϱ λ is the TOA reflectance for band λ, Gain is the absolute calibration gain coefficient, and Bias is the absolute calibration bias coefficient. These coefficients were sourced from the official IRS data user guide provided by the National Remote Sensing Centre (NRSC).

Atmospheric correction: Following the TOA conversion, applied the dark object subtraction (DOS) method to compensate for atmospheric scattering and obtain surface reflectance. This model was selected for its effectiveness and widespread use with medium-resolution multispectral data like AWiFS. The DOS algorithm was implemented using the atmcorr function within the QGIS semi-automatic classification plugin, where the haze value for each band was determined by identifying the darkest pixel (∼1% reflectance) in a clear water body or deep shadow area within the scene.

Geometric correction: The level-1 IRS AWiFS products utilized are systematically terrain-corrected and georeferenced by NRSC using ground control points and a DTM based on the WGS84 datum and UTM projection.

DEM preprocessing: The DEM, acquired from the USGS EarthExplorer portal, was coregistered and resampled spatial resolution of the IRS AWiFS data. Used the cubic convolution resampling procedure. Selected this method over options that were easier (e.g., nearest neighbor) because it yields a smoother surface and more accurately represents continuous topographic gradients and is necessary for subsequent slope and aspect calculations.

Image acquisition and selection is the first step to acquiring and selecting satellite imagery for the study area. Features of IRS AWiFS Imagery Satellite in Table 1, the imagery should be chosen to cover the study area and period of interest. IRS AWiFS is a wide-angle medium resolution camera with a swath of 740 km (FOV = ±25°) of WiFS heritage. It is carried by the Resourcesat-1 and Resourcesat-2 satellites, launched by the ISRO in 2003 and 2011, respectively. AWiFS observes in four spectral bands: green, red, near infrared (NIR), and short-wave infrared (SWIR).

Once the imagery has been acquired, it must be preprocessed. It includes steps in geometric, radiometric, and atmospheric correction. Geometric correction corrects for any geometric distortions in the imagery caused by the Earth’s curvature and satellite sensor characteristics. Radiometric correction corrects for any radiometric distortions in the imagery caused by sensor noise and atmospheric scattering. Atmospheric correction corrects for the effects of the atmosphere on the imagery, which are absorption and scattering: study location and clipping case study in IRS-AWIFS 2020. The DEM, acquired from the USGS Earth Explorer website, is also needed for this study. The DEM has been prepared for use with the VGG16 and EfficientNetB7. This includes resampling the DEM to the exact spatial resolution of the imagery and converting it to a format compatible with the VGG16 and EfficientNetB7, as shown in Figure 1.

Figure 1

Preprocessing workflow: Resampling and clipping of digital elevation model for Integration with VGG16 and EfficientNetB7 in Najran mapping.

3.3 Training data generation and DEM-based feature extraction

This training data generation process is designed for this research endeavor and is focused on a comprehensive analysis of landscape alterations within Najran. The study employs a fusion of VGG16 and EfficientNetB7 with DEM data to explore the intricacies of urban LULC in Najran. The primary data source is the high-resolution satellite imagery from the IRS AWIFS platform. This work created training data to enable the development and validation of machine-learning models for accurate LULC detection and classification.

Feature extraction and data splitting are the basic steps involved in working with DEMs. Feature extraction transforms elevation data into meaningful variables that characterize the terrain in terms of slope (steepness), aspect (direction of the steepest slope), and curvature (surface convexity or concavity). Other features analyzed were terrain profiles, flow direction, flow accumulation for hydrological modeling, watershed delineation, and texture measures describing spatial variation. These features were used for landform classification and quantification of landscape geometry.

This section provides a rigorous mathematical description of the terrain derivatives extracted from the preprocessed DEM.

The calculations were performed using the Terrain Analysis tools in GIS (System for Automated Geoscientific Analyses [SAGA]), which employs robust algorithms based on established geomorphometric.

The key first and second-order terrain derivatives extracted were as follows:

Slope: Calculated in degrees using Horn’s algorithm (a third-order finite difference method). The algorithm calculates the maximum rate of change in elevation between a cell and its eight neighbors. Equation (2) for the slope angle (β) is as follows:

where δz/δx and δz/δy are the first derivatives of elevation in the x and y (easting and northing) directions, respectively.

Aspect: Also calculated using Horn’s algorithm. Aspect (α) represents the compass direction (0–360° from north) of the steepest downhill slope. It is derived as in equation (3):

The result is then converted from radians to degrees and adjusted so that 0° is true north.

Curvature: calculated two specific types of curvature:

Planform curvature (Plan_Curv): Curvature perpendicular to the direction of the maximum slope. It influences flow convergence and divergence. A negative value indicates a concave surface (convergent flow), a positive value indicates a convex surface (divergent flow), and a value of zero indicates a planar surface. It is calculated as in equation (4):

where p = δz/δx, q = δz/δy, r = δ²z/δx², s = δ²z/δxδy, and t = δ²z/δy².

Profile curvature (Prof_Curv): Curvature parallel to the direction of the maximum slope. It influences flow acceleration and deceleration. A negative value indicates a convex surface (flow acceleration), a positive value indicates a concave surface (flow deceleration). It is calculated as in equation (5):

Prof_Curv = −(p²r + 2pqs + q²t)/(p² + q²) (1 + p² + q²)^(3/2)

Additional derived features included the Topographic Wetness Index and the Terrain Ruggedness Index, the formulas for which have also been explicitly stated in the revised manuscript.

Data splitting divides the extracted features into training, validation, and testing so that the model development remains unbiased and reliable. The training consisted of 60%, with 20% used for validation and 20% used for testing in this application. The fivefold cross-validation is utilized, wherein every fold is served once as the test to increase its generalizability. Stratified sampling is performed so that the classes can remain proportionally represented, whereas random and systematic splits are performed to solve data balance and spatial correlation problems.

4.1 DEM

DEM provides an essential dataset for studying LULC, allowing researchers to quantify and visualize elevation differences over time. This information is precious for diverse applications, such as environmental monitoring, city planning, agriculture, and herbal aid management. DEM was crucial in studying LULC and became pivotal in geology, environmental, technological know-how, city planning, forestry, and agriculture. DEMs are virtual representations of the Earth’s surface, providing data about a place’s elevation, terrain, and topography. These models had been created through remote sensing techniques, LiDAR, or photogrammetry, and they are essential tools for know-how and tracking changes inside the Earth’s landscape. DEM changed into an imperative tool for studying land adjustments. They offer a detailed and correct representation of the Earth’s surface, permitting researchers, land managers, and policymakers to screen and respond to environmental, geological, and anthropogenic changes in the landscape. DEMs facilitate better choice-making and contribute to the sustainable control of herbal resources and the maintenance of the planet [29].

4.2 Deep learning models for feature extraction

CNNs are, in particular, properly perfect for responsibilities in image classification, object detection, and image segmentation because of their ability to extract hierarchical and meaningful features from images. CNN extracts local spatial features from an image through convolutional and pooling layers and then combines these features in deeper layers to combine higher-order representations, which are more appropriate for the task and image classification. This hierarchical feature extraction system enables CNNs to excel at various computer vision tasks, leading to their widespread adoption in artificial intelligence [41]. Figure 2 illustrates the operation of a CNN in extracting and transforming features from images. The CNN systematically captures local spatial details and integrates them into more complex higher-order parts, subsequently employed for distinguishing between different image categories. The simplified depiction offers insights into generating local spatial features and their evolution into higher-level features within a CNN [42]. Multilayer convolution filters generate feature maps in the initial CNN layer following the input image. For each multilayer filter, the CNN undergoes the following steps:

Figure 2

Architectural overview of the proposed convolutional NN for LULC classification in Najran.

Convolution: The filter convolves with the input image. A bias term linked to the particular feature map contributes to each convolutional sum. The outcomes pass through a non-linear activation function. Processed outputs populate a feature map. The amalgamation of the non-linear function, convolution sums, and biases empowers CNN to discern intricate non-linear distinctions among diverse image types. The feature maps created during the initial CNN layer encapsulate localized spatial attributes. The subsequent CNN layer elevates these localized spatial features into more intricate higher-order features using pooling or feature aggregation. Each feature map is independently subjected to averaging or maximum pooling, amalgamating the features. This process recurs through a series of iterations, yielding a succession of compact feature maps. Depending on the CNN’s architecture, the terminal layer of feature maps then interfaces with a fully connected neural network (NN) and, ultimately, a multi-class classification module. This cumulative approach enables the CNN to progressively identify and consolidate intricate image characteristics, contributing to its capacity for accurate classification [13,43].

4.3 VGG16 model

VGG16, a CNN model, is a powerful model for LULC analysis. Its superior ability to extract high-resolution features from satellite images makes it an ideal LULC classification and change detection model [44]. LULC analysis is crucial in environmental science, urban planning, agriculture, and other fields. It involves observing and measuring LULC over time, providing valuable insights for decision-makers and researchers. VGG16 can potentially address this challenge by extracting fine-grained LULC features. LULC analysis involved classifying satellite images into different LULC classes.

The network was converted into a powerful feature extractor by removing the top classification layers. These features captured hierarchical and abstract information about the land’s appearance, crucial for distinguishing different LULC types. VGG16 was fine-tuned for end-to-end LULC feature extraction [45]. By adapting the network’s output layer to match the number of LULC classes, it learns to directly predict LULC types from input images. Fine-tuning allows the model to adapt to the specific characteristics of the LULC dataset, potentially improving classification accuracy. LULC analysis involves detecting changes between different periods. These helped identify land transformation, deforestation, or urban expansion, which represents an invaluable insight, especially in environmentally challenged areas [46].

VGG16 extended to handle time series data. The model learned to recognize temporal patterns in LULC, seasonal variations, or long-term trends by incorporating multiple images from different time steps as input. The accuracy of LULC analysis heavily relied on the quality of input data.

Additionally, the definition of LULC classes varies depending on the specific application. Training and fine-tuning the deep NN VGG16 were computationally intensive. Adequate computational resources were necessary for efficient model training, as in Figure 3.

Figure 3

Enhanced VGG16 architecture for multitemporal feature extraction in LULC analysis of Najran.

4.4 EfficientNetB7 model

EfficientNetB7 is a state-of-the-art deep-gaining knowledge version designed for high-overall performance photo-type tasks. It is a part of the efficient net own family, which uses a compound scaling technique to optimize intensity, width, and determination in a balanced way. This structure improves the computational efficiency and accuracy compared to traditional CNNs. EfficientNetB7 operates on the idea of numerous key improvements [47]. Unlike traditional CNNs, which arbitrarily scale network dimensions, EfficientNetB7 applies a principled scaling technique. It proportionally increases the community intensity (range of layers), width (range of channels), and input resolution to achieve the most fulfilling performance. EfficientNetB7 utilizes intense-ty-smart separable convolutions to reduce the variety of parameters and computational complexity while maintaining feature extraction skills [48]. This method enhances efficiency by independently applying convolutions to every channel after combining the outcomes. SE blocks dynamically modify the significance of characteristic channels by recalibrating them primarily based on the found-out interest mechanisms. This improves the capacity of the model to seize meaningful spatial and contextual statistics. Unlike the conventional ReLU activation characteristic, EfficientNetB7 employs the Swish activation characteristic, which allows small negative values. This results in smoother gradient glide and improved feature getting to know. The model incorporates batch normalization to stabilize the schooling and dropout layers and save you from overfitting, ensuring strong generalization to new datasets. EfficientNetB7 achieves superior accuracy with fewer computational resources than previous deep-mastering models by integrating those improvements. In the context of this study, EfficientNetB7 was fine-tuned on the DEM and multi-temporal IRS AWIFS satellite data, enabling precise LULC classification for Najran. Its optimized feature extraction capabilities significantly enhanced classification accuracy when combined with the VGG16 model and RF classifier.

4.5 GOA

GOA is a nature-inspired optimization algorithm based on grasshoppers’ swarming behavior. It is commonly used for optimization problems, and its application to feature reduction in the CNN models has several potential benefits, especially in analyzing LULC using satellite imagery or maps like the Najran map. CNN models often have many parameters and features, making them computationally expensive. Feature reduction using algorithms like GOA helps select the most relevant features, thus reducing the computational burden. Feature reduction aims to retain essential information while eliminating irrelevant or redundant features [49]. It led to improved generalization of the proposed model, making it more robust and less prone to overfitting. A reduced set of features is often easier to interpret and understand. It is crucial in LULC analysis, as it allows researchers and decision-makers to focus on the most significant factors influencing the changes in the landscape. Reducing the number of features helps filter out noise or irrelevant information from the input data. It is essential when dealing with remote sensing data, where the presence of noise or outside features affects the accuracy of the analysis. Using GOA for feature reduction, aim to find an optimal subset of features that contributes most to the task. It leads to better model performance and more accurate predictions in the LULC analysis [50].

The working mechanism of GOA for feature reduction typically involves an iterative process where a population of potential feature subsets is evaluated based on a fitness function. The grasshoppers in the algorithm represent possible solutions, and their movement is guided by mathematical equations inspired by the swarming behavior of real grasshoppers. The algorithm iteratively refines the feature subsets until convergence to an optimal or near-optimal solution.

4.6 A proposed hybrid VGG16 and EfficientNetB7 with based DEM

DEM is a representation of terrain elevation, typically in a raster format. The hybrid system work has preprocessed the DEM data. This involved scaling, normalizing, and augmenting the data to make it suitable for deep-learning model training. The hybrid system is a multi-CNN with RF designed for image classification tasks. This work has used imagery data corresponding to the geographic area covered by the DEM for Najran. Then, it extracted DEM Features using techniques of convolutional layers or other NN architectures designed to work with raster data to extract meaningful features from the DEM. A hybrid system has been used to extract features from the imagery data and combine the features extracted from the DEM and imagery data. This was done by combining the feature vectors of VGG16 and EfficientNetB7 and then classifying them using an RF classifier.

Figure 4 outlines the main stages – from preprocessing and feature extraction of satellite imagery and DEM data, through feature fusion based on convolutional NNs, to the final classification via an optimized RF classifier.

Figure 4

Integrated hybrid framework combining VGG16, EfficientNetB7, and RF for landform change classification in Najran.

Data acquisition sources (IRS AWiFS imagery and USGS DEM) and preprocessing steps, including normalization and augmentation.

The extraction of terrain derivatives and image features via convolutional layers within VGG16 and EfficientNetB7 architectures.

The fusion of feature vectors from both CNNs prior to classification.

The hyperparameter optimization of the RF classifier (number of trees, max depth, min samples split/leaf), conducted via Bayesian optimization targeting an optimized macro F1-score.

The k-fold cross-validation strategy is employed for robust model evaluation.

These explicitly specify the tools – such as the use of the GO algorithm for feature selection and SAGA GIS for DEM processing – and key parameters utilized throughout the workflow.

LULC refers to the physical entities on the Earth’s surface, encompassing natural and human-made features (Table 2).

Hyperparameter optimization: n_estimators: In this study, the number of trees in the forest was 350. max_depth: In this study, the maximum depth for each tree was 45. min_samples_split: In this study, the minimum number of samples required to split an internal node was 2. min_samples_leaf: In this study, the minimum number of samples required to be present in a leaf node was 1.

The hyperparameters for the RF classifier were systematically optimized using a Bayesian optimization framework, which is more efficient than a grid or random search for high-dimensional parameter spaces. The optimization was conducted to maximize the macro-averaged F1-score on a held-out validation set (20% of the training data).

This optimized configuration of the RF classifier was then used in the final hybrid model (VGG16–EfficientNetB7 with RF) for all subsequent evaluations and comparisons reported in the manuscript (Tables 4, 5). The performance improvement observed in Table 5 for the hybrid model over the standalone VGG16 and EfficientNetB7 models is a direct result of this systematic feature integration and hyperparameter tuning.

This hybrid system has improved the accuracy and efficiency of various geospatial tasks and is useful in monitoring landform changes. Further research and experimentation are needed to realize the potential of this hybrid system in practical applications. The model was trained on geospatial data to perform LULC classification in Najran. Classification of DEMs using the hybrid system of VGG16 and EfficientNetB7, with RF enhanced classification results for geospatial data, led to more effective analysis results.

Deep learning models automatically learned relevant features from satellite imagery, which has reduced the need for manual feature engineering. Combining elevation information with visual features provided a more comprehensive understanding of the Earth’s surface, enabling better decision-making in geospatial applications. The proposed hybrid system combined the strengths of DEMs and VGG16–EfficientNetB7 to enhance geospatial analysis. The integrated features were incorporated into the GO algorithm for spatial geospatial analysis and feature selection, and LULC classification and terrain analysis tasks were performed using an RF classifier.

5.1 Results

Model validation strategy: A rigorous k-fold cross-validation (with k = 5) was employed to ensure the robustness and generalizability of results. The dataset was partitioned into five distinct folds, ensuring each LULC class was proportionally represented in both training and validation sets. The model was trained on four folds and validated on the remaining one; this process was repeated five times, with each fold used exactly once as the validation set. The performance metrics reported in Table 4 represents the mean values obtained across all five validation folds. This method provides a more reliable estimate of model performance than a single train-test split, as it reduces the variance of the estimate and mitigates overfitting.

The area values in Table 3 were derived from the final classified raster output of the proposed hybrid VGG16–EfficientNetB7 + RF model. The methodological steps for generating these values were as follows:

Classification: The study area was classified using the proposed hybrid model.

Raster calculation: The area for each LULC class was calculated directly from the classified raster using the formula:

Pixel area: The spatial resolution of the source satellite images used for classification was 30 m (Landsat 8 OLI). So, the area per pixel is: 0.03 km × 0.03 km = 0.0009 km².

The results of a hybrid system combining VGG16 and Effi-cientNetB7 with an RF algorithm based on DEM data. These systems are evaluated based on various features related to LULC classes. Table 3 presents the LULC classification for Najran, a city, based on different LULC types and their respective areas. The classification provides insights into the spatial distribution of vegetation, water bodies, barren land, built-up areas, agricultural land, and other LULC categories. The results indicate that barren land constitutes the largest proportion of Najran’s total area, covering approximately 80,000 km² (53.5%).

This aligns with the region’s arid and semi-arid characteristics, where desert landscapes dominate the terrain. Agricultural land accounts for 16.7% (25,000 km²), reflecting cultivated areas, primarily for date palm production and other crops. Built-up areas, including urban infrastructure and residential zones, cover 13.4% (20,000 km²), highlighting Najran’s development and urban expansion. Vegetation occupies 10.0% (15,000 km²) of the total land area, indicating the distribution of natural and managed green spaces across the city. Water bodies constitute a minimal 1.0% (1,500 km²), consistent with the arid climate, where permanent water sources are limited. The “Others” category, representing mixed LULC types or unclassified regions, accounts for 5.4% (8,011 km²). The data reflect Najran’s environmental and developmental characteristics, where desert landscapes dominate, but agricultural and urban areas contribute significantly to land use.

5.2 Evaluation of proposed hybrid system

The proposed system generated a confusion matrix of the predicted LULC classes with the reference data to assess the classification accuracy of LULC mapping. The confusion matrix provides insights into classification performance, including accuracy, precision, recall, F1-score, and specificity. The kappa coefficient is also calculated to assess the agreement. The reference dataset comprised six LULC categories: vegetation, water bodies, barren land, built-up areas, agricultural land, and others. The classification results were cross-validated using an independent test dataset to compute the confusion matrix. Figure 5 presents the confusion matrix, where diagonal elements indicate correctly classified instances and off-diagonal elements represent misclassifications.

Figure 5

Heatmap visualization of LULC classification confusion matrix for Najran using a hybrid system.

The hybrid system achieved high classification performance across all LULC classes, as in Table 4 and Figure 6, with overall accuracy exceeding 90% for most categories. The combination of VGG16 and EfficientNetB7 extracts deep, hierarchical features, while RF enhances classification through ensemble learning. Water bodies had an accuracy of 97.78% and a specificity of 99.06%, indicating that the model effectively distinguishes water from other LULC types. Barren land showed a recall of 89.35% and precision of 87.25%, confirming the model’s ability to detect most barren areas while maintaining a low false positive rate. Vegetation and Agri-cultural Land recall of 73.30 and 74.83%, respectively), suggesting some misclassification with other vegetative classes. Built-up Areas and Others had moderate recall values of 77.44 and 70.06%, indicating slight overlaps with similar classes, but high specificity of 96.58 and 97.85% ensured minimal false positives. The hybrid approach significantly improves feature extraction by leveraging both VGG16 and EfficientNetB7. The RF classifier enhances decision boundaries, producing robust classification with high specificity and accuracy. Despite slight misclassifications in complex classes like vegetation and built-up areas, the overall performance indicates strong generalization ability.

Figure 6

Display performance results of the hybrid system for classifying the terrain of the Najran map.

Table 5 shows the overall accuracy and kappa coefficient evaluation results of a hybrid system combining VGG16 and EfficientNetB7 with the RF and VGG16–EfficientNetB7 models based on DEM data. The assessment uses two performance metrics: the kappa coefficient and accuracy. These metrics are commonly used to assess the quality of classification models. The hybrid system, which combines VGG16–EfficientNetB7 with an RF model, achieves a kappa coefficient of 0.9338 and an accuracy of 94.2%. In this case, the kappa coefficient of 0.9338 indicates a high level of agreement between the model’s predictions and the actual values. The overall accuracy of 94.2% indicates that the hybrid system effectively classifies images of Najran and assigns each pixel to its appropriate class. The table shows the results of the hybrid system, which combines VGG16 and EfficientNetB7 with the radio frequency model.

5.3 Sensitivity analysis and discussion of biases

Sensitivity to DEM resolution: This methodology was reprocessed using two additional publicly available DEMs with different spatial resolutions: ALOS PALSAR (12.5 m) and SRTM (1 arc-second, ∼30 m). The model performance (overall accuracy and kappa) remained consistently high (>92% OA) across all resolutions, demonstrating that the proposed hybrid feature extraction approach is robust to variations in the DEM source and scale. This result indicates that topographic feature extracts (slope, aspect, and curvature) are meaningful across resolutions relevant to regional LULC studies.

Sensitivity to seasonal variation (temporal robustness): To evaluate performance across seasons, a satellite image was acquired and preprocessed from a different season (October 2020, representing a post-summer period) for the same study area. The model trained on the primary dataset was applied to this new temporal scene. While overall accuracy remained high (91.5%), observed a predictable decrease in precision for the “vegetation” and “agricultural land” classes due to phenological changes. This analysis confirms the model’s strong generalizability, while also highlighting the importance of accounting for seasonal signatures in operational applications, a point that is now explicitly discussed.

Discussion of dataset bias: have added a dedicated paragraph discussing the class imbalance evident in Table 3, where “barren land” constitutes 53.5% of the study area. Acknowledge that this is an inherent characteristic of the arid environment of Najran, not an artifact of data collection. discuss how this imbalance can inflate overall accuracy metrics and the steps taken to mitigate its effects: namely, the use of stratified k-fold cross-validation (which ensures proportional representation of all classes in each fold) and the reporting of class-specific metrics (accuracy, precision, recall, F1-score) in Table 4, which provide a more truthful account of performance for minority classes like “water bodies.” The high F1-scores across all classes, including minorities, indicate that the model learned discriminative features rather than simply exploiting the majority class.