Heatwaves can greatly impact societies, underscoring the need to extend current heatwave prediction lead times. This study investigated multiple machine-learning (ML) model approaches for heatwave occurrence prediction with long lead times of one to five months based on explanatory atmospheric and land surface features. Five ML classifiers were built using Google Earth Engine remote sensing datasets to predict heatwaves at national scale (Sweden) based on 16 features referring to the period of 1989–2019. Extreme Gradient Boosting performed best for lead times of one month (F1-score = 0.63, accuracy = 0.81) and four months (F1-score = 0.54, accuracy = 0.79), while K-Nearest Neighbour was best for lead times of two, three and five months (respective F1-score = 0.63, 0.65, 0.49, accuracy = 0.77, 0.79, 0.78). When applying the SHapley Additive exPlanations technique for model interpretation, land surface features emerged as more impactful heatwave predictors than atmospheric features at longer lead times. More frequent heatwave occurrence was associated with places in Sweden characterized by lower values of geopotential height, latitude, topographical slope, evaporation, precipitation and cropland area, and higher values of average temperature, mean sea level pressure, and southerly and westerly winds. The study also concretely exemplifies how use of this multi-model ML method can enhance predictions and further step-wise improve them, thereby facilitating earlier warning in support of better planning of measures to mitigate adverse heatwave impacts, up to several months ahead of their possible occurrence.