Algorithms

Overview

Atmospheric Correction

The surface reflectance derivation from Landsat and Sentinel-2 images uses Land Surface Reflectance Code (LaSRC) originally developed by Eric Vermote (NASA/GSFC) (Vermote et al., 2016, 2018) and adapted by USGS Earth Resources Observation and Science (EROS) Center for operational use. The LaSRC performance has been evaluated in Doxani et al., 2018, 2023. A detailed description of the operational Landsat surface reflectance and its ancillary layers are provided by the U.S. Geological Survey.

LaSRC derives surface reflectance with the aerosol optical thickness retrieved from the images assuming a continental aerosol model and using ancillary data from multiple sources. The ancillary atmospheric ozone and water vapor data for gaseous absorption correction are from the MODIS 0.05° Climate Modeling Grid (CMG) data (before May 2024) or from VIIRS 0.05° CMG (after May 2024). Standard sea level atmospheric pressure is adjusted for 0.05° resolution local elevation in molecular (Rayleigh) scattering correction. The aerosol optical thickness retrieval uses the spatially explicit ratios of the red to blue surface reflectances, also in the 0.05° CMG grid, that are characterized from 10-year MODIS observations; the spectral ratio data are adjusted by the observed finer-resolution Landsat and Sentinel-2 data in spectral bands less sensitive to the aerosol effects (Vermote et al., 2016, 2018). The HLS L30 product also includes the two thermal-infrared bands from the Landsat 8/9 TIRS sensor. The data values for these bands are not atmospherically corrected, i.e., are brightness temperature.

A true color composite (bands 4-3-2) of a Sentinel-2B image with top-of-atmosphere (TOA) reflectance (left) and LaSRC atmospherically-corrected surface reflectance (right). The same contrast stretch was applied to both images. Acquired on November 17, 2017 over NASA GSFC and surrounding areas (tile 18SUJ).

Cloud Masking and Quality Assessment

HLS provides per-pixel mask of cloud, cloud shadow, snow, water, and aerosol optical thickness levels. In earlier versions of HLS, the cloud and shadow masks were a union of results from multiple sources, but in HLS V2.0, they were generated exclusively by Fmask 4.7 (Zhu and Woodcock, 2012, Zhu et al., 2016), an update of the previous Fmask version reported in Qiu et al. (2017) and Qiu et al. (2019). Fmask was found to be one of the better performing algorithms when assessed on a variety of test data sites (Skakun et al., 2017, 2022). The Fmask cloud and cloud shadow output is dilated by 150 m during HLS processing and the dilated area is labeled as “adjacent to cloud/shadow.” The data user can treat the dilated area as clear if data loss is a concern. The mask labels and the dilation (cloud, shadow, snow/ice, water, adjacency) are stored in bits 1-5 in the per-pixel quality assessment (QA) band. The QA band also records, in bits 6 and 7, the aerosol optical thickness level (low, medium, or high) derived during atmospheric correction.  The QA band is named after the cloud masking algorithm Fmask, but contains both Fmask output and the LaSRC aerosol optical thickness level information. The use of pixels with high aerosol optical thickness is not recommended.

Spatial Co-registration

HLS has adopted the Universal Transverse Mercator (UTM)-based Military Grid Reference System (MGRS) into which ESA grids the Sentinel-2 L1C data.  Although the USGS Landsat data are also in UTM projection and the HLS output pixel size is 30 m, the input 30-m Landsat data can not be simply copied over because for the Landsat data the UTM coordinate origin corresponds to the center of a 30-m pixel, while in the MGRS system it corresponds to a pixel corner. That is, even when the input Landsat data and an MGRS tile are in the same UTM zone, there is a 15-m shift between the reference points. The HLS resampling of Landsat spectral data uses the cubic convolution method. If the input Landsat scene and the MGRS tile are in two adjacent UTM zones, reprojection of Landsat is needed before resampling. The resampling of the accompanying categorical QA values only considers the inner 2 x 2 pixels in the cubic convolution kernel, and the presence of any label there will turn on the corresponding output QA bit. Because the 2 x 2 pixels may have different cloud mask labels, an output HLS pixel may have multiple bits set to 1 in its QA byte; for the aerosol optical thickness level, the label of the highest aerosol level in the inner 4 pixels is chosen for output.

The 10-m, 20-m Sentinel-2 spectral data are resampled to 30 m with an area-weighted average, and pixels in the 60-m bands are evenly split into four 30-m pixels. The cloud mask labels (derived at 20 m) and the categorical LaSRC aerosol thickness levels (derived at 10 m) are resampled to 30 m from the overlapping input pixels with a rule similar to the Landsat categorical data resampling.

BRDF Normalization

Landsat OLI and Sentinel-2 MSI are narrow field-of-view sensors, but the view angle effect is still noticeable in the contrasting backward and forward view image pairs acquired on the same day or a few days apart. HLS processing normalizes this view angle effect by predicting a nadir-view reflectance with a Bidirectional Reflectance Distribution Function (BRDF)-based c-factor technique from Roy et al. (2016, 2017). The end result is 30-m Nadir BRDF-Adjusted Reflectance (NBAR) – for Landsat-derived HLS this is the final L30 product, whereas Sentinel-2 NBAR undergoes a subsequent bandpass adjustment to produce the final S30 product.

The HLS BRDF normalization uses a set of constant BRDF coefficients, derived from 12-month MODIS 500-m global BRDF product (MCD43A1 Version 6; Schaaf et al., 2002; Shuai et al., 2011). The derived BRDF coefficients are applied to OLI and MSI bands equivalent to MODIS ones. For the normalization of the Sentinel-2 MSI red-edge bands (705, 740, and 783 nm) that have no MODIS equivalents, the linearly interpolated BRDF coefficients from the enclosing MODIS red and NIR wavelength bands are used (Roy et al 2017). The technique has been evaluated using off-nadir (i.e. in the overlap areas of adjacent swaths) ETM+ data (Roy et al. 2016) and MSI data (Roy et al. 2017). The BRDF coefficients for the three kernels (isotropic f_iso, geometric f_geo, and volumetric f_vol) of the Ross–Li semiempirical model are shown in Table 1. The kernel definitions are described in the ATBD of the MOD43A1 product (Strahler et al. 1999) and the operational implementation for MODIS BRDF/albedo product are described in Schaaf et al., 2002.

Table 1: BRDF coefficients used in the HLS V2.0 c-factor normalization (Roy et al. 2016, 2017; Ju et al., 2025).

Band name	HLS L30 Band Code	HLS S30 Band Code	Equivalent MODIS band	fiso	fgeo	fvol
Coastal/Aerosol	B01	B01	3	0.0774	0.0079	0.0372
Blue	B02	B02	3	0.0774	0.0079	0.0372
Green	B03	B03	4	0.1306	0.0178	0.0580
Red	B04	B04	1	0.1690	0.0227	0.0574
Red-Edge 1	–	B05	–	0.2085	0.0256	0.0845
Red-Edge 2	–	B06	–	0.2316	0.0273	0.1003
Red-Edge 3	–	B07	–	0.2599	0.0294	0.1197
NIR Broad	–	B08	2	0.3093	0.0330	0.1535
NIR Narrow	B05	B8A	2	0.3093	0.0330	0.1535
SWIR 1	B06	B11	6	0.3430	0.0453	0.1154
SWIR 2	B07	B12	7	0.2658	0.0387	0.0639

The L30 and S30 NBAR are derived as:

where

and are the LaSRC-derived surface reflectance and the BRDF-normalized surface reflectance respectively; denotes the angle set of the observed solar zenith, view zenith, and the relative azimuth, and the angle set of the prescribed solar zenith, a 0° view zenith, and accordingly a 0°relative azimuth used in normalization. The c-factor is a ratio of the modeled surface reflectance for the normalization angle set to the modeled surface reflectance for the observation angle set.

During BRDF normalization, the solar zenith angle is slightly adjusted to account for the ~30 minute overpass time difference between Landsat and Sentinel-2. For an image of a given day gridded to a given MGRS tile, the solar zenith angles for that day at the tile center’s latitude are calculated for the Landsat and the Sentinel-2 overpass time respectively, based on their nominal equatorial overpass time (Li et al, 2018).  The mean of the two solar zenith angle values is used for all the pixels in the tile on that day for NBAR derivation. The amount of solar zenith angle adjustment is small.  And the adjustment has the additional benefit of normalizing the local solar zenith difference for the swath sidelap from the same type of sun-synchronous orbits. This is clear in the Sentinel-2 A/B tandem fly available since March 2025. In this configuration, 10 minutes after Sentinel-2B overpasses, Sentinel-2A overpasses in the west, with substantial swath sidelap between the two, forming an image pair of contrasting backward and forward views. The observations in the sidelap are made at slightly different solar zenith angles (at slightly different local solar time) because they are not at the nadir. This solar zenith angle setting for NBAR contrasts with the use of a latitude dependent but temporally constant solar zenith angle in earlier versions of HLS.

Bandpass Adjustment

The bandpass adjustment in HLS V2.0 is intended to correct for the effects of small differences in the relative spectral response (RSR) functions between Landsat OLI and Sentinel-2 MSI. The OLI/MSI common spectral bands use the same spectral filters, but there are subtle RSR differences. Although the TOA data calibration by NASA/USGS and ESA has minimized the difference between OLI and MSI data based on measurements on a few sites,  the difference may be greater over other distinctive land cover types and propagate into surface reflectance. HLS uses the OLI spectral bandpasses as references, to which the MSI NBAR data are adjusted. The bandpass adjustment on the Sentinel-2 NBAR uses a linear relationship derived for data in the common bands simulated from hyperspectral reference data. A total of 500 hyperspectral surface reflectance spectra from 160 globally distributed EO-1 Hyperion scenes covering diverse land cover types (vegetation, desert, urban, and snow) were selected and used to simulate the MSI and OLI bands by convolution with the respective OLI/MSI RSRs. The formula to derive the OLI-like MSI NBAR (ρ_OLI) is given in (3), where ρ_MSI is the Sentinel-2 derived NBAR before adjustment, and band-specific slope (a) and intercept (b) adjustment coefficients are given in Table 2. The 30-m ρ_OLI  is the final S30 product.

Separate band-specific coefficients for the formula above are derived for Sentinel 2A, 2B and 2C. Once corrected, the spectral difference between the MSI and OLI bands is less than 2% and a standard deviation of residual error less than 0.005 reflectance units (Claverie et al., 2015, 2018).

Table 2: Bandpass adjustment coefficients

			Sentinel-2A	Sentinel-2A	Sentinel-2B	Sentinel-2B	Sentinel-2C	Sentinel-2C
HLS Band name	OLD Band number	MSI Band number	Slope (a)	Intercept (b)	Slope (a)	Intercept (b)	Slope (a)	Intercept (b)
CA	1	1	0.9959	-0.0002	0.9959	-0.0002	1.0030	-0.0000
Blue	2	2	0.9778	-0.0040	0.9778	-0.0040	0.9851	-0.0027
Green	3	3	1.0053	-0.0009	1.0075	-0.0008	1.0038	-0.0009
Red	4	4	0.9765	0.0009	0.9761	0.0010	0.9718	0.0011
NIR 1	5	8A	0.9983	-0.0001	0.9966	0.0000	0.9995	-0.0003
SWIR 1	6	11	0.9987	-0.0011	1.0000	-0.0003	0.9994	-0.0007
SWIR 2	7	12	1.0030	-0.0012	0.9867	0.0004	0.9910	0.0004

References:

• Vermote and Kotchenova, 2008: https://doi.org/10.1029/2007JD009662 
• Vermote et al., 2016, 2018: https://doi.org/10.1016/j.rse.2016.04.008, https://doi.org/10.1109/IGARSS.2018.8517622 
• Doxani et al., 2018, 2023: https://doi.org/10.3390/rs10020352, https://doi.org/10.1016/j.rse.2022.113412 
• Zhu and Woodcock, 2012, https://doi.org/10.1016/j.rse.2011.10.028 
• Zhu et al., 2016: https://doi.org/10.1016/j.rse.2019.05.024 
• Qiu et al. (2017) https://doi.org/10.1016/j.rse.2017.07.002 
• Qiu et al. (2019) https://doi.org/10.1016/j.rse.2019.05.024 
• Skakun et al., 2017, 2022 https://doi.org/10.1016/j.rse.2022.112990 
• Gao et al., 2009
• Roy et al. (2017) https://doi.org/10.1016/j.rse.2017.06.019 
• Roy et al 2017 https://doi.org/10.3390/rs9121325 
• Vermote et al., 2009
• Li et al., 2018: https://doi.org/10.1109/TGRS.2018.2885967 
• Zhang et al., 2016: https://doi.org/10.1109/TGRS.2015.2480684 
• MCD43A1 Version 6; https://doi.org/10.5067/MODIS/MCD43A1.006 
• Schaaf et al., 2002; https://doi.org/10.1016/S0034-4257(02)00091-3 
• Shuai et al., 2011): https://doi.org/10.1016/j.rse.2011.04.019 
• Strahler et al. 1999: https://lpdaac.usgs.gov/documents/97/MCD43_ATBD.pdf 
• Claverie et al., 2015, 2018: https://doi.org/10.1016/j.rse.2018.09.002, https://doi.org/10.1016/j.isprsjprs.2023.03.011 
• Ju et al., 2025: https://doi.org/10.1016/j.rse.2025.114723