Soil aggregate stability (SAS) is a critical parameter of soil quality and its mapping can help determine erosion hotspots. Despite this importance, SAS is less documented in available literature due to limited number of analyzes besides being a time consuming. For this reason, many researchers have turned to alternative methods that often use readily available variables such as soil parameters or remote sensing indices to estimate this variable. In that framework, the aim of the present study focused on the investigation of the feasibile use of adapted Leo Breiman’s random forest algorithm (RF) to mapping different mean weight diameter (MWD) tests as an index of SAS (mechanical breakdown (MWDmb), slow wetting (MWDsw), fast wetting (MWDfw) and the mean of the three tests (MWDmean)). The model was built with 77 samples distributed in the three watersheds of the study area located at Settat Ben-Ahmed, in Morocco and with the use of several environmental variables such as soil parameters (organic matter and clay), remote sensing indices (band 2, band 3, band 4, band 5, normalized difference vegetation index (NDVI) and transformed normalized difference vegetation index (TNDVI)), topography (elevation, slope, curvature plane and the topographic wetness index (TWI)) along with additional categorical variables as geological maps, land use and soil classes. The results showed a good level of accuracy for the training phase (75% of samples) for the different tests (R2 > 0.92, RMSE and MAE < 0.15) and were satisfactory for the testing phase (25% of samples, R2 > 0.65, RMSE and MAE < 0.31). Also, organic matter, topography and geology were the most important parameters in the spatial prediction of SAS. Finally, the maps build during this study could be of great use to identify areas of less stable soils in the perspective for taking the necessary measures to improve their quality.