Snow cover extent and distribution over the years have a significant impact on hydrological, terrestrial, and climatologic processes. Snow cover mapping accuracy using remote sensing data is then particularly important. This study analyses Landsat-8 NDSI snow cover datasets over time and space using different NDSI-based approach. The objectives are (i) to investigate the relation snow-NDSI with different environmental variables, (ii) to evaluate the accuracy of the common NDSI threshold of 0.4 against in-situsnow depth measurement and (iii) to develop a method that optimises snow cover mapping accuracy and minimises snow cover detection errors of omission and commission. Landsat-8 snow cover datasets were compared to ground snow depth measurements of climate stations over Switzerland for the period 2014–2020. It was found that there is a consistent relationship between NDSI values and land cover type, elevation, seasons, and snow depth measurements. The global NDSI threshold of 0.4 may not be always optimal for the Swiss territory and tends to underestimate the snow cover extent. Best NDSI thresholds vary spatially and are generally lower than 0.4 for the three snow depth threshold tested. We therefore propose a new spatiotemporal NDSI method to maximize snow cover mapping accuracy by using a generalized linear mixed model (GLMM). This model uses three environmental variables (i.e., elevation, land cover type and seasons) and raw NDSI values and improves snow cover mapping accuracy by 24% compared to the fixed threshold of 0.4. By using this method omissions errors decrease considerably while keeping a very low value of commission errors. This method will then be integrated in the Snow Observation from Space (SOfS) algorithm used for snow detection in Switzerland.