In C2, the strategy from C1 of modifying the DN originally captured by the sensor to new values has continued. These new values are referred to as Q_cal, or ‘Quantized and Calibrated Standard Product Pixel Values’. Previously, and still today for most sensors with similar characteristics, data captured by the sensor were scaled to DN at the sensor through a known linear function. In other words, radiances (L, typically in W/(m²·sr), or in W/(m²·sr·μm) if the spectral magnitude, L_λ, is provided) were written as integer values, using the appropriate radiometric resolution to avoid precision losses. Consequently, the inverse linear function, with slope a and Y-intercept b, could also be provided to retrieve L:

Usually, a and b were kept constant over time, until a further recalibration.

However, in the new Landsat collections this is no longer the case except for Landsat 1-7 level 1 (https://www.usgs.gov/landsat-missions/using-usgs-landsat-level-1-data-product): DN are rewritten to Q_cal to be able to provide rescaling factors from Q_cal to exoatmospheric reflectance, or “Top of the Atmosphere” (TOA) reflectance (R_TOA hereafter) in the level 1 product (containing Q_cal level 1), or at ground level, or “Bottom of the Atmosphere” (BOA) reflectance (R_BOA) in the level 2 product. Obviously, Q_cal level 2 values written in the level 2 product are different from the Q_cal values in the level 1 product.

The rescaling factors to reflectance (not to radiance), REFLECTANCE_MULT_BAND_# (the slope) and REFLECTANCE_ADD_BAND_# (the Y-intercept) (where # is the spectral band number), are adjusted to be always the same over time and even for all spectral bands in level 1 for Landsat 8-9 and for all platforms that are processed to level 2 (Landsat 1-3 are excluded: see section 3.3). It is important to note that, in the MTL files, different rescaling factors are documented with the same names, but in different “sections” (called groups): to TOA in the LEVEL1_RADIOMETRIC_RESCALING group (USGS, 2020, 2023), and to BOA in the LEVEL2_SURFACE_REFLECTANCE_PARAMETERS group (USGS, 2023).

For instance, under the section tagged as:

GROUP=LEVEL2_SURFACE_REFLECTANCE_PARAMETERS

rescaling factors are written as follows:

However, the drawback of this strategy in level 1 is that the rescaling factors to L (RADIANCE_MULT_BAND_# and RADIANCE_ADD_BAND_#) are no longer constant between dates, as they were in previous Landsat products.

To understand the meaning of this change, level 1 product is considered. Although the reasoning is analogous for the level 2 product, conceptual differences inherent to level 2 must be considered.

In Equation 1 by Pons et al. (2014), R_TOA can be computed in the same way as USGS applies in its product, following:

where d is the Earth-Sun distance in astronomical units and E₀ is the exoatmospheric solar irradiance, typically in W/m². Moreover, for simplification, this computation ignores the cosine of the solar zenith angle, assumes transmittances (τ) equal to 1 and discards any atmospheric correction [atmospheric radiance, L_a=0] since this is a TOA calculation.

The data thus calculated from R_TOA (which will theoretically be between 0 and 1, both included) are scaled between the minimum and maximum value defined by the USGS and, consequently, Q_cal can be obtained through a new linear transformation:

By isolating R_TOA we can deduce the rescaling factors that will allow converting Q_cal to R_TOA:

From Q_cal, it is also possible to obtain L by calculating a new slope and a new Y-intercept. However, this calculation will be dependent on d², which varies for each date. Note that this dependency is not the usual practice in most remote sensing data products.

Indeed, combining the Equation 2 and Equation 3 Q_cal can be rewritten as:

Then L can be rewritten as:

And finally, the rescaling factors can be computed as follows:

It can be observed that both the slope and the Y-intercept for obtaining radiances exhibit an annual cyclic shift depending on the date, as the Earth-Sun distance squared (d²) changes for each date. Figure 1 illustrates this behavior, reflecting the well-known variation of the square of the Earth-Sun distance by approximately ±3.3% throughout the annual cycle.

The existence of a simple transformation using the REFLECTANCE_MULT_BAND and REFLECTANCE_ADD_BAND parameters to obtain R_TOA from Q_cal values may seem practical. However, users must be aware that the same perfectly invariant surface under identical atmospheric conditions on different dates will yield different R_TOA values. This is because the coefficients provided in the MTL files do not account for the Sun angle, as explained in Equation 2 and stated in the USGS Landsat Level-1 Data Product documentation (https://www.usgs.gov/landsat-missions/using-usgs-landsat-level-1-data-product). Remember that these constant values are not available in Landsat 1-7 level 1 and a classical processing must be applied. Regarding level 2, BOA coefficients are constant in all spectral bands in TM, ETM+ and OLI sensors.

When analyzing an overlapped area in consecutive image pairs, it is common to find pixels that, with the same geographic position, exhibit different radiometric values between the two scenes. This inconsistency is already present in the level 1 product and will propagate to all subsequent products (level 2 and beyond). For instance, on June 10th, 2004, differences in pixels of interest from a Landsat 5 TM image of path 197 and rows 030 and 031 (WRS-2) showed discrepancies greater than 1 Q_cal in multiple spectral regions simultaneously. As examples, this pattern was found in the following pixels (longitude, latitude): snow (2.468796, 42.510419), cloud (2.519669, 42.495735), water (2.102762, 42.571561), coniferous vegetation (2.550657, 42.499102) and urban area (2.873622, 42.416712). Another example can be found in Landsat 9 OLI-TIRS at (1.165870, 42.510689) in Figure 2, in this case in path 198 and the same rows than the previous examples.

As no source of evidence to support the usage of the upper or lower scene is known from the authors of this paper, our suggestion is to produce the mosaic by averaging the pairs of values. This can be done, for example, using a layer calculator producing NoData when operating two cells if one of them, or both, are NoData, and having a special cover function, CVR(arg1, arg2) that properly manages NoData and data values in that way: the cell from arg1 is selected, unless that cell is NoData, then the cell from arg2 is selected. Through the expression CVR(([A]+[B])/2, CVR([A],[B])), the resulting cell will be the average if it can be computed, or, on the contrary the cell of the first or second argument (NoData is only obtained if the cell is NoData in both layers). This can be done, for example with the MiraMon CalcImg module (https://www.miramon.cat/help/eng/msa/CalcImg.htm) that also has a more high-level function CVRM(arg1, arg2), producing the same result with the more compact syntax CVRM([A], [B]).

When analyzing pixels from the same area captured a few days apart within the same year, and with no evidence of sudden changes in cover, it is expected that discrepancies should only be due to atmospheric causes, given that the illumination and observation angles are highly similar. For instance, when examining a pseudo-invariant fully concrete pixel in a flat area from an image of path 197 and row 031 (longitude: 2.226218, latitude: 41.410353), pixel differences reach 1.9 units of % R_TOA in the Near Infrared (NIR) between December 5th and January 3rd. This might be acceptable due to their TOA nature. However, this difference reaches 6.2 units of % R_BOA in the level 2 product, which is unexpected. Furthermore, when calculating the extreme differences in R_TOA and R_BOA for this pseudo-invariant pixel in the Blue, Green, Red, NIR, and Shortwave Infrared (SWIR) bands of 1.6 and 2.2 μm throughout the year, differences up to 16.1% (SWIR-2) in R_TOA, with an overall average of 13.0% across all bands, and up to 6.7% in NIR for R_BOA, with an overall average of 4.7%, can be observed (Table 1). Therefore, these discrepancies, along with those noted in subsection 2.6, suggest the potential benefits of employing atmospheric adjustment methods based on pseudo-invariant zones (Pons et al., 2014) or representative field radiometry data obtained through the use of UAV (Padró et al., 2018), as alternatives to the official product.

The task of radiometric homogenization between Landsat sensors that has been carried out up to the C2 is very interesting. However, when comparing the L values of identical images (same dates) from ESA-CEOS products, one can observe significant discrepancies in images from some instruments (Table 2). The pertinent question is whether these differences are justified by the temporal homogenization performed by the USGS. The largest differences were found in TM, while in ETM+ they were small. In our exploration, the coefficients in ETM+ were very similar or identical, while in Landsat 5 TM they were different. In ETM+, CEOS DN values were practically identical to C2 Q_cal, but this was not the case in TM scenes. Since the L_λ differences were almost always positive, it implies a systematic bias in which the C2 values are almost always higher than those of ESA-CEOS. While addressing this question is beyond the scope of this article, it is essential to highlight this potential issue to the Landsat user and producer communities.

Estimating surface temperature using the new image processing techniques in Landsat Collection 2 faces several challenges and limitations. Currently, data from the ASTER Global Emissivity Dataset (GED) are used to estimate the emissivity required for calculating surface temperature. However, various issues in its accurate estimation, such as blockiness artifacts due to nearest neighbour resampling in the GED product, missing time series data, or vegetation emissivity anomalies, result in the failure to estimate ST for certain pixels or images in the C2 ST product. Additionally, the surface temperature products may exhibit a “blocky” appearance over small surface targets with contrasting thermal infrared signatures. This artifact arises from the nearest neighbour interpolation of the 100-meter ASTER GED products to the finer 30-meter Landsat grid. Another challenge is the geometric mis-registration between the ASTER NDVI and Landsat data, which can cause visual artifacts in the surface temperature products, particularly in regions with significant thermal contrast. Furthermore, the ASTER GEDv3 dataset used for emissivity is static and represents average emissivity from 2000-2008, necessitating adjustments for current land surface conditions at the time of the Landsat overpass. These adjustments can introduce inaccuracies, making it difficult to achieve precise surface temperature estimates. To address these issues, alternative algorithms can be employed to calculate vegetation or bare soil emissivity from Landsat solar spectrum data, such as the method proposed by Sobrino et al. (2008). These and other problems are documented at USGS Landsat Collection 2 Known Issues. These limitations underscore the complexities involved in accurately estimating surface temperature with the new image processing methods in Landsat Collection 2 (https://www.usgs.gov/landsat-missions/landsat-collection-2-known-issues#ST).

As a new feature in C2, it is noteworthy that two layers are provided to indicate saturated pixels, which is beneficial for excluding them in precise radiometric approximations. Additionally, several issues related to the radiometric correction algorithm are documented at USGS Landsat Collection 2 Known Issues (https://www.usgs.gov/landsat-missions/landsat-collection-2-known-issues#SR), although a detailed discussion of these issues is beyond the scope of this brief article. It is, however, disappointing that, as acknowledged in the referenced source, the coding of NoData values is still conflated with zero radiometric values, which is unrelated to the inherent challenges of radiometric correction in certain situations. Nevertheless, the wide unsigned integer data type [0, 65535] used for distributing the imagery allows for a well-defined flag value (e.g., 65535), which is distinct from reasonable radiometric values.

USGS distribution platforms are excellent for discovering and downloading data, but they have ceased offering the classic products that are closest to raw data. While it is true that Landsat users rarely had access to completely raw data, less processed products allowed for improvements in georeferencing (using finer Digital Elevation Models, DEM) and, importantly, avoided the use of bicubic interpolation as the resampling method. This method is notably the least respectful of the original radiometry, which is particularly problematic in categorical thematic mapping applications. Mixing the radiometry of 16 pixels, even when done in a weighted manner, creates a suboptimal situation for classifiers, especially in areas with high radiometric contrast (e.g., sandy beaches and water).

On other occasions, very poorly georeferenced images are encountered (e.g., LT05_L1GS_198031_20060927_20200831_02_T2 in the Iberian Peninsula, as well as numerous images of Greenland, due to the poor quality of the DEM used). These images are labelled as ‘tier 2’ as a warning. However, in some cases, such as the aforementioned image, proper correction is entirely feasible, as demonstrated by previous processing carried out by the authors of this paper. Unfortunately, this is no longer possible as the original image is no longer available. While it is understood that the USGS cannot geometrically correct each image with the utmost precision, the current policy prevents users who are capable of doing so from achieving a higher level of accuracy in specific areas.

Additionally, in some regions, users have access to DEM and/or Digital Surface Models (DSM) that are superior to those used in C2. This is particularly critical for radiometric corrections: computing the cosine of the incident angle benefits from more detailed surface modelling, and using the same DEM/DSM for both geometric and radiometric corrections is crucial to avoid overcorrections and undercorrections in ridges and valleys.

Finally, in countries where the UTM projection is not the official one, or where a datum significantly different from WGS84 is used (and is relevant for a spatial resolution such as that of Landsat), image resampling can cause undesirable issues. This criticism is, in fact, acknowledged by the USGS: regarding the level 2 products of most of the United States, which are in the Albers conic projection, the USGS states: ‘Raw Landsat data are calibrated and reprojected directly into the U.S. ARD Albers projection, which means the U.S. ARD are not resampled twice’ (https://earth.esa.int/eogateway/documents/20142/0/Landsat-Collection-1-vs-Collection-2-Summary.pdf).

Numerous Quality Assessment (QA) metadata are still provided in C2 in image (band) format. However, the format remains non-trivial for many users due to the mixture of bits with different meanings within the same byte and the fact that some properties require the consideration of more than one bit at a time. Consequently, some software has developed ad hoc modules (e.g., ArcGIS: https://www.usgs.gov/landsat-missions/landsat-quality-assessment-arcgis-toolbox) or have separated specific bands for each quality assessment (e.g., MiraMon: Pons, 2025). These bands support significant compression, so distributing them in more user-friendly raster values would benefit the user community. It would also be preferable to codify pixels where the QA flags have no meaning with an explicit NoData value (e.g., all bits set to one).

Another issue relates to bands containing angles (sensor or Sun) (USGS 2022, p. 53). Problems were reported in Landsat Sensor Viewing Angle Bands Known Issues (https://www.usgs.gov/landsat-missions/landsat-sensor-viewing-angle-bands-known-issues), which have fortunately been resolved. However, it would be better to avoid coding NoData angles with 0, because this is a valid angle value, and to avoid coding angle values outside the sensed image area. Indeed, profile figures on the USGS Solar Illumination and Sensor Viewing Angle Coefficient Files webpage (https://www.usgs.gov/landsat-missions/solar-illumination-and-sensor-viewing-angle-coefficient-files) show these 0 values. Some software (e.g., MiraMon: Pons, 2025) recodes cells affected by NoData (deduced from other bands) to an appropriate value (e.g., –32768, as the angle files have units of hundredths of degrees, positive or negative) and explicitly documents this in the user metadata file. While angle files are included in C2 products, they are only for band 4 (near the center of the OLI focal plane). The angle tool has options to create angle bands for every band in the product, including the TIRS bands. The Solar Illumination and Sensor Viewing Angle Generation for Landsat 4-7 can be obtained at https://landsat.usgs.gov/sites/default/files/documents/LANDSAT_ANGLES_15_3_0.tgz. Detailed instructions, including how to obtain the source code are at https://d9-wret.s3.us-west-2.amazonaws.com/assets/palladium/production/s3fs-public/atoms/files/L4-7_AngleBands_README-20200520.txt.

For Landsat 8/9 the tool is at https://landsat.usgs.gov/sites/default/files/documents/L8_ANGLES_2_7_0.tgz and the instructions and code at https://www.usgs.gov/media/files/landsat-8-angles-creation-tools-readme. Landsat TM and ETM+ Solar and View Angle Generation Algorithm Description Document (ADD) is available at https://d9-wret.s3.us-west-2.amazonaws.com/assets/palladium/production/s3fs-public/atoms/files/LSDS-1691_L4-5TM-L7-ETM_Solar-View-Angle-Generation_ADD-v1.pdf. Landsat 8 Solar and View Angle Generation Algorithm Description Document (ADD) is available at https://d9-wret.s3.us-west-2.amazonaws.com/assets/palladium/production/s3fs-public/atoms/files/LSDS-1928_L8-OLI-TIRS_Solar-View-Angle-Generation_ADD-v2.pdf.

Level 2 products are not available for Landsat 1-3 (USGS Landsat Collection, https://www.usgs.gov/landsat-missions/landsat-collection-2, and Crawford et al., 2023), resulting in a significant loss of information in the time series (lack of temporal completeness) during the years between Landsat 1 (1972) and Landsat 4 (1982).

C2 is distributed in Cloud-Optimized GeoTIFF (COG) format, which is a significant advancement. This recent Open Geospatial Consortium (OGC) specification (Masó, 2023) enables images to be accessed directly from servers without the need for prior preparation, while still allowing the same files to be read efficiently from local environments.

In the level 2 product, rescaling factors for L and R_TOA are included (see section 2.1). Upon considering their meaning, it can be deduced that they are not applicable to the image they document. Given the non-reversibility of the process of writing the new Q_cal proportional to the R_BOA using simple coefficients, their inclusion seems highly confusing and prone to errors, making them unsuitable for instructional purposes. For instance, it is not easy to understand why, in the same file, there are two values for the same key (e.g., REFLECTANCE_MULT_BAND_4=), or a value (e.g., RADIANCE_MULT_BAND_4=) that should not be applied to that product.

In other words, the appropriateness of the “GROUP=LEVEL1…” in level 2 metadata files and the avoidance of the inclusion of radiance and reflectance rescaling factors, as they cannot be applied, could be reconsidered. Although the USGS (2023) states that ‘these fields are for provenance’, the risk of confusion is significant