Simplified Automatic Atmospheric Correction for THEMIS Infrared Data

Thermal infrared spectral data from the Thermal Emission Imaging System (THEMIS) on Mars Odyssey can provide diagnostic information to characterize Martian surface compositions. Accurately isolating the atmospheric contribution from the thermal radiance of THEMIS data is essential for surface mineralogy analysis. In this work, we present improved semi‐automatic and fully automated tools to simplify the traditional atmospheric correction method for THEMIS data. Both methods use the Thermal Emission Spectrometer surface emissivity data but in different ways. These two methods are applied to several locations that have been well documented in previous studies, to assess their capability of determining surface spectral features. The fully automatic method demonstrates its simplicity, viability, and robustness by which more consistent results have been generated for the same area across multiple THEMIS images and for areas with large topographic gradients. The tools developed in this work will broaden the utilization of THEMIS infrared spectral data to the wider planetary science community.

olivine-rich materials (Figure 1a). The emissivity spectrum derived from radiance data without the atmospheric transmission removal (Figure 1b) is significantly different from the emissivity spectrum after complete atmospheric correction (Figure 1c), which displays a clear olivine absorption feature at ∼11 μm. The emissivity spectra after the atmospheric correction will help to investigate the surface mineralogy correctly.
Martian radiance measured in the TIR is the sum of several factors: surface emission that is attenuated by the atmosphere, atmospheric emission, downwelling atmospheric emission reflected from the surface, and atmospheric scattering. Among them, radiation reflected from the surface as well as atmosphere scattering can be neglected because the surface is usually warmer relative to atmosphere for THEMIS daytime data. Then the observed radiance by THEMIS for a single pixel can be simplified as: where obs( ) is the measured radiance for each spectral channel ( ) and pixel ( ) ; ( ) is a multiplicative term consisting of surface emissivity and atmospheric attenuation of surface radiation; [ surf , , ] is the Planck radiance at the highest brightness temperature ( surf ) of the surface throughout the THEMIS bands; and ( ) is the radiance due to atmospheric emission (Bandfield, Rogers, et al., 2004).
A method for atmospheric correction of THEMIS daytime IR data to remove atmospheric attenuation and emission was described in detail by Bandfield, Rogers, et al. (2004). The first step is the calculation and removal of the atmospheric emission term. The radiative contribution of the atmospheric emission can be calculated from regions within a THEMIS image that are spectrally uniform with variable temperatures, such as the sunlit (warm) and shaded (cold) slopes of a surface of equivalent composition, because the relative contribution of the atmospheric emission is greater over cold surfaces.
The second step of THEMIS atmospheric correction is the calculation and removal of the atmospheric attenuation term. A typical correction method relies on the surface emissivity of a small training area within a THEMIS image derived from TES data to determine atmospheric attenuation. The determined atmospheric attenuation contribution can then be removed on a pixel-by-pixel basis from the entire image or a portion of the image, with the major assumption being that the atmosphere is invariant across the area being corrected. Variable atmospheric conditions are expected, however, in areas where large elevation variations (>1 km) are present in a THEMIS image and also over large spatial scales (e.g., 100s of km). Although the constant atmosphere assumption usually spans over a small subset of full THEMIS image (∼1,000-2,000 lines), inaccuracies will result if applied over a large portion of or even the entire THEMIS image using only one atmosphere condition. As a result, the atmospheric correction must be applied to small regions of similar elevation in a THEMIS image.
In addition, the determination of the atmospheric attenuation component over a small area is not an objective or automatic process; presently, the user must manually choose the training area that is used to obtain the atmospheric condition for THEMIS data correction. Usually, a spectrally uniform area within the THEMIS image that is near the area of interest in elevation and spatial distance would be selected. Potential spectral changes would be produced if different training areas are selected by different users. The need for user experience and/or intervention has likely reduced the accuracy of quantitative comparisons between THEMIS images of different locations.
In this work, we present an improved semi-automatic tool to simplify the processing of the traditional THEMIS atmospheric correction (Bandfield, Rogers, et al., 2004) and a fully automatic technique for extracting THEMIS surface emissivity at a larger spatial scale. This work will benefit the wider planetary science community to investigate the science of Martian surface mineralogy more easily using THEMIS data.

Methodology
THEMIS data are provided as Reduced Data Records through the Planetary Data System (PDS), which are radiometrically calibrated radiance data derived from the THEMIS Experimental Data Record . Before atmospheric correction, we applied a series of previously-described data processing steps to remove noise and artifacts, which include the correction for time-dependent focal plane temperature variations, low frequency temperature drift from beginning to end of image acquisition, temperature variations across the calibration target, row and line correlated noise, as well as the projection of the images (Bandfield, Rogers, et al., 2004;Edwards et al., 2011) ( Figure 2). As mentioned above, the first step of THEMIS image atmospheric correction is the calculation and removal of the atmospheric emission. This process is documented in Bandfield, Rogers, et al. (2004) and has been further refined as described by Edwards et al. (2011). It is part of a conventionally used THEMIS processing pipeline (Edwards et al., 2011), and has been successfully applied in many previous studies (e.g., Bandfield et al., 2013;Pan et al., 2021;Rogers & Nazarian, 2013). The whole process has been written with Davinci, which is an open-source software package developed and maintained by Arizona State University (ASU) (https://davinci.mars.asu.edu/), and these steps requires no additional modifications and these data are archived as a part of the NASA/THEMIS PDS (https://static.mars.asu.edu/pds/).
With the availability of THEMIS data after standard processing stored on the ASU PDS Imaging Sub-Node, we developed a semi-automatic tool and a fully automatic technique written with Davinci to conduct the atmospheric attenuation correction.

Semi-Automatic Method
We developed a semi-automatic tool that provides a simplified way to select a TES-covered region as the training area for determining the atmospheric attenuation component.

Search High-Quality TES Data Overlying a THEMIS Image
For the targeted study area within a THEMIS image (Figure 3a), the tool first acquires all high-quality TES data over the THEMIS image through Davinci via a Simple Query Language Postgres database interface (a geospatially enabled database of all TES data acquired, including geometry, telemetry, and spectra, which is leveraged by JMARS and our tool to provide community access to TES data). The quality criteria for TES data include Orbit Counter Keeper (OCK) prior to 7,000, surface temperature >255 K, emission angles <10°, albedo values <0.2, atmospheric dust opacities <0.15, and water ice opacities <0.08 (Table 1). Then the program converts the coordinates in latitude/longitude of each TES footprint to the corresponding column/row values in the THEMIS image based on the image-provided projection. A blended topography map derived from High-Resolution Stereo Camera (HRSC) and Mars Orbiter Laser Altimeter (MOLA) at 200 m/pixel (Fergason et al., 2018) is used to provide the elevation of each pixel in the THEMIS image ( Figure 3b).

Provide Geographic and Topographic Constraints to Select the Training Area
The ideal TES pixels for atmospheric correction should be less than 1 km topographic difference and within 1,000 lines from the targeted study area. The TES pixels which meet these criteria are then provided (Figures 3c  and 3d). Then users need to select a small spectrally bland region (a compositionally homogeneous region, indicated by the region having the same hues in a range of THEMIS DCS stretches) from these TES covered areas to be used as the training data ( Figure 3e).   (Edwards et al., 2011). Improvements to the atmospheric attenuation correction step (red box) are the focus of this work.

Atmospherically Correct THEMIS IR Image With Atmospherically Corrected TES Spectra
The surface emissivity of each selected TES pixel is calculated and averaged to provide the surface emissivity of the training area. The calculation of TES surface emissivity follows the most common technique for TES atmospheric correction. Since the Martian surface and atmospheric spectral components have been demonstrated to combine in a linear manner , a linear least squares unmixing algorithm with a library of atmospheric components  and surface minerals (e.g., Rogers & Aharonson, 2008) can be used to retrieve the spectral response caused by the atmospheric components, and then their contributions are removed. The spectral library of atmospheric components and surface minerals is given in Table 2. With the availability of surface emissivity of the training area, it is converted to THEMIS 10-band spectral resolution. This averaged THEMIS emissivity is divided from the radiance of the training area to retrieve the atmospheric component, and this atmospheric component is applied to the targeted study area.
Overall, this semi-automatic tool only requires users to provide the coordinates in longitude/latitude or x/y pixels for the targeted study area and the training area. The output results of this tool are THEMIS surface emissivity with removal of atmospheric components for a small area within 1,000 lines from the center of the targeted area. For cases that may not be possible to select the TES data with similar elevation (the difference of elevation smaller than 1,000 m) of targeted area within THEMIS scene, the users may choose to give up the atmospheric correction, relax some of the elevation and distance constraints, or select the TES data within an area with more variation of elevation to the targeted THEMIS scene, accepting likely higher uncertainties in the derived surface emissivity.   (2000) Dust, high CO 2 Bandfield, Christensen, and Smith (2000) Water ice cloud (high latitude) Bandfield, Christensen, and Smith (2000) Water ice cloud (low latitude) Bandfield, Christensen, and Smith (2000) CO 2 Bandfield, Christensen, and Smith (2000) H 2 O Bandfield, Christensen, and Smith (2000) Surface Note. TES, Thermal Emission Spectromete.

Fully Automatic Method
In contrast to the traditional correction method that applies one atmospheric condition to a THEMIS image, this fully automatic method generates a variety of atmospheric contributions, and each atmospheric condition is only applied to the local area that is directly covered by the TES footprint within a THEMIS image to minimize the effects of variable atmospheric conditions.

Find All Acceptable TES Data Overlying a THEMIS Image
The fully automatic method first finds all relative high-quality TES data overlying a THEMIS image. In order to maximize the spatial coverage of TES data over a THEMIS image, lower TES quality criteria is allowed (OCK <10,000, albedo values <0.25, atmospheric dust opacities <0.25, and water ice opacities <0.15) ( Table 1) by default but the user may further restrict these parameters as desired.

Atmospherically Correct THEMIS IR Image in TES Covered Areas
In this step, the THEMIS image is masked to the TES footprint covered areas (Figure 3f). Then an atmospheric component is generated for each TES footprint one by one and removed from the underlying THEMIS pixels within each footprint. In other words, each TES footprint only acts as the training area for the exact same region itself, without the need to apply the atmospheric condition beyond the local training area. For areas in a THEMIS image that are targeted by multiple TES footprints, the final surface emissivity will be averaged.
This fully automatic method allows THEMIS-based quantitative spectral studies to a larger spatial scale without subjective selection of TES covered training area. However, a major difference in the overall process is that the training area is not spectrally bland, and incorporates all surfaces in the TES footprints, including the regions of interest. Given the near-global coverage of high-quality, clear atmosphere TES data when examining a 100 m/ pixel map and the global coverage of warm, clear atmosphere THEMIS data, a majority of the THEMIS images will have some section that is possible to automatically correct with no user intervention. However, given that TES data is discontinuous and typically covers sections within a THEMIS image, the semi-automatic method may be required for some locations on the planet.

Results and Discussion
The methods described above were applied to six locations (ranging from olivine-rich units to quartz-bearing terrain) that have been well analyzed by previous studies. For each location, we examined two THEMIS images over the same area to cross validate the viability and robustness of our fully automatic method. At the end of this section, we also discuss the uncertainties in our derived THEMIS surface emissivity spectra.
A THEMIS image I38360010, presented as the example above to illustrate the procedures of atmospheric correction methods, and another image I35989015 covering Jezero crater have been corrected with the traditional and automatic methods, respectively, and the results are shown in Figure 4. For the traditional method, we chose two different training areas to conduct the atmospheric correction for the image (Figures 3c-3e). Both training areas meet the selection criteria (within 1,000 lines and 1 km elevation difference from the target area) and have high-quality TES data available. While the derived surface emissivity from different training areas (black and blue spectra) are comparable with each other, small variations are noticeable that could lead to different interpretations. THEMIS emissivity spectra of the same area in two images (solid and dashed spectra) show large differences using the traditional method. Unit 1 represents an olivine-rich area, and the spectra of this unit should have lower emissivity value at band 7 (∼11 μm). Spectra of Unit-1 from image I35989015 indeed show an emissivity minimum at band 7 while same unit spectra from image I38360010 have an emissivity minimum at band 5 (∼9.3 μm) when corrected using training area 1. Therefore, the selection of training area is an important factor for the traditional THEMIS correction method.
In contrast, THEMIS surface emissivity calculated from the automatic method (red spectrum) shows more consistent results between images. Additionally, the surface compositions can be well characterized using the automatic method regardless of the training area selection (Figures 4b, 4d and, 4f). Olivine-rich materials represented by Units 1 and 3 show a significant absorption at band 7 for the automatic method compared to those from the traditional method. Spectral variations from Units 2, 3, and 4 within a TES pixel are also more prominent using the automatic method compared to those of the traditional method. Together, the automatic correction approach shows great capability to accentuate spectral features to well characterize the surface composition without the need of training area selection.

Gale Crater
Gale crater, the landing site of Curiosity rover, preserves a ∼5 km tall central mound of layered deposits changing from phyllosilicates/sulfates to sulfates/oxides from bottom to top of the mound . The central mound was formed and eroded by a combination of fluvial, lacustrine, and aeolian processes (e.g., Kite et al., 2013;Thomson et al., 2011). The stratigraphic trends suggest a paleoclimate transition from wet to dry throughout Mars' history. A large dune field composed of olivine-basaltic sand Olivine spectral features from the dune field are likely better characterized by the automatic method compared to the traditional method ( Figure 5). Similar to previous work by Rogers and Bandfield (2009), surface emissivity spectra of the dune field extracted from the traditional method (Unit 1, black spectrum) in this work show a slightly deeper absorption at band 7 while the emissivity minimum is located at band 5. However, spectra of the same area of the dune field corrected by the automatic method (Unit 1, red spectrum) show deep absorption and emissivity minimum at band 7, which unambiguously indicates the presence of olivine. Spectra from two small areas (Units 2 and 3) along the rover traverse have also been characterized. These two units are well covered by two TES pixels, respectively. Unit 3 shows a stronger olivine spectral feature relative to Unit 2 in image I34877001, but this difference is not apparent in image I50292001. The spectral variations between Units 2 and 3 are more prominent when using the automatic method for both images. For all three areas, spectra from the automatic method show better consistency between two images compared to those of the traditional method.
Spectra of the olivine-rich unit have already been shown to have a significant absorption at band 7 based on the traditional atmospheric correction from both THEMIS images ( Figure 6). However, THEMIS spectra of the same unit are more consistent between the two images using the automatic method compared to the traditional method. In addition, the olivine spectral feature is more evident using the automatic method, expressing a relatively higher emissivity value difference between band 6 (∼10.21 μm) and band 7 (∼11 μm).

Northeast Syrtis Crater
The presence of crystalline quartz has only been found in a few isolated exposures in the north/northeast of Syrtis Major region to date (Bandfield, 2006;Bandfield, Hamilton, et al., 2004). The spatial association of quartz with hydrated silica and nearby phyllosilicates indicates the quartz may have formed as a diagenesis product of amorphous silica instead of an igneous material (Smith & Bandfield, 2012;Smith et al., 2013). One of the quartz exposures is located in a crater in northeast Syrtis Major. Two THEMIS images I10770004 and I11082005 covering the crater have been analyzed here.
Spectra of the quartz-bearing unit show comparable results from the traditional and automatic atmospheric correction methods between the two images ( Figure 7). The relatively smaller emissivity value difference between band 4 (∼8.5 μm) and band 5 (∼9.3 μm) from the automatic method is still noticeable compared to the traditional method, which indicates stronger feature showing a higher bulk content of silica.

Nili Patera
Nili Patera is located at the summit of the low-relief Syrtis Major volcanic complex, which preserves a wealth of volcanic history and compositional diversity. The caldera floor contains basaltic rocks and dunes, as well as a volcanic lava flow composed of dacite . Several hydrated silica deposits have been identified on the flanks of a volcanic cone associated with the dacite unit using CRISM data (Skok et al., 2010). Small exposures of feldspathic rock have also been detected within the caldera floor bedrock by CRISM data (Rogers & Farrand, 2022;Wray et al., 2013).
THEMIS spectra of the dacite unit match well with each other from both the traditional and automatic methods ( Figure 8). We note that the emissivity value at band 7 from the automatic method is slightly lower than that of the traditional method, indicating the presence of a small amount of olivine, probably originating from the nearby olivine-basaltic caldera floor .
For the feldspathic rocks, THEMIS spectra are more consistent between the two images using the automatic correction method compared to the traditional method ( Figure 8). The band minimum position shifts to the lower wavenumber, indicating lower silica content compared to the nearby dacite. This result agrees with findings by Rogers and Nekvasil (2015) and Rogers and Farrand (2022), suggesting a silicic composition is not likely for these plagioclase-dominant units. In addition, the band contrast from the automatic method is much deeper than that of the traditional method. It supports the idea that these feldspathic units are composed of large plagioclase crystals (Rogers & Farrand, 2022;Rogers & Nekvasil, 2015).

Valles Marineris
Olivine-enriched bedrocks identified in Ganges and Eos Chasma in Valles Marineris are likely one of the largest continuous compositionally uniform units found on Mars by characterizing its lateral and vertical occurrence Edwards et al., 2008). These olivine-enriched outcrops are generally exposed on the canyon floor close to the canyon walls, representing in-place rocky materials as determined from morphological and thermal inertia analysis. The spatial and temporal distribution of olivine-enriched bedrocks in Ganges and Eos Chasma potentially indicate that a significant volume of volcanic lavas erupted to form extensive olivine-enriched flood basalts in early Martian history (Edwards et al., 2008). The olivine-enriched unit on the canyon floor identified by Edwards et al. (2008) shows consistent spectral features with a relatively deeper absorption at band 7 when corrected with the traditional method using training area 1 for both THEMIS images, while it exhibits remarkable variations when corrected using training area 2 between two THEMIS images ( Figure 9). Training area 1 shares a similar elevation with the olivine-rich unit on the canyon floor, so the atmospheric condition can be assumed to be constant for the atmospheric removal ( Figure 9e). However, a large elevation difference (∼4,000 m) between the canyon floor and training area 2 (Figure 9e) induces higher uncertainties in the derived surface emissivity among THEMIS images although spectra from THEMIS image I47751001 indicate a strong olivine feature (Figure 9g). In contrast, for the fully automatic method, spectra of the olivine-enriched unit on the canyon floor show a clear absorption at band 7 and are nearly identical between both THEMIS images (Figures 9b, 9d and 9g).
This olivine-enriched unit on the canyon floor has been characterized as materials carried down by mass wasting based on the morphological features, however, no olivine spectral feature was found on the associated canyon walls previously (Edwards et al., 2008). Here a small unit from the canyon walls shows a consistent spectral feature with an evident absorption at band 7 from the automatic method for both THEMIS images, representing olivine-rich materials (Figure 9). The spectral signature of the wall unit corroborates the land slide scenario for this olivine unit on the canyon floor. The emissivity spectra calculated from the traditional method express significant variations between two THEMIS images using both training areas because of the large elevation difference from the wall unit to either training area (Figures 9e and 9g).
For a targeted study area that shows large topographic gradients, such as the canyon walls in this example, high uncertainties are created in the derived surface emissivity when using the traditional method. However, the automatic correction method is able to produce more consistent results across THEMIS images regardless of the topographic gradients.

Uncertainties
The fully automatic correction algorithm shares the same processes with the traditional method in terms of the standard processing of THEMIS data. Therefore, the errors associated with the THEMIS data are the same with those of the traditional method, which have already been well described by Bandfield, Rogers, et al. (2004), and Edwards et al. (2011). The fully automatic method also uses TES surface emissivity data to determine atmospheric components but in a different way from the traditional method, so the uncertainties of the corrected THEMIS surface emissivity highly depend on the accuracy of TES surface emissivity retrieval. The intrinsic errors associated with TES spectral data have been described by Smith et al. (2000) and Bandfield, Hamilton, and Christensen (2000), and the uncertainties of TES surface emissivity retrieval from non-negative linear unmixing model have also been explicitly characterized by Rogers and Aharonson (2008). Another potential error that has never been mentioned in the previous works is the manual selection of the training area in the traditional method. As we demonstrate in this work, different selection of the training area could impact the surface emissivity derivation. Unlike the traditional method to choose the training area subjectively, the fully automatic method directly identifies the overlapping TES and THEMIS data and provides the leverage of the TES correction to the specific location where TES data is available without extrapolation, it results in THEMIS surface spectra more consistent with TES surface emissivity as opposed to situations where the atmospheric properties are extrapolated beyond the extent of the TES data, so it could avoid the potential error caused by different training area selection from different users. Finally, it is challenging to estimate the uncertainties of the corrected emissivity values because there are no "actual" values of the surface emissivity. In summary, the uncertainties of the corrected surface emissivity from the fully automatic method are mainly associated with the uncertainties of THEMIS and TES data, as well as errors related to TES surface emissivity calculation. Given the more consistent surface emissivity values derived from the fully automatic method across multiple images, it is likely that the uncertainties of the fully automatic method are relatively low.
From all the examples shown above, the automatic atmospheric correction method for THEMIS data we developed in this work illustrates the significant utility of this methodology to better characterize surface spectral features when compared to the traditional method. The surface emissivity spectra extracted from this automatic method exhibit great consistency across multiple images and large topographic gradients. The fully automatic method also generates more accentuated spectral features, which likely represents more accurate surface compositions. Given the advantages of the fully automatic method as shown in this work, the fully automatic method is preferred to the semi-automatic method, while the semi-automatic method should be used instead when the fully automatic method cannot be applied to some areas that are not covered by acceptable TES observations.

Conclusions
In this work, we have developed a semi-automatic tool to simplify the traditional atmospheric correction method for THEMIS infrared spectral data. This tool provides options of quality control and image visualization for users to select TES covered training areas more easily. The users only need to provide the coordinates in longitude/latitude or x/y pixels of the target area and the training area, and the program can conduct the atmospheric correction for the THEMIS data in an automatic way.
We have also developed a fully automatic technique to correct the atmospheric effects for THEMIS data. This method treats each TES pixel as the training area to isolate the atmospheric component for the local area that is directly covered by the TES footprint within a THEMIS image. Without the need to select a training area for the entire or a portion of THEMIS image, this method is able to automatically calculate the surface emissivity for each TES covered area throughout a THEMIS image. For areas that are not covered by any TES pixels, the semi-automatic tool can be used instead.
We have applied the semi-automatic and fully automatic methods to several THEMIS images at six locations representing a wide range of volcanic materials. The results demonstrate the advantages of the fully automatic method undoubtably in three aspects: (a) more consistent surface emissivity can be extracted for the same area across multiple images; (b) more consistent surface emissivity can be extracted for areas across large topographic gradients; (c) a more easily implemented tool can be used to benefit the wider planetary science community. The automatic method also shows great capability to accentuate spectral features at both large and small regional scales, indicating that the fully automatic method possesses the potential to characterize surface compositions more accurately.

Data Availability Statement
All command tools have been developed by Davinci and archived at Zenodo.org https://doi.org/10.5281/ zenodo.6648739 (Ye et al., 2022). The scripts are also available at Davinci website (http://davinci.asu.edu) as a part of the standard distribution maintained by Arizona State University.