Disastrous natural hazards, such as landslides, floods, and forest fires cause a serious threat to natural resources, assets and human lives. Consequently, landslide risk assessment has become requisite for managing the resources in future. This study was designed to develop four ensemble metaheuristic machine learning algorithms, such as grey wolf optimized based artificial neural network (GW-ANN), grey wolf optimized based random forest (GW-RF), particle swarm optimization optimized based ANN (PSO-ANN), and PSO optimized based RF for modeling rainfall-induced landslide susceptibility (LS) in Aqabat Al-Sulbat, Asir region, Saudi Arabia, which observes landslide frequently. To obtain very high precision and robust prediction from machine learning algorithms, the grey wolf and PSO optimization algorithms were integrated to develop new ensemble machine learning techniques. Subsequently, LS maps produced by training dataset were validated using the receiver operating characteristics (ROC) curve based on the testing dataset. Based on the area under curve (AUC) value of ROC curve, the best method for LS modeling was selected. We developed ROC curve-based sensitivity analysis to investigate the influence of the parameters for LS modeling. The Gumble extreme value distribution was employed to estimate the rainfall at 2, 5, 10, 20, 50, and 100 year return periods. Then, the landslide hazard maps were prepared at different return periods by integrating the best LS model and estimated rainfall at different return periods. The theory of danger pixels was employed to prepare a final risk assessment of the resources, which have been exposed to the landslide. The results showed that 27–42 and 6–15 km2 were predicted as the very high and high LS zones using four ensemble metaheuristic machine learning algorithms. Based on the area under curve (AUC) of ROC, GR-ANN (AUC-0.905) appeared as the best model for LS modeling. The areas under high and very high landslide hazard were gradually increased over the progression of time (26 km2 at the 2 year return period and 40 km2 at the 100 year return period for the high landslide hazard zone, and 6 km2 at the 2 year return period and 20 km2 at the 100 year return period for the very high landslide hazard zone). Similarly, the areas of danger pixel also increased gradually from the 2 to 100 year return periods (37 km2 to 62 km2). Various natural resources, such as scrubland, built up, and sparse vegetation, were identified under risk zone due to landslide hazards. In addition, these resources would be exposed extensively to landslides over the advancement of return periods. Therefore, the outcome of the present study will help planners and scientists to propose high precision management plans for protecting natural resources, which have been exposed to landslides.