Background

The Python Hyperspectral Analysis Tool (PyHAT) is an open-source library developed by the U.S. Geological Survey (USGS) with the goal of providing a single free and easy-to-use source for tools and algorithms that scientists can use to analyze spectroscopic data. Although initially developed for planetary science, PyHAT includes tools that are broadly applicable for planetary and terrestrial remote sensing, as well as laboratory data. PyHAT includes a graphical user interface (GUI) that allows easy access to certain back-end capabilities of the library. The GUI is particularly useful for users who do not wish to develop their own code. The goal of PyHAT is to enable scientists to focus on their core scientific investigation rather than spend time and money on software development or on expensive closed-source software.

PyHAT development has been supported by three National Aeronautics and Space Administration (NASA) Planetary Data Archiving, Restoration, and Tools program grants and received support from the USGS in 2023. The first grant, led by Ryan Anderson, involved the development of analysis capabilities for point spectra and the GUI. A second grant, led by Lisa Gaddis, focused on the analysis of data from orbital imaging spectrometers, such as the Moon Mineralogy Mapper (M³; Pieters and others, 2009) or the Compact Reconnaissance Imaging Spectrometer for Mars (CRISM; Murchie and others, 2007). A third grant, led by Itiya Aneece, expanded on PyHAT capabilities by implementing additional dimensionality reduction, endmember identification, unmixing algorithms, and incorporating these capabilities into the GUI. Funding through the USGS, led by Travis Gabriel, included consolidating orbital and point spectra code; developing documentation, examples, and this user guide; establishing protocols for continuous integration, continuous delivery, and incorporating external contributions; and releasing PyHAT v. 0.1.0. PyHAT was originally designed for the analysis of laser-induced breakdown spectroscopy (LIBS) spectra from the ChemCam and SuperCam instruments on Mars rovers and has received support from those projects. However, updates to the code now ensure that these techniques are broadly applicable to other types of spectra as well, including orbital data used in remote sensing applications.

Motivating an Open-Source Spectral Analysis Tool

Spectroscopic data require specialized analytical techniques, which can be a barrier to their broader use by the scientific community. Scientists may be domain experts, such as having in-depth knowledge of reflectance spectroscopy phenomena, but may not be expert programmers. Furthermore, they may not have funding to pay for expensive software packages and licenses, limiting their ability to interpret spectral data. In the cases where the scientist is also a well-versed programmer, their tools are commonly closed source or are only available upon request, which can compromise reproducibility of their work, set up barriers to access, or create quid pro quo situations.

PyHAT GUI

Alongside the back-end set of modules in PyHAT, which involve a variety of relevant spectral data processing and analysis methods, we have developed the PyHAT GUI, which provides easy access to back-end functionalities for nonprogrammers. Commonly, spectroscopists have a strong visual intuition for their datasets. Therefore, the GUI includes several plotting tools for visual interpretation of spectra and analysis results. The back-end and GUI are installed as a single standalone package.

PyHAT Access

PyHAT is an open-source project hosted and managed by the USGS. The open-source nature ensures these capabilities are readily available to the entire scientific community. In the development of PyHAT, open-access Python libraries are leveraged where possible. When not available in Python, we have translated open-access scripts from other languages. The GUI uses PyQt5 to facilitate user interfacing, and data manipulation relies on the NumPy (van der Walt and others, 2011) and pandas (McKinney, 2010; pandas development team, 2023) libraries. Many of the underlying algorithms used for the analysis of spectra are sourced from the scikit-learn machine learning library (https://scikit-learn.org; Pedregosa and others, 2011) and the PySptools library (Therian, 2018). This guide provides a brief introduction to different methods and algorithms available in PyHAT and the GUI. Readers are referred to the official documentation of scikit-learn, PySptools, and (or) other underlying packages for more detailed information.

Data Format in PyHAT

PyHAT works with data in a simple comma separated value (CSV) tabular format (fig. 1) with each spectrum and its associated metadata stored together in a row. The CSV file has two-level column labels with the top row indicating broad categories of data (by default, “wvl” for wavelength, “meta” for metadata, or “comp” for composition¹ ) and the second row indicating more specific categories. For the “wvl” category, the underlying row would contain the wavelength values of the spectrometer and every row thereafter would contain spectral intensities at each wavelength. For “meta” this could include “target name” or other categories. For “comp” this could be “SiO2” or “Quartz” compositional categories for which the user has data. The tool uses the pandas library to efficiently read data in this CSV format into a data frame, making use of the multi-index capability to handle the two levels of column names. Using multi-indexed columns enables the software to easily access large blocks of data—for example, all the spectra in the data frame (all rows and only columns with “wvl”). Another example includes selecting specific metadata columns for each spectrum, such as SiO₂ content of all spectra, which is relevant to the analysis of laboratory spectra where the composition is known and included in the metadata. The use of pandas data frames is also helpful for users that access the PyHAT back-end in custom Python scripts, as it provides a familiar and widely used framework for data storage and manipulation.

Note that the column names must be unique. The tool does not handle spectra with overlapping spectrometers that have identical spectral channels; that is, for two columns, the first row reads “wvl” and the second row reads an identical numerical value for the spectral wavelength. Users should resolve nonunique wavelengths or other repeated column names outside of the PyHAT tool prior to loading data.

Users may specify the top-level column labels used to indicate spectral data, compositional data, and metadata. By default, spectral data are assigned the top-level label of “wvl,” metadata (if there are any) are assigned the top-level label of “meta,” and compositions (or other metadata types containing continuous numerical values to be used in regression) are assigned the top-level label of “comp.” Composition columns are not required if the user is not regressing the spectral data to compositional data.

Interfacing with Hyperspectral Cube Data

PyHAT code is designed to work with data stored in tables as described in the “Data Format in PyHAT” section. However, many planetary (for example, CRISM) and terrestrial (for example, Airborne Visible/Infrared Imaging Spectrometer [AVIRIS]) datasets are collected as hyperspectral cubes, rather than as point spectra. PyHAT includes tools for “flattening” a cube into the native PyHAT spectral object. This converts each pixel in the hyperspectral cube to a row in the PyHAT data table. The row includes the intensity for each spectral band (the spectrum) as well as the x and y pixel coordinates from the original cube. These flattened PyHAT objects can then be used to reconstruct the original cube or to produce output of a single column from the table as a two-dimensional image. Information about the geographic projection of the original data cube can also be stored in the same object with the spectra. PyHAT uses the Geospatial Data Abstraction Library (GDAL) for this functionality, allowing header and metadata information to be retained so that reconstructed cubes or images can inherit the correct projection. This has been tested for cubes from CRISM and M³ datasets. Other datasets may require separate dedicated input/output functionality so that they can be flattened into PyHAT format and retain projection information.

PyHAT SpectralData Object

PyHAT stores data in a Python class named SpectralData. This class is initialized by providing a pandas data frame in the format discussed in the previous section along with keywords to indicate the strings used to label metadata, compositional data, and spectral data as well as a name for the dataset and geographic projection data where relevant. The SpectralData object includes numerous methods that provide a simple way to apply certain functionality to the spectra stored in the object. This simplifies writing custom scripts that utilize PyHAT functionality, and the SpectralData object is also used where possible in the interface between the GUI and the back-end code.

Tool Capabilities

PyHAT’s capabilities are grouped into several categories: data management, preprocessing, classification, regression, and visualization. These categories are also reflected in the top-line menu of the GUI. The categories are briefly introduced here and are discussed in detail in their respective sections.

Data management determines what data will be analyzed and how they will be organized. This includes reading and writing data; looking up metadata; combining, dividing, and organizing datasets; identifying and removing outliers; and identifying endmembers.

Preprocessing includes all steps involved in preparing the spectra for analysis once the data have been ingested into the tool. Whereas data management deals with what data will be used, preprocessing makes changes to the values of the data. For example, spectral masking, baseline removal, calibration transfer, normalization, dimensionality reduction, and unmixing are preprocessing steps.

Classification refers to methods that are used to group similar data within a dataset. It is used in various terrestrial and planetary studies, including mapping land cover, plant species, and minerals. The ability to separate classes depends on spectral variability within and across classes. Classes can be defined by the user, which is necessary when performing supervised classification when the user already knows the features of interest. PyHAT does not currently have supervised classification routines implemented. Unsupervised classification, or clustering, is implemented in PyHAT and is useful when the user does not know which classes to expect in their dataset. With unsupervised classification, spectral samples are divided into a (typically user-defined) number of classes based on some metric that assesses their similarity. Currently the GUI includes the k-means and spectral clustering algorithms.

A major focus of the GUI is to enable users to perform regression analysis of their data. Regression is the prediction of a continuous numerical quantity (for example, the abundance of SiO₂ in a target) based on observed spectra. It is a supervised technique that relies on a training dataset for which the spectrum and the quantity in question, SiO₂ content in the above example, are both known. This “known” information is used to generate a statistical model that can accurately predict the quantity in question when presented with new spectra. The GUI includes numerous regression algorithms and the ability to run cross validation to ensure that the model parameters are tuned appropriately and that the model is not over-trained so that it performs well on novel data. The tool also allows users to blend the results of multiple submodels to ensure accurate results across a wide variety of targets. We provide an in-depth description of the submodel blending in subsequent sections.

Finally, visualization is essential to allow users to interpret their data using results from the GUI’s statistical analysis. The GUI includes interfaces to plot either rows or columns of the data and to visualize the results of dimensionality reduction, classification, and regression analyses.

Tool Interface

PyHAT can be used without the GUI by importing the back-end library into any Python script. For users who do not wish to write their own scripts, the GUI provides access to many PyHAT capabilities. The GUI is designed around the concept of “workflows”—sets of individual data processing and analysis steps that are applied in a specific sequence to achieve a result (fig. 2). Each step within the workflow is contained in a separate “module,” and these modules can be arranged in any sequence specified by the user, giving the user flexibility to make changes to the logical steps in their analysis. There are logical limitations to this flexibility. For example, principal component analysis (PCA) must be run before PCA results can be visualized, data must be loaded before they can be normalized, and so on.

The GUI is centered around a main window that contains the workflow modules. Every workflow begins with a module that is automatically added in the main window when the program opens asking the user to set the default path for output. The top-level menu items (or toolbar items for Macs) are workflow, data, preprocessing, classification, regression, and visualization. Selecting an entry from one of these menus adds the corresponding module to the workflow. There are several buttons along the bottom of the GUI: rerun, stop, insert after, delete, and run. These buttons are used to make changes to the workflow and to run the loaded modules. Above these buttons is a progress bar, and above the progress bar is a console window that provides output from running modules and reports any errors. Detailed descriptions of each part of the interface are provided below. Although the exact workflow needed will depend on the data and the task at hand, the end of the guide includes workflow examples that demonstrate how to use the GUI.

Version and Installation

This user guide is for version 0.1.2 of the PyHAT repository, which is available at https://code.usgs.gov/astrogeology/pyhat/-/releases/0.1.2.

Instructions for installation are included in the README.md file in the repository. For updated versions of the code, visit the master branch of the repository at https://code.usgs.gov/astrogeology/pyhat.

Community Support

We encourage users to become part of the PyHAT community through the submission of GitLab Issues. Using this functionality, users can report bugs, suggest improvements, and request assistance. Access to this functionality requires an account on https://code.usgs.gov and a USGS employee must make an account for you. You may contact PyHAT developers to start the process. Install instructions, external developer instructions, and guidelines are included in the repository.

Documentation

PyHAT implements an automated documentation generator that produces HyperText Markup Language (HTML)-based documentation. Documentation is generated based on the markdown files in the repository as well as in Python docstrings, which are large comment blocks at the top of most functions. Docstrings and the HTML documentation provides the user with a description of inputs and outputs, notes, and example use cases. Efforts to improve the HTML documentation are ongoing. Users are encouraged to use GitLab Issues to suggest documentation improvements, or users can assist in documentation directly by becoming a contributor to the repository.

The documentation is hosted at https://astrogeology.code-pages.usgs.gov/pyhat/build-docs/libpyhat.html#module-libpyhat. Alternatively, users can generate a local copy of the documentation using the following commands.

pip install -U sphinx nbsphinx sphinx-rtd-theme
sphinx-apidoc -o docs .
sphinx-build -b html docs public

About This Guide

This user guide is intended to provide an overview of the capabilities of the PyHAT GUI. The GUI serves as a front end for spectral data management, preprocessing, analysis, visualization, and a variety of statistical and machine learning algorithms enabled by the PyHAT library back end. This document describes the conceptual approach of these algorithms and their strengths and weaknesses, but it is beyond the scope of this guide to describe the algorithms in detail. In addition, some of the algorithms rely on abstract mathematical concepts, which are challenging to describe in general terms. Users should (1) refer to the references in this guide for more details on the algorithms and (2) take advantage of PyHAT’s capabilities to experiment with different algorithms and parameters to determine what works best for the task at hand. There is no guarantee that the available algorithms are appropriate for the user’s data, and it is the user’s responsibility to validate and properly interpret the results. This software is under continuous development and the screenshots in this guide may differ slightly from the most up-to-date PyHAT version.

Save Workflow

The “Save Workflow” option allows the user to save a workflow for future use. When this option is selected, a dialog window will pop up for the user to navigate to the location in which they want to save the workflow. Workflows are saved as JavaScript Object Notation (JSON) files, so they can be opened in a text editor, if necessary. Doing so is not required in standard PyHAT use.

Restore Workflow

If the user has previously created and saved a workflow of modules within the GUI, it can be restored using the “Restore Workflow” option. When selected, an “Open Workflow File” window pops up with which the user can browse to the workflow file. To avoid conflicts between multiple restored workflows, this option is disabled after one workflow has been restored. The user must close the PyHAT GUI and reopen it to load a second workflow.

New Window

The “New Window” option opens another instance of the GUI window. This can be useful when the user wants to work with multiple workflows or wants to experiment with a module before adding it directly to their main workflow.

Exit

The “Exit” option allows the user to exit the GUI, which can also be done by clicking on the “X” on the window. Workflows are not automatically saved, nor are the datasets or regression models. To save datasets, write to a CSV file using the “Save Data to CSV” option under the Data menu. Workflows and regression models can also be saved as described above and in the “Regression” section, respectively.

Load CSV Data

The Load CSV Data module allows the user to load a CSV file into the GUI. The user can either enter the path to the file directly into the File Name field or browse by clicking on the button to the right of the field. The user can then enter a name for the dataset. This name will be used to identify the dataset in subsequent modules of the workflow. The dataset name defaults to “Data.”

If the input file contains header rows that need to be skipped, the user can specify the number of rows to skip. Three checkboxes allow the user to indicate whether the file being read contains spectral data, metadata, and compositional data. When checked, a field appears where the user can specify the string used to label these columns in the dataset. If duplicate columns of spectral data (that is, identical wavelength values) are encountered while loading the data, the tool will keep the first (leftmost) column and will notify the user that a column has been removed in the loaded data product.

Save Data to CSV

The Save Data to CSV module allows the user to write output data and (or) results to a CSV file. The user first chooses which dataset to export. The module dynamically reads the top-level column labels from the selected dataset and uses them to populate the list of variables to write. By default, all variables are selected, but the user can change the selection to customize the export. For example, to reduce the file size it may be useful to export just the metadata and predictions for a dataset but not the actual spectra. The user can also specify the name of the file.

Load ChemCam Data

The Load ChemCam Data module allows the user to use the GUI to work with ChemCam (Maurice and others, 2012; Wiens and others, 2012) data, which are available on NASA’s Planetary Data System at https://pds-geosciences.wustl.edu/missions/msl/chemcam.htm. The first field in the module allows the user to specify the search string, and the second field specifies the search directory on the user’s computer. When run, this module will search recursively within the specified directory for files matching the search string and read them into the standard PyHAT format. The default search string “cl5*ccs*.csv” is set to identify CSV-formatted cleaned calibrated spectra (CCS) files with names like those on the Planetary Data System. ChemCam performs most LIBS observations in a series of shots over the same point, such that several points form a raster of points for a given target. The “Averages” and “Single Shots” radio buttons allow the user to choose whether to ingest the average spectra (one spectrum per analysis point) or the single shot spectra (one spectrum per laser shot per point).

After reading the data, the resulting data frame can optionally be written to the output directory as a CSV file. Whether the data are saved to CSV or not, they are also stored in the GUI so that they are available for analysis. The user can specify both the name of the dataset that will be used through the rest of the workflow and the name of the CSV file to which the dataset will be written.

Load SuperCam Data

The Load SuperCam Data module is nearly identical to the Load ChemCam Data module but reads Planetary Data System data from the SuperCam instrument instead. SuperCam data are available at https://pds-geosciences.wustl.edu/missions/mars2020/supercam.htm.

Rename Data

If the user wants to change what a dataset is called within the GUI, the name can be changed using the Rename Data module. Simply choose an existing dataset from the drop-down list and enter the new name into the text field.

Rename Column

Like the Rename Data module, the Rename Column module allows the user to select a new name for an existing column. This is particularly useful to clarify the meaning of a column before subsequent analyses automatically append more columns to the dataset.

Look Up Metadata

The Look Up Metadata module (fig. 5) provides the capability to match metadata from one dataset to metadata in another dataset. For example, this can be useful to match the target name in two different datasets, one that contains spectra for each target and the other that contains their corresponding composition. The spectral and compositional datasets can then be combined to be used in, for example, regression analysis. The user first chooses their datasets in the top-left and bottom-left fields. The user then selects the “Variable to Match” in both datasets; the drop-down list is populated with the existing metadata column names in the corresponding datasets. The columns selected in “Variable to Match” should contain unique strings or numbers that allow the tool to map from the metadata dataset to the primary dataset. It is the user’s responsibility to ensure that the values in the dataset column exist in the specified column of the metadata file. If they are not there, the program will leave those metadata fields blank. Likewise, it is the user’s responsibility to ensure that the identifiers are unique. If they are not, the module will populate each row in the dataset with the first match encountered in the metadata file, which could result in incorrect duplicated metadata and potentially compromise subsequent analyses.

Remove Rows

The Remove Rows module allows the user to remove rows from a loaded dataset by evaluating logical operations based on the values of certain variables (fig. 6). The interface for this module begins with a drop-down menu that allows the user to select which dataset to work with. Below the dataset selection menu are one or more rows of drop-down menus that represent individual logical operations. These logical operations dictate which rows to remove from the data. The leftmost drop-down menu in each row is populated with the metadata (“meta”) and composition (“comp”) column names from the dataset selected. Note that for large datasets, this module can be slow to update when a dataset is selected. This is because the module needs to read and parse the selected dataset to populate the leftmost drop-down menu. The next drop-down menu to the right contains a selection of mathematical operators. The third drop-down menu contains the unique values contained in the column selected in the leftmost drop-down menu. Each entry in this third drop-down menu includes unique column values, followed by the number of rows in the dataset that have the identical value. The fourth drop-down menu in each row is used to determine whether an additional filter should be applied to the data. This drop-down menu can be null (containing just a dash symbol) or can contain the word “and.” If “and” is selected, another row of drop-down menus will appear. This allows the user to combine a series of logical filters to remove specific rows in the dataset. When this module is run, the console window will print out the dimensions of the dataset before and after row removal. To examine the modified dataset, the user can then add the Save Data to CSV module to examine the CSV file.

Combine Datasets

With the Combine Datasets module, the user can combine two datasets at a time. The datasets are concatenated via the pandas library using the outer method, such that if one dataset has a column that the other does not have, that column will be retained in the concatenated dataset. For simplicity, the indices (the numerical values assigned to rows of data by pandas when creating data frames) in the original datasets will be ignored and a new index column will be created. If the indices contain information that the user wants to preserve, such as the original order of spectra in the dataset, that information should be stored as a metadata column instead. For more information on this function within pandas, the user can refer to the pandas user guide (pandas development team, 2023).

Copy Dataset

This module creates a duplicate of a dataset with a new name specified by the user.

Identify Outliers

With the Identify Outliers module, the user can identify outliers using the isolation forest algorithm or the local outlier factor algorithm. Outliers may be indicative of poor-quality data, such as in the case of instrument malfunction, anomalous pixels, no signal, and so on. Outliers may also represent valid data that stand out in some way, for example, spectra collected of material standards and compositional extremes. This module only provides the means to identify potential outliers; users should use their expertise to determine whether outliers should be removed from the dataset. The user can use the Split Data or Remove Rows modules within PyHAT or use an external program to remove the rows identified as outliers. The scikit-learn library offers several outlier detection algorithms, including robust covariance, one-class SVM, isolation forest, and local outlier factor. In scikit-learn documentation (https://scikit-learn.org/stable/auto_examples/miscellaneous/plot_anomaly_comparison.html), the local outlier factor and isolation forest algorithms are shown to outperform other outlier identification methods and thus were included in PyHAT.

Isolation Forest

The isolation forest algorithm (Liu and others, 2012) randomly selects features and a threshold value for splitting data into isolated observations (fig. 7). The length of the tree, that is, the number of random splits, needed to isolate the observation is an indicator of whether it is an outlier. Outliers have much shorter trees that are easier to isolate. For more details on this algorithm, refer to the scikit-learn documentation, available at https://scikit-learn.org/stable/modules/generated/sklearn.ensemble.IsolationForest.html. Within PyHAT, the user needs to define the number of estimators and the proportion of the data to treat as outliers. When this algorithm is run, two metadata columns are added to the dataset. One column contains the outlier score, which is a floating-point number that indicates how strongly a given spectrum registers as an outlier. Lower values are more likely to be outliers. The second column is the binary classification of the data into inliers (values of 1) and outliers (values of −1). The fraction of inliers and outliers is specified by the user, and the most appropriate value varies depending on the application.

Local Outlier Factor

The local outlier factor algorithm (Breunig and others, 2000) calculates the local outlier factor score for the k number of nearest neighbors around certain observations (fig. 8). The score compares local density with global density of observations. Observations in very low-density areas compared with the global density are potential outliers. Although no single selection for k is best for all situations, 20 works in many cases. In very noisy datasets, a higher value may be more appropriate. This algorithm performs well with high-dimensional datasets. For more information on local outlier factor, refer to the scikit-learn documentation, available at https://scikit-learn.org/stable/modules/generated/sklearn.neighbors.LocalOutlierFactor.html. For this method, the user needs to define the number of neighbors, leaf size, distance metric, and percentage of outliers. Leaf size refers to the number of points at which the algorithm changes its strategy for constructing the tree. Changes to this variable do not affect the results but can have a substantial effect on the execution speed and memory required by the process.

Split Data

With the Split Data module, the user can split a dataset by unique values of a certain variable, creating a new dataset for each value. This can be useful for manipulating data as an alternative to removing rows when the user wants to keep both sets of data. For example, the user may want to identify outliers, but retain them as a separate dataset rather than removing them entirely. The Split Data module is also useful when the user wants to isolate individual groups of spectra from the dataset, such as spectra that have been flagged with unique cluster identifications by a clustering algorithm or splitting off an individual spectrum.

Stratified Folds

The Stratified Folds module is used to divide a dataset into N number of groups, also called folds, that can be used to cross validate, train, and test a regression model. This is required to ensure that regression models are generalizable to new data and to determine the accuracy of a model. Refer to the sections below for more detail on cross validation, regression, and the importance of test and train set design. To ensure that the folds have similar distributions of the variable of interest, the data are stratified—that is, they are sorted on a user-specified variable of interest, and then each unique value of that variable is assigned to folds sequentially. Thus, if the data are sorted on SiO₂ content, for example, and three folds are desired, the lowest SiO₂ value will be assigned to fold 1, the next lowest to fold 2, and so on. Optionally, the user may also specify a secondary tiebreaker variable. This is useful in cases where there are multiple samples with identical values in the primary stratification variable but different properties overall (for example, if there are many different samples with a composition of zero in the primary variable). Specifying a tiebreaker can avoid lumping all of these samples together in one fold. The user can select which fold to use as a test set. This module creates two new datasets, one for the training folds and one for the test fold.

Identify Endmembers

The Identify Endmembers module allows users to identify certain rows as endmembers, which typically correspond to pure spectral signatures. In spectral analysis, users may want to understand how much of a target that has a mixed composition, such as a pixel that includes both rocks and soils, is represented by the individual components, rock and soil in this example, also known as endmembers. In algorithms such as pixel purity index, data are projected in such a way that the most spectrally extreme, and thus spectrally “pure,” pixels are detected and identified as endmembers. Once endmembers are identified, they can be used in unmixing classification methods that estimate the relative contribution of endmembers to a particular spectrum. PyHAT includes several endmember identification algorithms (described below and in table 1): pixel purity index (PPI), N-FINDR, automatic target generation process (ATGP), fast iterative PPI (FIPPI), and sequential maximum angle convex cone (SMACC). To run these algorithms in PyHAT, the user must define the number of endmembers. It is important to note that many endmember identification algorithms at their core are functionally similar yet may nevertheless yield different results; this is because of differences in preprocessing steps, such as the randomization of initial conditions performed in some algorithms (Chang and others, 2016).

An example of the results of each PyHAT endmember identification algorithm is shown in figure 9. This example uses the Salinas Scene dataset (available at https://paperswithcode.com/dataset/salinas, accessed June 26, 2025), which is a high-resolution Earth remote sensing dataset of agricultural lands in the Salinas River valley, California. Although the dataset contains ground truth information, where the endmember compositions are known (that is, each pixel contains a numerical ground truth variable “gt” that refers to a specific type of vegetation) we ignore this information in the example and allow an endmember identification algorithm to find a select number of spectra that represent the endmembers of the dataset. Some algorithms overlap in the identified spectra and others result in unique identifications of endmembers. This figure is purely for illustration purposes because only the raw spectra were used to compute endmembers; preprocessing steps, such as dimensionality reduction, were not performed, likely resulting in poorer performance of certain algorithms. The PPI algorithm, for example, identified two rather similar spectra as endmembers.

Spectral Parameters

Particularly when working with reflectance spectra, it is useful to calculate “spectral parameters” that summarize some aspect of a spectrum into a single value that can then be plotted or mapped to aid in interpretation. Reflectance of a single band, reflectance ratios, slopes, band depths, and band asymmetry are all common spectral parameters. PyHAT implements many predefined spectral parameters that are commonly used with the CRISM instrument on Mars and the M³ instrument on the Moon. CRISM parameters are based on those defined by Viviano and others (2014) and the Interactive Data Language (IDL) code released for the CRISM Analysis Toolkit (CAT; Morgan and others, 2017). Where the two CRISM references differ, we defer to the formulation in the CAT. M³ parameters were provided by Lisa Gaddis (M³ team, written comm., 2023). For convenience, table 2 lists the parameters for M³ and table 3 lists the parameters for CRISM. For details of the parameter calculation, refer to the documentation of the individual spectral parameter functions in the PyHAT code. When calculated, the parameters are added as a new column in the spectral data frame, with a top-level label of “parameter” and a second-level label of the parameter name.

Resample

In workflows involving multiple different instruments that have different spectral resolutions and overlapping spectral overage, it is useful to resample data from one instrument to the wavelength spacings of another.

The Resample module has a simple interface with two drop-down menus. The first allows the user to choose which data to resample from among the available datasets. The second drop-down menu allows the user to specify which dataset will be used as the reference dataset. When this module runs, it uses the SciPy function scipy.interpolate.interp1d (https://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.interp1d.html#scipy.interpinter.interp1d, accessed July 9, 2025) to linearly interpolate the data being resampled onto the wavelengths of the reference dataset. Note that this module will only interpolate between existing data points and will not extrapolate spectra beyond the maximum or minimum wavelength of the original dataset. In situations where the dataset extends beyond the wavelength range of the other, both will be truncated to where they overlap. For example, if the data to be resampled had a wavelength range of 400 to 800 nanometers (nm) and the reference data had a wavelength range of 500 to 900 nm, resampling would result in both datasets having an extent of 500 to 800 nm. When truncation occurs, the program will warn the user by listing the wavelengths of the columns that will be removed in the console window.

Remove Baseline

Spectral data typically consist of three components: baseline, signal, and noise (Shen and others, 2018). Baseline removal, also known as continuum removal, is a common preprocessing step in spectroscopy that seeks to remove signal that contains little to no diagnostic information about the spectral features of interest (Giguere and others, 2015). For example, when working with nongated LIBS data, such as ChemCam data, the spectra contain broad continuum emission owing to bremsstrahlung and ion-electron recombination processes (Giguere and others, 2015). Additionally, baseline removal can be important when working with visible and near-infrared reflectance data, like those from CRISM, especially for comparisons across platforms (Giguere and others, 2015). In the Remove Baseline module, PyHAT includes multiple algorithms that can be used for baseline removal (table 4), which are discussed in detail in the following sections. The code for the algorithms was originally written by C.J. Carey and was provided by T. Boucher (both at the University of Massachusetts at the time). It has been modified to work within PyHAT.

It is challenging to optimize the baseline removal routine because there is no background-free spectrum to compare against. Users should first determine whether continuum removal is appropriate and necessary for their application. As with other algorithms, users should try a variety of approaches and parameters, plot the results, and visually inspect the identified baselines to qualitatively determine which method reliably yields appropriate results. Figure 11 shows an example comparing several different baseline removal algorithms on a LIBS spectrum.

Apply Mask

With the Apply Mask module, the user can apply a mask to a dataset to exclude specified parts of the spectrum from analysis. The range of wavelengths to mask are specified in a simple CSV file format (fig. 12). The first column of the file contains the name of the feature to be masked. This column is primarily to make the CSV file more human readable and is not used by the masking algorithm. The second and third columns specify the minimum and maximum wavelengths to mask, respectively. To mask multiple ranges, list them in separate rows of the CSV file. The module reads the CSV file and then steps through each wavelength range, changing the top-level column name for wavelengths from “wvl” (or the user-specified spectral data column label) to “masked.” The data are not removed from the dataset; their labels are simply changed so that the masked data are not used in subsequent calculations. The masking is inclusive (that is, it uses the ≥ and ≤ operators) so if the wavelength range is 200 to 210 nm and a spectral channel has a wavelength of precisely 200 nm, it would be masked. It is generally best practice to choose minimum and maximum wavelength values that fall in between the spectral channels of the data. This ensures that the columns that were masked will not change if a computer rounding error changes the precise value of the wavelength labels for the columns.

Peak Areas

The Peak Areas module implements a version of the method used by Clegg and others (2020) for emission spectroscopy applications to consolidate the data from the adjacent spectral channels on an emission line into a single bin. It is a fast and simple alternative to more precisely fitting individual emission lines in a spectrum and can be useful for collecting data from broad and (or) weak features for use with multivariate methods. It can also make spectra more robust to small wavelength shifts in the peaks (for example, caused by varying instrument temperature) and reduces the size of spectra to enable faster calculations (that is, several spectral bands across a single peak are reduced to a single value).

The method relies on a set of wavelengths corresponding to local minima and maxima, typically derived from an average spectrum (fig. 13A). For each individual spectrum (fig. 13B, colored lines), the data values between each pair of local minima are summed and assigned to the wavelength value of the local maximum within that range (fig. 13B, brackets). Note that because the minima and maxima are based on the average spectrum, they may not correspond to exact minima and maxima for a given individual spectrum (for example, the red spectrum in fig. 13B has peaks that are one pixel away from the peak locations defined by the average spectrum). The average spectrum is used to predefine the bins so that the data are summed into the same wavelength channels for each spectrum.

The default option is for the module to calculate the average spectrum and define the local minima and maxima, but the user can choose to load the values from a file. If the user chooses to calculate the values, they will be written out to a file named peaks_mins.csv. This file can be used to apply the same peak area binning to multiple different datasets. The user can also edit this file to have greater control over the binning ranges.

Multiply by Vector

The relatively simple Multiply by Vector module allows the user to multiply spectra by a user-provided one-dimensional array of values (vector) from a file. This is useful for applying corrections to the spectra, such as an instrument response function. The vector file has a simple format; the first column contains the wavelengths of spectral channels, and the second column contains the value by which to multiply that spectral channel.

Normalize

The Normalize module divides one or more user-defined spectral ranges by their respective totals, such that for each spectrum, the signal in each range sums to 1. This type of normalization can be used to correct for differing signal levels from different spectrometers and is commonly used in LIBS preprocessing to mitigate fluctuations in plasma emission intensity (Wiens and others, 2013).

Standardize

The Standardize module allows a user to standardize their data such that the spectra in a dataset have a mean of zero and a standard deviation of one (fig. 14). The mean and standard deviation correction factors can be derived from a different dataset if desired (for example, to standardize a test set based on a training set). Standardization is important if the input data have different units and (or) if the user wants all channels to have equal weight. This can reduce the influence of strong spectral features, for example, deep absorption features or strong emission lines, and help to emphasize the variability of weaker spectral features. The user should exercise caution when applying this step because it can also amplify noise—noise can have similar variability as weak spectral features. When this module is run, a new dataset containing the variance and mean vectors is created. To reverse the standardization, the standardized spectra can be multiplied by the square root of the variance vector and added to the mean vector.

Dimensionality Reduction

With the Dimensionality Reduction module, data dimensionality can be reduced using various algorithms, which are discussed in the following sections and in table 5. Dimensionality reduction consolidates information from numerous variables, like wavelengths, into a few components. This allows for faster processing and more interpretable results. It can also help with endmember identification and spectral unmixing (Ghamisi and others, 2017). Within dimensionality reduction, there are supervised and unsupervised methods (Shahdoosti and Mirzapour, 2017). In unsupervised methods, such as principal component analysis (PCA) or independent component analysis (ICA), samples do not need to be labeled. Supervised methods, like linear discriminant analysis (LDA), on the other hand, use labeled samples and select dimensions that maximize the distance between classes (Ghamisi and others, 2017).

Dimensionality Reduction Algorithms

Table 5.    

Description of dimensionality reduction algorithms available in PyHAT.

[#, number; 2D, two dimensional; 3D, three dimensional]

Table 5.    Description of dimensionality reduction algorithms available in PyHAT.

Algorithm	Description	Parameters
Principal component analysis (PCA)	Algorithm that treats spectra as points in a high-dimensionality cloud and identifies projection onto a small number of components that capture most variation.	# of components
Independent component analysis (ICA)	Works similar to PCA, but seeks to identify components that are independent from one another, rather than maximizing variance.	# of components
t-distributed stochastic neighbor embedding (t-SNE)	Nonlinear, projects components in 2D or 3D map, results are easy to visualize but the meaning of plot axes is difficult to interpret.	# of components, perplexity (# of nearest neighbors), learning rate, # of iterations, # of iterations with no progress (determines when to stop)
Locally linear embedding (LLE)	Uses a series of local PCAs to find a low-dimensionality nonlinear projection that preserves local distances.	# of neighbors, # of components, regularization (to simplify model)
Nonnegative matrix factorization (NNMF)	Dimensionality reduction and unmixing, no a priori information needed, no assumption of number of endmembers (or pure samples), able to handle highly mixed cases	# of components, whether to add a constant
Linear discriminant analysis (LDA)	Efficient supervised dimensionality reduction and classification. Seeks to maximize ratio of between-class variability to within-class variability.	# of components, class column (labels for supervised method)
Minimum noise fraction (MNF)	Dimensionality reduction and noise removal. Uses two PCAs to rank components by signal-to-noise ratio.	# of components
Local Fisher discriminant analysis (LFDA)	Dimensionality reduction and classification. Similar to LDA, but tries to preserve data structure while maximizing class separability. Does not require unimodal Gaussian data distribution.	# of dimensions, class column, type of metric, # of neighbors

Principal Component Analysis (PCA)

Principal component analysis (PCA) is a linear dimension reduction technique that uses singular value decomposition (SVD). It can be visualized by treating each spectrum as a point in a high-dimensionality space, where each wavelength is a dimension. PCA identifies the axis through the cloud of points in that space along which there is the most variability (fig. 15 illustrates the two-dimensional case). This axis is the first principal component. Subsequent orthogonal axes explain progressively less of the total variation in the dataset. The result is that information from many variables (wavelengths) can be conveyed in a much smaller number of variables (principal components). The values of each spectrum when projected onto the principal components are called the scores, and the vector that is used to project the spectrum and derive the score is called the loading.

The PCA implementation in PyHAT adapts the scikit-learn LAPACK (linear algebra package) method using a full SVD or a randomized truncated SVD depending on the shape of input data and the user-specified number of components. For many datasets in spectroscopy, the first few components explain almost all variability in the data. For more information regarding PCA, refer to the scikit-learn user guide, available at https://scikit-learn.org/stable/modules/generated/sklearn.decomposition.PCA.html.

Independent Component Analysis (ICA)

Independent component analysis (ICA) is a dimensionality reduction method intended to separate d individual sources from mixed signatures. The basic ICA form is X = AS, where X is an n×m matrix such that n is the number of measured signals (spectra) and m is the number of variables (spectral channels), A is an n×d matrix, and S is a d×m matrix of source signals. PyHAT uses two ICA methods. The first is FastICA, as implemented in scikit-learn (Hyvärinen and Oja, 2000). FastICA uses weight vectors to determine one independent component at a time, which is more efficient than other approaches (Wang and others, 2008). Within PyHAT, the user defines the number of components for FastICA. For more information on FastICA, refer to the scikit-learn user guide, available at https://scikit-learn.org/stable/modules/generated/sklearn.decomposition.FastICA.html.

The other ICA method in PyHAT is joint approximate diagonalization of eigenmatrices (JADE; Cardoso and Souloumiac, 1993). PyHAT uses an implementation of JADE translated from Cardoso’s MATLAB version into Python by Gabriel Beckers (https://github.com/gbeckers/jadeR). JADE ICA does not need gradient searches and so is able to converge easily but can be slow for larger datasets (Wang and others, 2008). JADE ICA performs well in extracting the signal from individual chemical elements from LIBS spectra (Forni and others, 2013). As with PCA and FastICA, the number of components is defined by the user.

t-Distributed Stochastic Neighbor Embedding (t-SNE)

Whereas PCA performs linear dimension reduction, t-distributed stochastic neighbor embedding (t-SNE; van der Maaten and Hinton, 2008) is a nonlinear method for representing high-dimensionality data in a two- or three-dimensional map. Whereas the standard SNE uses a Gaussian distribution, t-SNE uses the Student’s t-distribution when it is computationally faster (van der Maaten and Hinton, 2008). The outputs of t-SNE are easier to visualize than other dimension reduction methods because it allows for the visualization of data points that are very close together, and thus similar, as well as those that are very dissimilar. PyHAT uses the scikit-learn implementation of t-SNE, available at https://scikit-learn.org/stable/modules/generated/sklearn.manifold.TSNE.html. The user needs to define several parameters: the number of components, perplexity, learning rate, number of iterations, and number with no progress. The number of components specifies the dimensionality of the embedding and should not exceed three. Perplexity controls the effective number of nearest neighbors used and generally needs to be increased for larger datasets. The learning rate controls the rate of convergence and is typically in the range of 10 to 1,000. Number of iterations sets the maximum number of iterations, and number with no progress allows the algorithm to abort before reaching the maximum number if no progress is being made. Number with no progress is rounded to the nearest multiple of 50.

Locally Linear Embedding (LLE)

Locally linear embedding (Roweis and Saul, 2000) is a commonly used nonlinear dimensionality reduction algorithm owing to its ability to find global optima and its computational efficiency (Qing and others, 2013). PyHAT uses the scikit-learn implementation of LLE, available at https://scikit-learn.org/stable/modules/generated/sklearn.manifold.LocallyLinearEmbedding.html. The user-defined parameters are the number of neighbors, number of components (the dimensionality of the manifold to be used), and regularization. It should be noted that the number of neighbors can have a strong influence on the results and can be challenging to optimize (Qing and others, 2013).

Nonnegative Matrix Factorization (NNMF)

Nonnegative matrix factorization (NNMF) is an additive linear model that is used for dimensionality reduction and unmixing (Gillis and Plemmons, 2013; Wen and others, 2013; Liu and others, 2016a). NNMF does not require a priori knowledge of what is present in the input dataset, and, by enforcing nonnegativity constraints, it has an advantage over other dimensionality reduction techniques in that its results are more physically meaningful (Huang and Zhang, 2008; Liu and others, 2016a). This algorithm can estimate the relative contribution of endmembers to the observed spectra even in scenarios where a spectrum is composed of multiple other spectra (highly mixed), without assumptions of the presence of pure spectra (Ceamanos and others, 2011; Liu and others, 2016a). This is an advantageous attribute for analyzing orbital hyperspectral data when spatial resolution is limited and pure pixels containing a single mineral species are rare. Note, however, that especially in visible and near-infrared spectra, there are complex nonlinear effects and, in general, the results of unmixing such as NNMF should not be interpreted as abundance.

The algorithm deconstructs the input data nonnegative matrix V into two nonnegative matrices H (endmembers) and W (weights, or abundances), whose product approximates V, with W and H remaining nonnegative, and the columns in W summing to unity (Huang and Zhang 2008; Casalino and Gillis, 2017). This algorithm has been used for dimensionality reduction for terrestrial land cover studies (Huang and Zhang, 2008) and for unmixing endmembers for the classification of vegetation and land cover types (Wen and others, 2013, 2016; Zhu and Honeine, 2016), minerals (Liu and others, 2016a), and dusty features on Mars (Ceamanos and others, 2011). Within PyHAT, the user needs to define the number of components and choose whether to add a constant to ensure results are nonnegative.

Linear Discriminant Analysis (LDA)

Linear discriminant analysis (LDA) is a computationally efficient and easily interpretable supervised dimensionality-reduction and classification method (Shahdoosti and Mirzapour, 2017). LDA is also faster than methods like support vector machines (SVM) that require parameter optimization (Zhang and others, 2015a). The algorithm reprojects data into a K−1 dimensional space, where K is the number of classes being separated. It does this in a way that maximizes the ratio of between-class variability to within-class variability, making the resulting dimensions ideal for use in classification.

LDA has been used for a wide variety of applications, including distinguishing subtly different growth stages of insect puparia (Voss and others, 2017) and identifying mung bean varieties (Xie and He, 2018). This algorithm is used in a variety of applications, from land cover studies (Bratsch and others, 2016; Calvino-Cancela and Martin-Herrero, 2016; Shahdoosti and Mirzapour, 2017) to determining chicken embryo gender (Gohler and others, 2017). In addition to these terrestrial applications, LDA has been used for planetary studies, including classification and quantification of analytes in a simulated Martian atmosphere (Benkstein and others, 2014) and mineral classification using CRISM data (Bue, 2014). In PyHAT, the user needs to define the number of components and identify the class column because LDA is a supervised method.

Minimum Noise Fraction (MNF)

Minimum noise fraction (MNF) is commonly used for dimensionality reduction (Das and Seshasai, 2015), especially before endmember identification (discussed above). MNF is also used for noise reduction (Bjorgan and Randeberg, 2015; Hamzeh and others, 2016; Houborg and others, 2016). An MNF transformation consists of two consecutive principal component analyses (PCA), the result of which is a set of MNF components ranked by signal-to-noise ratio (Kale and others, 2017). Removing the noisiest MNF components and reconstructing the image allows for the removal of noise with the retention of information (Parente and others, 2009). This algorithm has been used in a planetary science context to remove noise from Mars lander and rover multispectral observations (Bell and others, 2002; Parente and others, 2009) and orbital hyperspectral data (Douté and others, 2007). To run MNF in PyHAT, the number of components must be defined.

Local Fisher Discriminant Analysis (LFDA)

The local Fisher discriminant analysis (LFDA) algorithm is like the LDA but uses a locality preserving projection (LPP) to preserve data structure while maximizing class separability (Bue, 2014). LPP is a dimensionality reduction technique that finds a linear map that preserves input data structure of neighboring sample points in input data dimensions (Ye and others, 2017) using k-nearest neighbor (Luo and others, 2015). As opposed to other popular dimensionality reduction techniques, such as PCA and LDA, LFDA does not require a unimodal Gaussian data distribution (Luo and others, 2015; Ye and others, 2017). This is particularly advantageous for spectroscopic data analysis because these data are commonly multimodal (Luo and others, 2015). Also, LFDA allows for dimensionality reduction to greater than K−1 dimensions, where K is the number of classes (Bue, 2014). In PyHAT, the user must define the number of dimensions (components), class column, type of metric (plain, orthogonal, or weighted), and number of neighbors.

Unmixing

The Unmixing module within PyHAT provides several algorithms (table 6) that help identify the contribution of endmember spectra to a given spectrum where multiple endmembers may be present; for example, multiple materials within a single pixel of remote sensing data. These algorithms tend to involve two constraints: (1) nonnegativity, which guarantees a zero or positive contribution of each endmember and (2) sum to unity, which guarantees the determined contribution of endmembers sums to one. Unmixing algorithms can also be linear or nonlinear, allowing the abundance of each endmember to affect the composite spectra linearly or nonlinearly (see Schmidt and others, 2014, Goudge and others, 2015, and Liu and others, 2016b, 2018, for examples).

Spectral Derivative

The Spectral Derivative module calculates a simple approximation of a spectral derivative by finding the difference between each value in the spectrum and the adjacent value, using the pandas diff function. Spectral derivatives are useful for identifying the location of peaks and isolating features from a smoothly varying baseline (fig. 16). When run, this module creates a new dataset containing the derivative spectra. The new dataset has the same name as the original dataset, with the string “ - Derivative” appended.

Wavelength Shift

The Wavelength Shift module allows the user to generate sets of artificially shifted spectra. This can be used to evaluate the robustness of certain analytical steps, such as in cases where the wavelength calibration of data may vary slightly (for example, owing to temperature fluctuations in the instrument). In the text field of the GUI, the user specifies one or more shift amounts separated by commas. When run, this module will add the shift amount to the wavelength values in the column names, then resample the original spectra onto the new wavelengths. This will be done for each of the shifts specified in the user interface, and the resulting data are concatenated and stored in place of the original data. Note that the shifting process results in sections of no data at the end(s) of the spectra, which are automatically cut off.

Calibration Transfer

Scientists commonly need to work with data collected by more than one instrument or under varying experimental conditions (for example, training data collected in a laboratory on Earth being used to interpret data collected by an instrument on a Mars rover). In these cases, simply resampling the wavelengths of the data may not be sufficient—the nature of the spectral response may be different even after instrument response corrections have been applied. However, it can be expensive or impossible to collect training data on both instruments and at all possible conditions. It is therefore desirable to determine a mathematical transformation that can correct spectra collected under one set of conditions so that they can be used with data collected under a different set of conditions. This process is called calibration transfer, and PyHAT provides several different algorithms that can be used for this purpose. The original code for these algorithms was kindly provided by T. Boucher (University of Massachusetts, written commun.). The GUI includes a Calibration Transfer CV (cross validation) module and a Calibration Transfer module.

Calibration transfer relies on a set of spectra of the same targets collected under different conditions. The greater the number of common spectra that are available, the better the results are likely to be. For example, the ChemCam and SuperCam instruments carry calibration targets with them on the rover. Calibration data collected with the flight instrument on Mars can be compared with spectra of identical targets observed under laboratory conditions on Earth and used to derive a calibration transfer model. Calibration transfer can then be used to transform laboratory data to be more Mars-like for calibration with the Mars spectra.

Calibration Transfer CV Module

Typically, the optimal settings for calibration transfer for a given dataset and algorithm are not known a priori. Instead, it is necessary to evaluate a wide range of settings to determine which works best. The Calibration Transfer CV module allows the user to do this by running cross validation on the data. It iteratively removes each spectrum from the data, generates a calibration transfer model using the remaining spectra, and evaluates the performance by calculating the root mean squared error (RMSE). A low RMSE indicates that the transformed spectra are more like the spectra collected under different conditions.

The interface for this module is shown in figure 17. At the top of the interface, the user can specify the two datasets to use. These datasets should contain spectra of the same targets, observed under two different conditions. Dataset A will be transformed to be as similar as possible to dataset B. The drop-down menus on the top right specify the metadata that will be used to match spectra of a given target in dataset A to the spectra of the same target in dataset B. The module will automatically remove any rows that do not have a corresponding entry in the other dataset. It will also automatically average multiple spectra of the same target together such that the calibration transfer algorithms are working with data frames containing a single average spectrum for each target in both datasets. When the module is running, the console window will inform the user as these preparation steps are being applied.

Beneath the drop-down menus, there is a check box that allows the user to choose whether to keep the transformed spectra as new datasets, or whether just the cross-validation results should be retained. Keeping the spectra can be informative, as they can be plotted to see how a given algorithm and set of parameters transform the spectra, but this option significantly slows down the cross-validation calculation. For better performance when evaluating many algorithms and parameter settings, it is recommended to uncheck this box.

Below the “keep transformed spectra” option there are several check boxes, one for each algorithm to be evaluated in the cross validation (see below for descriptions of these algorithms). Choosing each box will reveal the interface for that algorithm (if there are any parameters). Text fields allow the user to enter multiple comma-separated parameter values and list widgets allow the selection of multiple options so that the user can explore the full parameter space for the selected algorithms. Typically, these fields are populated with the default parameters, but users are encouraged to vary the parameters to suit the calibration transfer task at hand. It is commonly useful to perform a wide-ranging, coarsely gridded search through the parameter space as a first pass, followed by more densely spaced parameters around the settings that give the best initial results.

K-Means

K-means clustering is one of the most used clustering methods. It seeks to divide the data into a user-specified number of clusters (k) of equal variance. It works by starting with k cluster centroids, randomly distributed in the same space as the spectra to be classified (fig. 18). Each point is assigned to the nearest cluster centroid, and then a new centroid is calculated for each cluster, based on the spectra assigned to that cluster. This process is repeated until the centroid no longer changes with each iteration, within the user-specified tolerance value. It is recommended (though not required) that dimensionality reduction be applied to the spectra and k-means be run on the reduced dimensionality values rather than the full spectra. This significantly speeds up the calculation and can be used to ensure that clustering is based on meaningful values rather than noise (for example, by using the first several PCA components, which explain most of the variance in the data). K-means is susceptible to finding local minima, so typically the algorithm is run numerous times and the result with the minimum sum of squared distances in each cluster is used. K-means assumes that the clusters in the data are convex and isotropic, so it can give poor results when the actual cluster in the data are elongated or have irregular shapes.

Spectral

Spectral clustering is an alternative to k-means clustering that performs well when clusters have complex shapes that are poorly described by the cluster center and spread. PyHAT uses the scikit-learn implementation of spectral clustering, available at https://scikit-learn.org/stable/modules/generated/sklearn.cluster.SpectralClustering.html. Like k-means, the user can specify the number of clusters and the number of runs of the algorithm to use. Spectral clustering works by generating an affinity matrix between the data points and embedding it in a low-dimensional space where simple clustering like k-means is more effective. The user can select from several kernels to be used to generate the affinity matrix: radial basis function, nearest neighbors, sigmoid, polynomial, linear, and cosine. These kernels are described in the scikit-learn documentation, available at https://scikit-learn.org/stable/modules/metrics.html. In general, most kernels define a symmetric function that determines the similarity between points in such a way that closer points are given greater weight than more distant points. The parameters associated with the kernels adjust this falloff in sensitivity. In practice, the user can experiment with the parameters to get acceptable results.

Overfitting and Accuracy

Overfitting occurs when model complexity is increased (or other model parameters are tuned) such that the model performs well on the training set but fails to predict new data accurately. This is analogous to using a high-order polynomial to fit noisy data (fig. 20); the polynomial may fit the training data perfectly (black points in fig. 20), but it does not accurately capture the underlying trend and may give inaccurate results for new x values (red point in fig. 20).

Similarly, calculating the root mean squared error (RMSE) of the model’s predictions of the training data (also known as RMSE of calibration [RMSEC]) can give an incorrectly optimistic estimate of the model’s accuracy. In the case of the overfitted polynomial, the RMSEC would be zero, suggesting that the model is perfect, when the performance on novel data shows that it is not.

To avoid overfitting, cross validation (described below) should be used to tune the model parameters, and an independent test set should be used to assess model accuracy. The RMSE of the test set is also known as the RMSE of prediction (RMSEP). The test set is data that the model has never seen before—it should not be used to tune the model parameters, only to assess the performance of the tuned model.

Cross Validation

Cross validation is an essential step in optimizing a multivariate regression model and ensuring that the model is not overfitted to the training data. In cross validation, a subset of the training data is held out and treated as unknown. A model is then trained on the remaining data and used to predict the held-out data. This process is repeated iteratively so that every piece of training data is held out as an unknown once. The result is an estimate of how a model trained on the full training set will perform on novel data—the root mean squared error of cross validation (RMSECV). By repeating this process for different parameter settings for a model, the optimal settings (those that minimize RMSECV) can be determined. RMSECVs can also be compared between models to determine which gives the best results for the task at hand.

PyHAT implements stratified k-fold cross validation. The Stratified Folds module, described above, is used to split the data into training and test folds with similar distributions of the dependent variable of interest. One fold is completely held out as a test set to evaluate the performance of the final optimized model (RMSEP). The remaining folds are used for cross validation to tune model parameters and compare different regression models.

In the interface for regression cross validation, the user can specify which of the regression algorithms to use. For each algorithm, several parameters can be specified. By entering comma separated values in text fields or selecting multiple entries in list boxes, the user can determine all the permutations to be evaluated. Note that for some models and parameter settings, cross validation can be time consuming. The console window will display updates on the progress of the cross validation.

Regression Train

If the user already knows the optimal parameters for training a regression model, they can bypass cross validation and train a single regression model. The same algorithms are available as for cross validation; the primary difference is that only one algorithm can be chosen at a time, and each parameter accepts only a single value rather than a list of values. As with cross validation, the range of values for the dependent variable can be adjusted to permit a specialized regression model to be trained on a restricted range of input data. Note that the Regression Train module trains a regression model but does not perform predictions.

Regression Predict

Once the regression model has been trained by either using cross validation or regression training, the user can apply the model using the Regression Predict module. This module allows the user to choose from a list of trained regression models and use the selected model to predict the results for one or more datasets. The predicted results are added to the selected datasets as a column, with a top-level label of “predict” and a second-level label specifying the regression model used.

Blend Submodel Predictions

The Blend Submodel Predictions module allows the user to implement the submodel regression approach described by Anderson and others (2017). This approach is useful when the unknown data to be predicted may exhibit a wide range of values for the dependent variable, such that a single regression model would have trouble accurately handling all possible cases. Instead, by training multiple regression models, each on a restricted but overlapping range of values, the results can be blended to yield a more accurate prediction across all values. This method uses a reference model trained on the full range of dependent variable values (for example, 0–100 weight percent for a regression model predicting the SiO₂ content in a target analyzed using LIBS) as an initial guess of the true value. This initial guess is then used to determine which more specialized model (submodel) should be used to get a more accurate result. If the reference model predicts a value that falls in between a specified overlap range between two submodels, the submodel predictions are blended using a linear weighted average. This overlap range is a free parameter that can be optimized; it does not necessarily have to correspond to the overlap in training data ranges for the submodels.

This module does not perform predictions; it combines predictions that have already been run using the Regression Predict module. When the user selects a dataset, the prediction columns are read from that dataset and used to populate the drop-down lists of predictions for the reference model and submodels. The interface begins with just a reference model and a low and high submodel. If the user wants to use predictions from additional submodels, they can click the “Add submodel prediction” button, which will insert a new row between the low and high submodels. For each submodel, the user must specify the range over which it should be used. Thus, for a low model trained on 0–50 weight percent SiO₂, you might choose 50 as the maximum value for that model. For a middle submodel trained on 40–60 weight percent SiO₂, those ranges should be entered as the minimum and maximum values. The ranges selected for the submodels must overlap for this method to work. If they do not overlap, it is possible for the reference model to predict values that do not correspond to any submodel.

In practice, most users will want to optimize their blending ranges to get the best possible results. This can be done by predicting the training data and then selecting it in the Blend Submodel Predictions module. Load the desired submodel predictions and then select the “Optimize” check box. This will bring up an additional drop-down list, where the dependent variable of interest can be selected. When the module runs, it will start with the minimum and maximum blending ranges provided in the user interface for each submodel and adjust them to minimize the error in predicting the variable of interest. The progress of the optimization is printed to the console window. Once the optimal ranges have been determined using the training data, they can be used for other datasets. Note that optimization should not be applied on test set data; doing this will tune the blended results to give the best possible test set results, and thus the model settings will no longer be independent of the test set.

Save or Restore Trained Model

If desired, trained regression models can be saved to a user-specified binary file. This file can be restored using the Restore Trained Model module. When a file is selected in this module, information about the model in that file is displayed in the console window.

Regression Algorithms

The following sections provide a brief overview of the available regression algorithms in PyHAT. PyHAT relies heavily on the scikit-learn library for these algorithms. It is beyond the scope of this guide to describe the details of each algorithm. Instead, the focus is on practical considerations when using the algorithms as implemented in PyHAT. Additional information and documentation on these methods can be found in the scikit-learn documentation, available at https://scikit-learn.org/stable/index.html.

In most cases, trained linear regression models result in a vector of regression coefficients (w) that, when multiplied by a spectrum (x), yield a predicted y value:

It can be informative to plot the regression coefficients as a function of wavelength to see which parts of the spectrum have the strongest influence on the predicted results. Thus, when a regression model is trained in PyHAT, the regression coefficients are added to a new dataset labeled “Model Coefficients.” Rows from this dataset can be plotted using the Plot Spectra module, or the data can be saved to a CSV file and analyzed using another program.

Table 8 provides a high-level summary of the available regression algorithms and their properties.

Partial Least Squares (PLS)

Partial least squares (PLS) is a linear regression method closely related to PCA that works well in situations where there is a high degree of correlation between the independent variables (x) and where there are more independent variables than there are observations. Both situations are common in spectroscopy.

Like PCA, PLS projects high-dimensionality data into a lower dimensionality space (fig. 21). The difference is that in PLS both the x and y variables are projected into a new space, and the axes through the cloud of x data points are chosen such that they correspond to the greatest variation in the y space.

There are several different algorithms for PLS; PyHAT uses the nonlinear iterative partial least squares (NIPALS) algorithm (Wold, 1975) as implemented in the scikit-learn PLSRegression function.

The key parameter for any PLS model is the number of components or latent variables. A model with too few components will not adequately capture the complexity of the data and will not result in accurate predictions. A model with too many components will be overfitted, performing well on the training set (low RMSEC) but failing to predict new data accurately (high RMSEP). Thus, it is important to use cross validation to tune the number of PLS components appropriately.

Explore and Calculate Local RMSEP Parameters

The root mean squared error of prediction (RMSEP) provides a single summary statistic to indicate the prediction performance of a regression model on an independent test set. However, it is common for models to have better performance in certain composition or parameter ranges and worse performance in others. For example, if the training dataset contains mostly low compositions, the model may perform better close to zero but poorly at high concentrations. In this case, the overall RMSEP may be overly pessimistic for low predictions and overly optimistic for high predictions. An alternative, developed for expressing the accuracy of ChemCam and SuperCam quantification, is to report a local RMSEP (Clegg and others, 2017; Anderson and others, 2022). This is achieved by considering a series of simulated predictions at different values and calculating the RMSEP of only those test set predictions within a certain window around the simulated prediction in question. This window can be defined using an absolute composition range (use test set predictions within 1 weight percent of the simulated prediction), a minimum number of neighbors (use the 40 nearest test set predictions), or a combination of the two (use all test set predictions within 1 weight percent but expand that range to include at least 40 predictions). The result is a series of local RMSEP values corresponding to the simulated prediction values. This series of values commonly has abrupt steps, particularly in areas where the distribution of test set predictions is sparse and the same test set predictions fall into the window used to calculate the local RMSEP for multiple simulated predictions. To give a more realistic estimate of the RMSEP at a given point, the local RMSEP values can be smoothed.

PyHAT provides two modules related to local RMSEP. The first is the Explore Local RMSEP Parameters module, which allows the user to provide a series of values for the minimum number of neighbors, the window size, the amount of smoothing, and whether to extrapolate the local RMSEP values beyond the range of the test set. When run, multiple plots of the local RMSEP as a function of predicted composition will be generated, allowing the user to visually inspect the results and identify those results that give a reasonable-looking local RMSEP curve. Once the parameters have been identified, the Calculate Local RMSEP module can be set to those parameters to calculate the local RMSEPs corresponding to a set of actual predictions.

Plot

The Plot module allows the user to create scatterplots of one column in a dataset against another column in the same dataset (individual spectra can be plotted using the Plot Spectra module described below). The user first must select a dataset to plot. The module will then read all the column labels from the dataset and use these labels to populate the list of variables that can be plotted. Columns containing the spectra, masked spectral channels, and peak areas are excluded to keep the number of variables to choose from manageable.

The user can choose to create a new figure or, if the plot module has been run previously in the same workflow, choose an existing figure to plot additional data on the same axes. If a new figure is being created, the user may specify the figure name, plot title, axis titles, and axis ranges. When creating a new figure, if the user selects an x or y variable from the list, the corresponding axis title will be set to that variable name, and the range will be set to the minimum and maximum values of that variable. The user may then manually adjust these fields as needed.

If an existing figure is selected, the titles and axis limits are inherited from the previous instance of the Plot module and cannot be changed. Thus, when plotting multiple datasets, the user should set the plot axis ranges in the first instance of the plot module so that the axes span the full range of all data that will be plotted. If this is not done, the subsequent datasets may plot outside the bounds of the original axes.

The user should enter a name for the dataset being plotted into the legend field so that the plot legend is informative. This is especially important when plotting multiple datasets on the same axes.

The user may also adjust settings related to the plot markers. The first option is whether to color the markers based on the value of another column from the dataset. Color coding is a way to visualize a third variable on a two-dimensional scatterplot. If color coding is selected, the options for marker color, opacity, and connecting line style are turned off. The shape of the marker may still be adjusted. The resulting plot will include a color bar on the right labeled with the name of the variable used for color coding. Plotting multiple datasets with color coding on the same figure is not recommended; this will result in multiple color bars on the right of the figure and a confusing plot. If the user selects “None” for color coding, then it is possible to choose the marker color, opacity, and the style of the lines connecting the markers (if any). The marker shape may also still be adjusted. The color of the connecting line will be the same as that of the markers.

There is also a check box labeled “one to one.” Checking this box indicates to the module that the data being plotted is a one-to-one plot, of the sort used to evaluate a regression model (fig. 22). This will cause the module to draw a thick black line from the origin marking where the x- and y-axis values are the same. On such a plot, a perfect prediction would plot along this line. Checking this box will also cause the module to automatically calculate the RMSE of the data being plotted, assuming that the y-axis represents predicted values, and the x-axis represents known values. This RMSE value will be appended to the text in the legend field, making it simple to compare the results of two different regression models plotted on the same axes. Finally, the user must specify a file name for the plot. Plots are saved as publication quality (1,000 dots per inch) portable network graphic (PNG) files.

Plot Spectra

The Plot Spectra module (fig. 23) allows the user to plot rows of data, typically spectra or similar (for example, peak areas, model coefficients, and so on). As with the Plot module, the user selects a dataset and specifies whether to create a new figure or plot additional data onto an existing figure. If creating a new figure, the axis titles and ranges may be specified. If plotting on an existing figure, these values are already set and are fixed.

The user then selects the x variable for the plot. The options for x variable are drawn from the top-level column names in the specified dataset (meta, comp, or wvl, as shown in fig. 1). After specifying the x variable, the range of x values is automatically set to span the full range of values for that variable. The user may adjust this range to create a plot that is zoomed in on a certain subset of the spectrum. The x-axis title defaults to the name of the x variable selected but may be edited by the user.

Next, the user chooses the column that contains the unique identifier that will be used to choose which row to plot (for example, Target Type from fig. 1). The choices in this drop-down list are drawn from the second-level column names in the dataset selected (Target Name, Target Type, Spectrum ID, and so on in fig. 1). When a column is chosen, the “Choose Rows” list box is populated with all the entries from that column. For example, it is common for a dataset to contain a metadata column containing target names or observation identifications that can be used to identify a spectrum, so the user would choose the “Target Name” column and then select the target desired from the “Choose Rows” list (for example, Olivine or Basalt from fig. 1). As with the x variable, when a row is selected, the y-axis range is automatically set to the minimum and maximum values from the spectrum but may be customized. The y-axis title is specified by the user.

The user may also choose the color of the line to plot and the line style (solid, dashed, dotted, and so on), opacity, and width. The legend for the plot will automatically use the text from the selected row to identify the line being plotted. Finally, the user may specify the file name.

Plot ICA/PCA

The Plot ICA/PCA module is used to visualize both the scores and loadings resulting from PCA or ICA dimensionality reduction (fig. 24). The user selects the dataset and the dimensionality reduction method of interest. The x and y variable lists are then populated with the names of the scores columns resulting from the selected dimensionality reduction. After selecting which scores to plot on the x- and y-axes, the user may choose to color code the points based on the value of a column from the dataset. This can aid in interpreting the scores by visually linking the points to values, such as composition. Finally, the user can specify a file name for the plot.

When run, the module creates a figure that has a scatterplot on the left showing the scores selected and two plots on the right showing the loadings used to calculate the scores as a function of wavelength (fig. 24).

Preprocessing Workflow

The preprocessing workflow begins, as all PyHAT workflows do, by specifying the output folder, followed by loading a dataset from a CSV file. This example uses spectra provided in the repository (libpyhat/data_sets/LIBS/lab_data_example.csv). Once loaded, baseline removal is applied using the airPLS algorithm (fig. 26).

To see how well baseline removal did, we then plot the original spectrum and the baseline identified by the (fig. 27). If desired, additional baseline removal steps that use different parameters or algorithms can be performed on the original dataset and plotted in a similar manner to determine the best choice for the dataset. For the purposes of this example, this tuning process is not included, and the default airPLS settings are used.

Dimensionality Reduction Workflow

In this example, we use the SuperCam spectral library, which is available for download at https://pds-geosciences.wustl.edu/m2020/urn-nasa-pds-mars2020_supercam/calibration_supercam/libs_spectral_library_reference.csv. These data have already had baseline removal applied as part of the SuperCam data processing pipeline. These data have also had wavelength shifts applied, which we do not want to use in this example, so the first step after loading the data is to use the Remove Rows module to remove all rows with shift values other than 0 (fig. 28). Note that because this is a large dataset, the Remove Rows module takes some time to load – be patient! We next normalize the spectra so that each sums to one.

We then use the PCA dimensionality reduction step and visualize the scores and loadings for the first two principal components, which are color coded by the known total calcium oxide (CaO) concentration in each sample (fig. 29).

Cross Validation Workflow

When training a regression model, cross validation should be used to tune the parameters of the algorithms of interest to determine the optimal settings. This workflow begins by loading the preprocessed data from the previous workflow. The data are then stratified on the variable for which the regression model is trained—SiO₂ abundance in this example (fig. 30).

Next, cross validation is performed, which requires choosing a group of algorithms and hyperparameter values for evaluation (fig. 31). For the purposes of this simple example, the OLS and PLS regression algorithms are compared. In practice, testing additional regression algorithms would be more thorough. However, note that running cross validation with all the available regression algorithms and doing so over their many hyperparameters can be time consuming, particularly if local regression is included. Therefore, the user should determine how extensive of a cross-validation routine is required for their purposes.

When cross validation is finished, the results are exported to a CSV file (fig. 32). The lowest RMSECV for OLS in the file is listed as 6.54 weight percent, whereas PLS achieves lower RMSECV in several cases. The lowest RMSECV is not necessarily the optimal choice for PLS, since using too many components can result in an overspecialized model that performs poorly on novel data. To choose the best model, the RMSECV versus number of components for PLS can be plotted (fig. 33). The best model is where the plot forms an elbow, and additional components do not significantly improve the RMSECV. In this case, the optimal number of components is seven.

Blended submodels get better results than a single PLS model can achieve alone, so this cross-validation workflow can be completed two more times, changing the range of compositions used for cross validation to 0–60 and 50–100 weight percent. These models are optimized with seven and eight components, respectively.

Regression Training and Prediction Workflow

In this example workflow, the results from cross validation are used to train the regression models, optimize the blending process, and then predict the test set. First, preprocessed data need to be loaded and the stratified folds method selected to separate the training and test data.

Regression models are then trained (fig. 34), using the parameters determined by cross validation, and those models are used to predict the train and test set compositions (fig. 35).

Once predictions are complete, the Blended Submodel Predictions module can be used with the training set predictions first (fig. 36). The 0–100 weight percent PLS model in this case is used as the reference model, the 0–60 weight percent model as the low model, and the 50–100 weight percent model as the high model. When optimizing the blending by comparing the results with the known SiO₂ compositions, the lowest RMSE is achieved with a blending range from 27.71 to 64.16 weight percent. This means, for reference model predictions below 27.71 weight percent, the predictions from the low model will be used, above 64.16 weight percent the high model will be used, and in between the results will be blended using a linear weighted sum. Next, the same blending range is applied to the test set predictions, and the test set results are plotted against the known compositions to evaluate the accuracy of the final blended predictions (fig. 37).

Endmember Identification and Spectral Unmixing Workflow

PyHAT is not only for analysis of planetary datasets, it is useful for terrestrial remote sensing as well. In this example, endmembers are identified from a subset of the Salinas Scene dataset (see the Endmember Identification section above for more information about the dataset). Note that because there can be considerable spectral variability for a single ground truth (for example, lettuce), the endmember identification algorithm may not always identify unique ground truths. Instead, it may include spectral extremes of a single ground truth (for example, old and young lettuce crop).

First, the Salinas Scene dataset multicolumn test file “labeled_Salinas_testfile.csv,” provided as part of the PyHAT repository, is loaded and named. For this example, the spectra are also normalized using the normalize functionality in the preprocessing drop-down menu. The variable to normalize is set to “wvl.” The spectra are normalized from a minimum wavelength of 0 to a maximum of 2,500 (the actual wavelengths of the data are 360 to 2,500 nm). Using the data drop-down menu, “Identify endmembers” is selected. For this example, the PPI algorithm is used, and the number of endmembers is set to three. When endmember identification is run, a column is added to the dataset.

To unmix the composite Salinas spectra, the “Unmixing” option in the preprocessing drop-down menu is used. In this case, the same dataset for both the data to unmix and the endmembers is used because the endmembers were added to the dataset. The unmixing method used in this example is FCLS. The “Normalize” box is left unchecked because the spectra were normalized in the previous steps. Figure 38 shows a screenshot of the unmixing workflow, and figure 39 shows a plot of the three identified endmember spectra.