The separation of fresh snow, exposed glacier ice and debris-covered ice on glacier surfaces is needed for hydrologic applications and for understanding the response of glaciers to climate variability. The end-of-season snowline altitude is an indicator of the equilibrium line altitude of a glacier and is often used to infer the mass balance of a glacier. Regional snowline estimates are generally missing from glacier inventories for remote, high-altitude glacierized areas such as High Mountain Asia. In this study, we present an automated, decision-based image classification algorithm implemented in Python to separate snow, ice and debris surfaces on glaciers and to extract glacier snowlines at monthly and annual time steps and regional scales. The method was applied in the Hunza basin in the Karakoram and the Trishuli basin in eastern Himalaya. We automatically partitioned the various types of surfaces on glaciers at each time step using band ratios combined with topographic criteria based on two versions of the Shuttle Radar Topography Mission elevation dataset. Snowline altitudes were extracted on a pixel-by-pixel basis using a “buffer” method adapted for each elevation dataset. Over the period studied (2000 to 2016), end-of-the-ablation season annual ELAs fluctuated from 4,917 m a.s.l. to 5,336 m a.s.l. for the Hunza, with a 16-year average of 5,177 ±108 m, and 5,395 m a.s.l. to 5,565 m a.s.l. for the Trishuli, with an average of 5,444 ± 63 m a.s.l. Snowlines were sensitive to the manual corrections of the partition, the topographic slope, the elevation dataset and the band ratio thresholds particularly during the spring and winter months, and were not sensitive to the size of the buffer used to extract the snowlines. With further refinement and calibration with field measurements, this method can be easily applied to Sentinel-2 data (5 days temporal resolution) as well as daily PlanetScope to derive sub-monthly snowlines.