Related code of the article (A generalized enhanced snow/ice index based on Landsat imagery for accurate snow/ice mapping in mountainous areas)

A minor issue was identified in the implementation of the terrain correction function. The code has been revised accordingly, and the corrected version is provided as follows: // Terrain correction ------------------------------------------function terrainCorrection(image) { // 1. Get solar zenith and azimuth angles from image metadata var solarZenith = ee.Number(90).subtract(image.get('SUN_ELEVATION')).multiply(Math.PI / 180); // radians var solarAzimuth = ee.Number(image.get('SUN_AZIMUTH')).multiply(Math.PI / 180); // radians // Convert solarAzimuth to an image to match the aspect dimension var solarAzimuthImage = ee.Image.constant(solarAzimuth); // 2. DEM data: slope and aspect var dem = ee.Image('USGS/SRTMGL1_003'); // SRTM DEM data var slope = ee.Terrain.slope(dem).multiply(Math.PI / 180); // slope in radians var aspect = ee.Terrain.aspect(dem).multiply(Math.PI / 180); // aspect in radians // 3. Calculate cosine of the incidence angle cos(i) var cos_i = slope.cos().multiply(solarZenith.cos()) // cos(θs)cos(θn) .add(slope.sin().multiply(solarZenith.sin()) // sin(θs)sin(θn)cos(φs - φn) .multiply((solarAzimuthImage.subtract(aspect)).cos())); cos_i = cos_i.where(slope.neq(0).and(solarZenith.neq(0)).not(), 0.0001); // assign value for invalid areas // 4. Process each band and correct non-QA_PIXEL bands var bandNames = image.bandNames(); // get all band names of the image var correctedBands = bandNames.map(function(bandName) { bandName = ee.String(bandName); // force bandName to ee.String // Check whether the band name is 'QA_PIXEL' var isQA = bandName.slice(0, 8).compareTo('QA_PIXEL').eq(0); // var isTIRS1 = bandName.slice(0, 5).compareTo('TIRS1').eq(0); // var cond= isQA.or(isTIRS1) // If it is QA band, return the original band return ee.Algorithms.If(isQA, image.select([bandName]), (function() { // Apply terrain correction for non-QA bands var band = image.select([bandName]); var regression = cos_i.addBands(band).reduceRegion({ reducer: ee.Reducer.linearFit(), geometry: image.geometry().intersection(dem.geometry(), 1), // limit to intersection of image and DEM scale: 30, maxPixels: 1e11 }); // Extract intercept a and slope b var intercept = ee.Number(regression.get('offset')); var slopeVal = ee.Number(regression.get('scale')); var C = intercept.divide(slopeVal).abs(); // calculate C factor and ensure non-negative // Apply SCS+C terrain correction formula return band.multiply((slope.cos().multiply(solarZenith.cos()).add(C))).divide(cos_i.add(C)); })() ); }); // Merge bands into a single image using ee.Image.cat var correctedImage = ee.ImageCollection(correctedBands).toBands() .rename(image.bandNames().map(function(bandName) { return ee.String(bandName); })); // Copy properties var finalImage = correctedImage.copyProperties(image, image.propertyNames()); return finalImage;}