1. Introduction
The Landsat program has a long heritage of acquiring, archiving, and distributing remotely sensed data [
1]. Starting in 2016, the US Geological Survey (USGS) adopted a tiered collection management structure for its Landsat data products that ensures a consistent set of processing for the Landsat archive within a given collection, while allowing a set of calibration updates to be performed between any two given collections [
2]. This collection philosophy provides the user community with a set of Level-1 and Level-2 products for the entire Landsat archive that is consistent not only in product format, but also in the radiometric and geometric algorithms used in their processing. A key part of integrating the Landsat 9 satellite and its sensors into this structure, from a radiometry perspective, involves the datasets collected during the Landsat 9 (L9) underfly of Landsat 8 (L8). The nearly time-coincident acquisitions of the imagery help alleviate time-varying components, such as atmospheric and temporal changes, that impact the ability to otherwise compare products from the two satellites more directly from a radiometry perspective. Similarly, this coincident imagery also provides better co-registration between the products generated from the two systems by eliminating the impacts of those same time-varying components from the differences that may arise during registration and measuring the systematic corrected imagery to the ground control. Once these temporal differences and changes are reduced, other aspects of the differences between the systems such as sun and sensor viewing angles can become more pronounced in the direct comparisons between the two systems. This paper describes some of the aspects of the underfly from an acquisition perspective, the co-registration between the Level-1 products produced from the two systems, and the driving forces in the differences between sensor viewing and sun angles associated with the products, giving some examples of these differences. Within this paper, these differences are only addressed for the multispectral data of the Operational Land Imager (OLI) instruments.
1.1. Operation Land Imager Architecture
The Operational Land Imager is made up of 14 sensor chip assemblies (SCAs), staggered in the along-track direction, where each band is aligned with respect to the even and odd SCAs in the across-track direction. This includes minor clocking of the SCAs to compensate for optical distortion.
Figure 1 shows the SCA staggering and band ordering within an SCA of the OLI instrument [
3].
This staggering in the along-track direction produces a noncontiguous geographic sampling across track within the odd and even SCA groupings, while remaining contiguous within an SCA across the 494 detectors in each band of a SCA. This along-track staggering can be represented as the number of nominal multispectral lines between the trailing and leading SCAs or the odd and even numbered SCAs, for each band in the along-track direction, and the corresponding timing delay when aligning the SCAs geographically.
Figure 2 shows this SCA staggering within an image, caused by the focal plane displacement of each band within each SCA. For
Figure 2, by not trimming the leading and trailing portions of the imagery acquired at the same time from the SCAs across bands during a given image acquisition in the output product, while also listing timing differences between the leading and trailing SCAs, the geographic displacement for a single frame of imagery acquired can be shown along with the timing differences between when the leading and trailing bands within an SCA acquire data within a collection. The timing in the table, therefore, represents the maximum relative timing delay of the odd and even SCAs for a given band. An item that is worth noting from
Figure 1 is that the staggering and placement of the SCAs indicate that the nominal nadir pointing of the instrument does not truly lie within any given SCA but lies within the grouping of the SCAs placed within the center of the focal plane.
This time delay between the sampling of the even and odd SCAs means that the sun angles and sensor angles will produce discontinuous values across track within the imagery. For sun angles, these timing differences produce small differences. Even though the timing differences between two acquisitions for the L9 and L8 underfly imagery are larger than those within scenes due to the staggering of the SCAs, the sun angle differences are still considered small between the two acquisitions. For the sensor viewing angles, especially for the azimuth viewing angles, which are dependent on the focal plane SCA placement over the geographic area acquired, discontinuities within an image can be large, as demonstrated in the analysis results presented within this paper [
4].
1.2. Sensor and Sun Angle Viewing File
Within this paper, a Level-1R (L1R) refers to a dataset that has been radiometrically corrected but does not have any geometric corrections applied to the imagery. Level-1T (L1T) refers to imagery that is radiometrically and geometrically corrected and has corrections applied for the terrain parallax effects within the imagery using a digital elevation model (DEM). This L1T referenced within this paper can be either a Level-1 Terrain and Precision (L1TP) or a Level-1 Systematic Terrain (L1GT) image. The L1TP image has terrain effects removed from the imagery and has biases present in the spacecraft telemetry data, position, and attitude, removed on the basis of measurements made between the imagery and the Landsat Global Survey (LGS) ground control library. The L1GT has the terrain effects removed from the imagery but the data do not have the spacecraft biases associated with its telemetry data removed. When speaking to Landsat 9 and 8 products within this paper, many can be part of both the L1TP and the L1GT products; where it is unnecessary to designate between the two product types, L1T is used, whereas, in cases where it is necessary to differentiate, and a subject can only be referred to with either the L1TP or the L1GT nomenclature, that designation is used.
Each Landsat product is accompanied by a solar illumination and sensor viewing angle coefficient file that contains sensor viewing and sun angle information for that product. The procedure for using the file contents to compute the angular relationships is performed through the following steps: (1) calculating an image time (per pixel) by determining a Level-1R line number, which can be directly related to time, using a rational polynomial relating Level-1T product coordinates to the corresponding Level-1R line number; (2) evaluating an ephemeris model that provides spacecraft position in an Earth-centered Earth-fixed (ECEF) coordinate system using the Level-1R line number as a proxy for time as an input; (3) evaluating a sun direction model that provides solar coordinates as an ECEF unit vector given a Level-1R line/time.
The rational polynomial relationship between the Level-1T and Level-1R image coordinates is shown in Equation (1). The height in the rational polynomial equations is the height at the given Landsat product pixel location, extracted from a DEM. The equations listed below, Equations (1) and (2), are determined from the spacecraft and instrument telemetry along with the spacecraft and instrument characteristics such as the instrument focal plane design. The USGS provides the user community the tools, using the equations listed in this paper, that are needed to generate their own set of a sensor viewing and sun angles for Level-1 products. [
4]. For the equations listed below, the nomenclature of L1T refers the Level-1T image product, either L1GT or L1TP. The subscripts Line and Sample refer to the Level-1T product line and sample pixel locations within the output product. The MeanLine and MeanSample subscripts are the mean line and sample pixel locations within the L1T product. The rational polynomial coefficients given in Equation (1), a
n, b
n, c
n, and d
n, allow the user to calculate time by determining the L1R pixel location of a given L1T product line and sample pixel location. The determination of the L1R pixel location allows the rational polynomial coefficients e
n and f
n to be used to determine the satellite’s position or the sun’s position for a given L1T output pixel location. This produces a flow of first determining the L1R pixel location and time for a given L1T output pixel location. Once time is determined, the sun’s location and the satellite’s position can be found for a given L1T output pixel location. This information can then be used to determine the corresponding viewing and sun angles for that L1T output pixel location. The satellite and sun’s position calculated resides in an Earth-centered Earth-fixed coordinate system with units of meters.
where
L1
TL =
L1
TLine −
L1
TMeanLine,
L1
TS =
L1
TSample −
L1
TMeanSample,
L1
TLine is the output pixel line location in
L1
T,
L1
TSample is the output pixel sample location in
L1
T,
Hgt =
Height −
HeightMean,
L1
R is the Level-1R pixel location,
L1
T is the Level-1T pixel location,
Height is the elevation at the
L1
T pixel location,
L1
TMeanSample is the mean sample pixel location in
L1
T,
L1
TMeanLine is the mean line pixel location in
L1
T,
HeightMean is the mean elevation in
L1
T, and
a0 to
a4,
b1 to
b4,
c0 to
c4, and
d1 to
d4 are rational polynomial coefficients (RPCs) of the model.
Equation (1) provides rational polynomials defining Level-1 terrain corrected image (L1T) line and sample locations to input Level-1R geometrically raw (uncorrected) line and sample locations. By determining the Level-1R geometrically uncorrected pixel location, the time for a given L1T line and sample location can be calculated.
The equation defining the satellite position on the basis of the Level-1T and Level-1R image, with the Level-1R location determined by Equation (1), is shown in Equation (2). Equations for the other satellite coordinate axis locations (
Saty and
Satz) and the sun’s corresponding cartesian coordinates (
Sunx,
Suny, and
Sunz) follow the same format as below for the satellite’s X (
Satx) location.
where
SatX is the satellite Earth-centered Earth-fixed
X-coordinate,
SatXMean is the satellite Earth-centered Earth-fixed mean
X-coordinate,
L1
TL =
L1
TLine −
L1
TMeanLine,
L1
TS =
L1
TSample −
L1
TMeanSample,
L1
TLine is the output pixel line location in
L1
T,
L1
TSample is the output pixel sample location in
L1
T,
Hgt =
Height −
HeightMean,
L1
RL =
L1
RLine −
L1
RMeanLine,
L1
RS =
L1
RSample −
L1
RMeanSample,
L1
R is the Level-1R pixel location,
L1
T is the Level-1T pixel location,
Height is the elevation at output L1T pixel location,
L1
TMeanSample is the mean sample pixel location in
L1
T,
L1
TMeanLine is the mean line pixel location in
L1
T,
HeightMean is the mean elevation in
L1
T, and
e0 to
e9 and
f1 to
f9 are RPCs of the model.
Equation (2) provides rational polynomials defining satellite location given a Level-1R, an image that is radiometrically corrected but geometrically uncorrected, and a Level-1T, an image that has terrain compensation and then may or may not contain corrections based on the application of ground control, t line, and sample location. These relationships define the satellite X-coordinate position. The equations for the satellite Y- and Z-coordinate positions follow the same format, with a different set of rational polynomial coefficients. The satellite’s coordinate location at a given time can in turn define the viewing angle for a given L1T output pixel location.
1.3. Landsat 9 Underfly
The Landsat 9 underfly refers to the period during the commissioning of the Landsat 9 satellite and its instruments that allowed L9 to slowly drift under the orbital track of the Landsat 8 satellite. This was observed during the commissioning period when a series of orbital maneuvers (termed “burns”) were needed to raise the orbit of the satellite to its final destination within the Worldwide Reference System 2 (WRS-2). [
5]. The exact timeframe during which the L9 satellite’s sensors overlapped with those of its L8 counterpart depends upon what is used to define an adequate amount of overlap between the systems for a given application. Typically, the defined time-frame for the underfly resides across 5 days where there was at least a 10% overlap between the fields of view of the instruments on the two satellites. However, since the amount of overlap also varies with latitude, even this 5 day window can be somewhat misleading, as represented in
Figure 3.
In
Figure 3, it is demonstrated that, at higher latitudes, there is more overlap between the field of view of the two satellites and its sensors than at mid and lower latitudes. At the beginning and end of the underfly period, overlap may only be present at high-latitude locations. To demonstrate the amount of overlap and number of scenes acquired during the underfly period
Figure 4 was assembled to show the percentage overlap with respect to the geographic location of the acquisitions. Blue in
Figure 4 represents the scenes with a limited amount of overlap, while red represents scenes with the maximum amount of overlap present in the image pairs.
Along with the percentage overlap between the two acquisitions, another aspect that can be important in data analysis is that, at the beginning of the underfly, as well as at the end of the underfly, the time difference between the two acquisitions was greater than at the heart of the underfly, i.e., the time in the middle of the underfly date range. These time differences lead to small differences in the sun angles associated with the acquisitions. While those differences are not in the range of multiple days as would be the case in normal operations, but rather minutes, this change will nevertheless produce differences in the sun angles associated with the two acquisitions. The time differences between scenes relative to the L8 acquisition time are shown in
Figure 5 across the date range of the underfly. These time differences were based on times extracted from metadata information listing scene center times to the nearest whole second. The plot in
Figure 5 shows the difference in scene center times over the range of date of the underfly. Jumps and discontinuities in the plot are associated with the fidelity of the scene times being kept. Taking into account this granularity in measurements, the time difference between L9 and L8 changed by approximately 0.002209 s/s or 190.8623 s/day.
3. Results and Discussion
As discussed previously, the co-registration accuracy between the L9 and L8 imagery during the underfly is based on their common registration to the Landsat Ground Control Library GCPs. Updates to these GCPs throughout the collection phased releases of Landsat products and software have increased the accuracy of this reference dataset where an accuracy of better than 9.5 m (CE95) can be expected in areas of limited cloud cover and temporal changes between the search and reference imagery [
11,
12,
13,
14]. The co-registration between the L9 (search) and L8 (reference) image products should have essentially no differences due to temporal changes between the search and reference imagery, and any differences that may be present in the co-registration would be due to the regions in the imagery that are not part of the overlap between the fields of view, where ground control will be applied to a geographic area for one scene and not the other. These differences in the two scenes of ground coverage between the reference and search imagery lead to differences in the set of control points used and may impact the precision solution step of correcting for any systematic biases present in the spacecraft telemetry data. However, these differences between the two datasets produced in the precision step are expected to be small.
Table 1 shows the image-to-image (I2I) statistics generated on the IAS for each band of data with at least 50% land within the imagery, less than 20% cloud cover within the imagery, and at least 10% overlap between the L9 and L8 products. Statistics are listed in units of meters. A total of at least 388 scenes were used in the analysis. As individual band results were filtered for outliers, and as different bands will have differing image-to-image correlation results due to the spectral band differences and the land features acquired during imaging, the number of images used in calculating the results is not the same across all bands. The table shows the mean, standard deviation (STD), and root-mean-squared error (RMSE) given in units meters. From
Table 1, the imagery is co-registered to less than 2.2 m radial root-mean-squared error (RMSEr) for all bands.
Figure 6 shows the geographic distribution of the scenes from
Table 1. As shown in the figure, the scene list was well distributed geographically.
To investigate whether there is a dependency between the co-registration accuracy and the amount of overlap between the two images, the RMSEr based on the percentage overlap was analyzed, as shown in
Figure 7. Each image pair’s radial RMSE from I2I is shown along the
y-axis, while the percentage overlap between the images is shown along the
x-axis.
Figure 7 shows that the amount of overlap between the L8 and L9 images had essentially no impact on the co-registration of the imagery within the overlap region. This assessment shows that, regardless of the overlap within the two products, the registration between the two products will essentially be the same due to the use of a common, consistent control point reference.
Another indicator of the registration accuracy of the Landsat products is the results generated from the geodetic accuracy assessment. This step determines how well the imagery fits to the Landsat ground control both before (pre-fit) and after (post-fit) adjusting for any bias or rate errors present in the spacecraft’s telemetry for position or pointing knowledge. Although this indicator represents the relationship between the imagery and the ground control, rather than the previous I2I results, which demonstrates how well the imagery stacks together from a co-registration perspective, the geodetic accuracy results provide useful information to support the I2I results and can indicate the geometric accuracy of an L1GT image product, an image that has the effects of terrain removed but no ground control applied to the data. Using the Landsat 9 and 8 I2I scene list and determining the overall geodetic accuracy pre-fit and post-fit results along and across track for those scenes produced the results listed in
Table 2. The RMSE listed in the table is the RMS of the RMSE of all the scenes analyzed.
Table 2 shows that both sets of instruments (L9 and L8 OLI instruments) were very consistent with the Landsat Ground Control library. This is based on the low RMSEs for both the pre- and post-fit results. The pre-fit results also help indicate what could be expected for any imagery from the underfly time-frame which may not be able to achieve a product level of precision correction and is available in an L1GT state. A question regarding L1GT product registration accuracy, i.e., for scenes that do not go to an L1TP product, but fall back to an L1GT, is how well they will be registered to each other. Scenes that are distributed as L1GT products that either are in areas where no ground control points are available (e.g., Antarctica) or fail to reach a precision state due to a lack of good clear ground cover because of clouds, snow, or landscape changes with respect to the ground control used and the newly acquired data. Landsat 8 has shown very good registration accuracy without the ability to apply ground control, typically within the 18–22 m CE90 range [
15]. During the early commissioning of L9, the same time-frame as the underfly period, L9 tended to exhibit similar characteristics. To determine the co-registration of the L1GT underfly scenes, an analysis was performed by processing a group of L1TPs, 40 cloud-free scenes, to an L1GT state, and an I2I analysis was performed on these scenes produced with the two differing product types, comparing L1TPs to L1GTs. These 40 scenes, now in an L1GT state, were compared using the same I2I processes used to compare the L1TPs, i.e., the previous I2I analysis performed and shown in
Table 1. The L9 L1GTs were compared against the L8 L1TPs by performing I2I between these 40 image pairs. Those results were then compared against the pre-fit geodetic accuracy reports generated from the L9 L1TPs. This comparison was conducted to verify that the pre-fit geodetic results provided from the L9 L1TPs would be an indication of what the co-registration accuracy of an L9 L1GT would be to the L8 L1TPs. As the pre-fit geodetic accuracy results are a measure of the geometric offsets between the precision-corrected L1GT and the Landsat ground control, measuring the L9 L1GTs to the L8 L1TPs that were also registered to the same ground control should essentially give the same results. A plot of this relationship is shown in
Figure 8, where the L9 pre-fit geodetic accuracy radial offset (RMSEr) results are shown along the
y-axis, while the radial offset of the I2I results are shown along the
x-axis. As the reports for the geodetic accuracy are given in terms of the along- and across-track direction and the I2I results are in map projection, line and sample coordinates looking at the radial components will help alleviate differences that would arise from plotting along and across track with the map-projected
X and
Y offsets. These results should show good correlation between the measurements, as the radial offset in an ideal situation between the two measurements should be equivalent since we expect the geometric differences between any given L1GT and the Landsat ground control library to be equivalent to the difference between that L1GT and an L1TP over the same area. Differences in the sampling of the measurements made within the two scenes between the two processes, with differing outlier logic performed on the measurements, would introduce some differences between the I2I and pre-fit geodetic accuracy results. As can be seen in
Figure 8, the correlation between the pre-fit geodetic accuracy and I2I results for all bands was 0.93. This figure demonstrates that L9 scenes that may not achieve an L1TP state and fall back to L1GT are well represented by the pre-fit geodetic accuracy results shown in
Table 2. Along with the outlier logic differences stated previously, the I2I results would produce a better measurement in that the I2I characterization approaches an almost auto correlation between two datasets, with only small differences in radiometry and geometry being present between the imagery. Conversely, the geodetic accuracy involves the correlation between the ground control library and the L9 imagery acquired; this correlation step will also include any temporal and landscape changes. Although these temporal type differences should be minimized by the geodetic accuracy step filtering poor correlations based on the ability to produce a precision model, it is the I2I correlation step that will provide a purer measurement from a correlation perspective between the images.
The second analysis involving the L1GTs was to compare the L9 and L8 L1GTs themselves. The same 40 scenes used in the L9 L1GT to L8 L1TP comparison were used in the matching L1GT comparison.
Table 3 shows the I2I results between the L9 and L8 L1GTs.
With the co-registration being a little more than 2 m for bands in the underfly image L1TP pairs, making the underfly data well co-registered, one element that could play a prominent role in comparisons of the imaging pairs from a science perspective is the sun and sensor viewing angle differences. As stated previously, the sensor architecture of the Landsat 9 and 8 instruments makes the sun and sensor viewing angles for the images acquired complex. Due to the focal plane of the OLI instrument being made up of 14 independent linear arrays that together form the full across-track field of view, the sun and sensor viewing angles are not contiguous in time of data acquisition when one traverses from east to west, or vice versa, across an image product.
Figure 9 shows the L9 and L8 underfly scenes from WRS-2 037/037 acquired on 15 November 2021. The overlap between these two scenes was 64%, while the acquisition time difference between the scene center times was 179.759 s. The L9 sensor azimuth angles shown in the upper right of the figure demonstrate the SCA-to-SCA discontinuity between the set of linear arrays that make up the focal plane for a given band. This discontinuity, which is most visible in the figure for the sensor azimuth angles, is present in all angle bands generated. This discontinuity in angles is due to the angular differences for each band of each SCA and the timing differences between odd and even SCAs for when a given SCA acquires across-track geographic locations on the Earth surface. The SWIR-1 OLI band was used for calculating the angles shown in the figure.
As the two acquisitions would not be viewing the same geographic location on the Earth’s surface at the same time or with the same SCA or sensor angles, these SCA dependencies can cause differences in the two images’ sun and sensor viewing angles. These differences can be difficult to visualize when looking at the data within a given scene. This is due to the fact that, when viewing the sun angles for a given acquisition for a single scene, as shown in
Figure 9 and
Figure 10, using horizontal profiles of the data within a map-projected space in the angle band images, there is essentially a linear relationship across the acquisition with a small nonlinearity due to the viewing properties for each band of each SCA. To better show the nonlinear SCA-to-SCA dependencies for these sun angles for a given scene (WRS-2 037/037 acquired on 15 November 2021), the dominant linear relationship can be removed from the data, and the small angular differences between the SCAs can be shown for these sun angles.
Figure 10 shows both these horizontal profiles derived from the data, providing a single profile for each angle, which is dominated by the linear relationship across the data, and then a profile with the linear component removed, demonstrating the SCA-to-SCA dependencies.
To demonstrate how these SCA-to-SCA dependencies drive differences in the sun and sensor viewing angles for the underfly image acquisitions, the sensor azimuth angles are shown in
Figure 11 for path/row 037/037 acquired on 15 November 2021. This figure shows the sensor azimuth viewing angle bands for the two acquisitions projected to the output map projection space of the mosaic of the two images, the images of which were shown in
Figure 9a, and the profiles at the same location for the mosaic are shown next to the sensor azimuth images for both acquisitions (
Figure 11c). From
Figure 11, it can be seen that the SCA number within the focal plane for Landsat 9 and 8 OLI instruments was generally not the same for the same location acquired in the overlap region.
Figure 11 shows that, although the viewing geometry was the same for both instruments of the images and for both acquisitions, the difference in their field of view coverage produced differences in the satellite viewing and correspondingly the sun angles for the two acquisitions in the overlap region. In many ways, as the display of the profiles show in
Figure 11, however, the difference in these sensors’ viewing azimuth angles can be visualized as a simple translation of the angles across the west-to-east direction of the images (although this is not exactly the case) from one acquisition field of view to the other. Although the differences in sensor azimuth angle, as seen in
Figure 11, can be thought of as a simple shift between the two sensors, in reality, it is more complex due to the different SCAs viewing the geolocations on the ground at a slightly differing set of times.
Figure 12 shows the sun and sensor viewing angle differences within the overlap region of the two acquisitions over WRS-2 037/037 acquired on 15 November 2021. Statistics for the original single scene sun and sensor viewing angles and the difference in the sun and sensor viewing angles are shown in
Table 4. The image mosaic is shown in
Figure 12a, while the difference in the azimuth sun and sensor viewing angles is shown in
Figure 12b,c.
The difference images in
Figure 12 show the SCA-to-SCA dependencies and the differences that they cause between the two image’s sun and sensor viewing angles. As shown in
Table 4, the largest differences in magnitude occurred in the sensor azimuth angles. The differences in sun viewing angles were less than 1° between the two acquisitions.
Toward the end of the underfly time-frame, the Landsat 9 satellite was pointed off nadir during several descending passes to optimize the amount of geographic overlap between the L9 and L8 images, which would otherwise be minimal or even possibly nonexistent. On 16 November 2021, 57 descending rows were acquired as off-nadir viewing imagery. On 17 November 2021, 64 descending rows were acquired as off-nadir viewing imagery. These off-nadir angles ranged from −9.493° to −14.752°. This off-nadir pointing would impact the sun and sensor viewing angle geometry.
Figure 13 shows an L9 and L8 underfly pair, where the L9 spacecraft was pointed off-nadir to cover a larger geographic area of the L8 field of view. These images were acquired over WRS-2 path/row 172/041 on 17 November 2021. The sensor viewing and sun angles per pixel are also shown in the figure for the L9 image. The L9 azimuth angles are the most noticeably different from a nominal nadir viewing acquisition. As the L9 spacecraft was pointed off-nadir, it introduced an angular bias in the sensor viewing angles, which removed the sharp transition in the azimuth angles at the center of the field of view, when compared with nadir looking acquisitions, and removed the peak in the sensor zenith angles. This is demonstrated by the angle band profiles taken from the two azimuth and zenith angles produced for each acquisition, as shown in
Figure 14 and
Figure 15, respectively.
The sensor and sun viewing angle statistics in the overlap region for scenes shown in
Figure 15 are shown in
Table 5.