Add files via upload

Author: maggiemmcadam-arc

BD_asteroid_characteristics.xlsx

@@ -0,0 +1,173 @@
+"""
+Working Python file to import Bus-DeMeo spectral data, Bus-DeMeo asteroid 
+sekected orbital and physical properties, import balanced training dataset, 
+and classify the Bus-DeMeo asteroids using the methodology described in the paper.
+ 
+
+"""
+import os 
+import pathlib
+import pandas as pd
+import numpy as np
+from scipy import interpolate
+from sklearn import preprocessing
+from matplotlib import pyplot as plt
+from pySankey.sankey import sankey
+from specmap import *
+
+
+#%% load balanced training dataset
+training = pd.read_pickle('balanced_training_dataset2026_01_20_1620')
+C = 64
+
+#%% Load asteroids data
+ast_loc = 'DeMeo_csv_files' #provide the location of local folder where the asteroid data is saved.
+#%% Load Bus-DeMeo asteroid spectra into a dataframe, remove NANs
+infolder = list(pathlib.Path(ast_loc).glob('*.csv'))
+count = len(infolder)
+wave_ = training.columns.tolist()#bal.columns.tolist()#
+nir_start = np.where(np.array(wave_[2:]) >= 0.44)
+wave =wave_[nir_start[0][0]+2:]
+BDdata = np.empty([len(wave),count])
+BDdata[:] = np.nan
+BDdata_sm = pd.DataFrame()
+BDdata_err = np.empty([len(wave),count])
+ast_desig = []
+count = 0
+wave_col = 1
+spec_col = 2
+err_col = 3
+
+for i in range(0,len(infolder)):
+    path,file_name = os.path.split(infolder[i])
+    temp = pd.read_csv(infolder[i],skiprows = 60,header=None)
+    ast_desig.append(file_name[0:-9])
+    if temp.iloc[0,wave_col] > wave[0] and  temp.iloc[-1,wave_col] <= wave[-1]:
+        start_wave = np.where(wave >= temp.iloc[0,wave_col])
+        end_wave = np.where(wave >= temp.iloc[-1,wave_col])
+        
+        sp_interp = interpolate.interp1d(temp.iloc[:,wave_col], temp.iloc[:,spec_col])
+        err_interp = interpolate.interp1d(temp.iloc[:,0], temp.iloc[:,err_col])
+        wave_range = wave[start_wave[0][0]:end_wave[0][0]]
+        b_ = sp_interp(wave_range)
+        c_ = err_interp(wave_range)
+        BDdata[start_wave[0][0]:end_wave[0][0],count] = b_
+        BDdata_err[start_wave[0][0]:end_wave[0][0],count] = c_
+        
+    if temp.iloc[0,wave_col] <= wave[0] and temp.iloc[-1,wave_col] <= wave[-1]:
+        start_idx = np.where(temp.iloc[:,wave_col].round(2) >= wave[0])
+        end_wave = np.where(wave >= temp.iloc[-1,wave_col])
+        if start_idx[0][0] == 0:
+            sp_interp = interpolate.interp1d(temp.iloc[start_idx[0][0]:,wave_col], temp.iloc[start_idx[0][0]:,spec_col])
+            err_interp = interpolate.interp1d(temp.iloc[start_idx[0][0]:,wave_col], temp.iloc[start_idx[0][0]:,err_col])
+        else:
+            sp_interp = interpolate.interp1d(temp.iloc[start_idx[0][0]-1:,wave_col], temp.iloc[start_idx[0][0]-1:,spec_col])
+            err_interp = interpolate.interp1d(temp.iloc[start_idx[0][0]-1:,wave_col], temp.iloc[start_idx[0][0]-1:,err_col])
+                
+        wave_range = wave[0:end_wave[0][0]-1]
+        b_ = sp_interp(wave_range)
+        c_ = err_interp(wave_range)
+        BDdata[0:end_wave[0][0]-1,count] = b_
+        BDdata_err[0:end_wave[0][0]-1,count] = c_
+    
+    if temp.iloc[0,wave_col] <= wave[0] and temp.iloc[-1,wave_col] >= wave[-1]: 
+        start_idx = np.where(temp.iloc[:,wave_col].round(2) == wave[0])
+        end_idx = np.where(temp.iloc[:,wave_col] >= wave[-1])
+        if temp.iloc[start_idx[0][0],0] > wave[0]:
+            sp_interp = interpolate.interp1d(temp.iloc[start_idx[0][0]-1:end_idx[0][0]+1,wave_col], temp.iloc[start_idx[0][0]-1:end_idx[0][0]+1,spec_col])
+            err_interp = interpolate.interp1d(temp.iloc[start_idx[0][0]-1:end_idx[0][0]+1,wave_col], temp.iloc[start_idx[0][0]-1:end_idx[0][0]+1,err_col])
+        else: 
+            sp_interp = interpolate.interp1d(temp.iloc[start_idx[0][0]:end_idx[0][0]+1,wave_col], temp.iloc[start_idx[0][0]:end_idx[0][0]+1,spec_col])
+            err_interp = interpolate.interp1d(temp.iloc[start_idx[0][0]:end_idx[0][0]+1,wave_col], temp.iloc[start_idx[0][0]:end_idx[0][0]+1,err_col])
+        
+        b_ = sp_interp(wave)
+        c_ = err_interp(wave)
+        BDdata[:,count] = b_
+        BDdata_err[:,count] = c_
+        print(f'for {count}, condition 3')
+    count += 1 
+    print(f'count = {count}')
+
+BDdata = pd.DataFrame(BDdata)
+BDdata.columns = ast_desig
+
+BDdata.bfill(inplace = True)
+BDdata.ffill(inplace = True)
+#%% Classify Bus-DeMeo asteroids using the balanced training dataset and a logistic regression classifier.
+
+wave_ = training.columns.tolist()
+nir_start = np.where(np.array(wave_[2:]) >= 0.44)
+wave =wave_[nir_start[0][0]+2:]
+
+train_spec_start_idx = 2
+desig = ast_desig
+test_wave = wave
+test_data = BDdata
+norm_wave = 1
+
+BDpredictions = classify(training, train_spec_start_idx, desig, test_wave, test_data, norm_wave, C)
+#%% Import selected characteristics of asteroids and save the classification results
+bd_char_file = "C:\\Users\\maggi\\Dropbox\\Work\\Research\\specmap\\BD_asteroid_characteristics.xlsx"  
+bd_char =       pd.read_excel(bd_char_file)
+bd__ = pd.concat([bd_char,BDpredictions],axis = 1)
+
+#%% Defining confidence of predictions
+# this piece of code takes the pred_proba values which are generated in the classifier above
+# and puts them into a new list; these values are the logistic regression's predicted probability
+# that the spectrum being classified is in each group. To define a confidence of the prediction
+# we can assess how much each prediction differes from the others. If it's a strong
+# prediction, the standard deviation will be high: one very large value and 9 values close to zero
+# If the prediction is 'weak' or less confidence, the probabilities of the groups will be similar
+# and the standard deviation will be small.
+bd_std = []
+bd_pred_val = []
+for i in range(0,bd__.shape[0]):
+    bd_pred_val.append(np.max(bd__.iloc[i,19:]))
+    bd_std.append(np.std(bd__.iloc[i,19:]))
+    
+bd_std = pd.DataFrame(bd_std)
+bd_pred_val = pd.DataFrame(bd_pred_val)
+bd_ = pd.concat([bd__,bd_pred_val],axis=1)
+bd_ = pd.concat([bd_,bd_std],axis = 1)
+
+label_encoder = preprocessing.LabelEncoder()
+
+bd = bd_.sort_values(by = 'BD-complex',ascending = False)
+
+#%% Rename groups from numerical identifiers into words and sort dataframe by spectral complex
+
+
+
+label_encoder = preprocessing.LabelEncoder()
+
+bd = bd_.sort_values(by = 'BD-complex',ascending = False)
+
+bd.loc[bd['pred'] == 1,'pred'] = 'Hydrated CCs'
+bd.loc[bd['pred'] == 2,'pred'] = 'CO/CV'
+bd.loc[bd['pred'] == 3,'pred'] = 'CK/R/Brach'
+bd.loc[bd['pred'] == 4,'pred'] = 'H/L/LL/Ure'
+bd.loc[bd['pred'] == 5,'pred'] = 'EH/EL/Aub'
+bd.loc[bd['pred'] == 6,'pred'] = 'ACA/LOD'
+bd.loc[bd['pred'] == 7,'pred'] = 'Irons'
+bd.loc[bd['pred'] == 8,'pred'] = 'HEDs'
+bd.loc[bd['pred'] == 9,'pred'] ='CY'
+bd.loc[bd['pred'] == 10,'pred'] = 'Primitive CCs'
+
+#%%Create a Sankey diagram for all the endmember spectral groups: 
+others = ['T', 'R', 'Q', 'L', 'K', 'D', 'A']
+bd_t = bd[~bd['BD-complex'].isin(['C'])]
+bd_te = bd_t[~bd_t['BD-complex'].isin(['S'])]
+bd_tes = bd_te[~bd_te['BD-complex'].isin(['X'])]
+bd_test = bd_tes[~bd_tes['BD-complex'].isin(['V'])]
+
+bd_2 = pd.concat([bd_test['spec_B'],bd_test['pred']],axis = 1)
+
+sankey(bd_2['spec_B'],bd_2['pred'],fontsize = 7)
+plt.show()
+
+#%% Create a Sankey diagram for the C-complex asteroids
+
+bd_test = bd[bd['BD-complex'].isin(['C'])]
+
+sankey(bd_2['spec_B'],bd_2['pred'],fontsize = 7)
+plt.show()

Bus-DeMeo-results

@@ -1,2 +1,56 @@
-# SpectralMeteoriteMappingTool
-The Spectral Meteorite Mapping Tool is a machine learning algorithm designed to analyze asteroid mineralogy using archived near-infrared meteorite spectroscopy.
+# Overview: 
+
+This repository contains the functions, data and scripts that produced the results in two publications: 
+
+McAdam, M. M., Dotson, J., Dyar, M. D., Hidden in plain sight: a new Main Belt population of aqueously altered and thermally metamorphosed asteroids. _Submitted to Earth and Space Sciences Journal May 2026_
+
+McAdam, M. M., Dotson, J., Dyar, M. D., Spectral Meteorite Mapping Tool for Asteroid Composition Identification: uses in Planetary Defense and General Asteroid Science. (in preparation)
+
+Please contact Dr. McAdam with any questions regarding the code and data: margaret.m.mcadam at nasa dot gov
+
+### Dependencies: 
+The code and scripts were created using Python version 3.11.11. 
+I used conda: spyder-env integrated development environment to write, test and run the code and scripts presented here. 
+
+To run the functions and scripts, you'll need the following packages: 
+- Numpy
+- Pandas
+- SciPy 
+- SciKitLearn
+- MatplotLib 
+- PySankey
+
+  
+### Functions and scripts
+
+**specmap.py** - The Spectral Meteorite Mapping Tool are several functions that create balanced training datasets and use a supervised machine learning approach (logistic regression) to classify data. The logistic regression is trained using the provided balanced dataset and indicated hyperparameter. These functions are used to classify asteroid or meteorite spectra into one of 10 groups based on their spectral similarity.
+
+ _classify_ function is the function that classifies near-infrared spectroscopy into one of 10 groups of meteorites using a logistic regression machine learning algorithm trained on near-infrared archival meteorite spectra. This function can be used to classify any near-infrared dataset that has the same sampling and wavelength range of the Bus-DeMeo asteroid spectra saved in the _DeMeo_csv_files_ folder.
+
+_balanced_generator_ function. The training data has uneven classes and so we elected to balance our training dataset, over- and under-sampling the meteorite classification groups to create a balanced training dataset. We generated multiple balanced training datasets and selected the one that performed best for the iron and hydrated carbonaceous chondrite classification groups. The function that generates the balanced training dataset is called 'balanced_generator' in the specmap.py code. The balanced_generator function uses the data saved in trainTestSplitTraining.xlsx to created balanced training datasets.
+
+There are two other functions in this file, _normalize_ and _renormalize_, which are called in the function 'classify'. Data normalization is a rescaling of data. In the spectroscopy context, we divide the whole spectrum by one particular datapoint in the spectrum; that point becomes unity and the rest of the data have the same attributes (absorption bands, slope etc.) as the raw data but are now scaled. Spectral normalization allows us to directly compare spectral features of different spectra. The _renormalize_ function takes any data, including previously normalized data, and renormalizes it to the user-defined normalization point. The _normalize_ function takes unnormalized data and normalizes it at the user-defined wavelength.  
+
+**BusDeMeoReanalysis.py** - This file is the script used to ingest, classify and generate the results presented in the publications listed above. We have not included all the plots for the figures in the papers but two of the Sankey diagrams are demonstrated. 
+
+### Pickle files:
+
+**balanced_training_dataset2026_01_20_1620** - This file is the balanced training dataset that was selected based on its performance classifying the test dataset. We generated 30 balanced datasets then assessed their performance on the test data. This training dataset was selected based on how well it classified hydrated carbonaceous chondrites and iron meteorites. The C parameter for this training dataset is 64. This was obtained by a cross-validation analysis for this balanced dataset where we determined the best overall accuracy for this C parameter.
+
+**Bus-DeMeo-results** - This pickle includes the assigned classification and group probability for each asteroid in the dataset. Can be loaded independently of the BusDeMeoReanalysis.py and used as needed.
+
+### Data files
+Folder **DeMeo_csv_files** 
+The folder containing the Bus-DeMeo asteroid spectra. These data are loaded into a dataframe in the working.py file and classified using the balanced training dataset included in this repo.
+
+**trainTestSplitTesting.xlsx** - This file contains 20% of the spectral dataset (including the original 1422 meteorite spectra from Dyar et al., 2023 as well as the new data). We stratified the dataset (ordering by unnormalized reflectance value at 1.25-µm) and selected every 5th spectrum. These spectra we reserved as a test cohort to validate the accuracy of our classifier.
+
+**trainTestSplitTraining.xlsx** - This file contains 80% of the dataset. These data are used to generate a balanced training dataset then to train the logistic regression classifier.
+
+**Dyar_etal(2023)_dataset.xlsx** - This file is the original 1422 meteorite spectra. It is a duplicate of the supplemental data file from Dyar et al., 2023. Citation: Dyar, M.D., Wallace, S.M., Burbine, T.H. and Sheldon, D.R., 2023. A machine learning classification of meteorite spectra applied to understanding asteroids. Icarus, 406, p.115718.
+
+**new_RELAB_spectra_resampled.xlsx** - This file includes the ~600 new spectra that were put on the RELAB archive after Dyar et al., (2023) created their dataset. These spectra augment the original dataset.
+
+**BD_asteroid_characteristics.xlsx** - This file contains selected orbital and physical parameters of the 371 asteroids that we reanalyze in the paper. These are all publicly available data that were gathered from the JPL Solar System Dynamics Small Bodies search query (ssd.jpl.nasa.gov). 
+
+**Bus-DeMeo-results_vfinal.xlsx** - This file contains the physical characteristics as well as the predicted class and probability of each class for all 371 asteroids in the study.

Bus-DeMeo-results_vfinal.xlsx

@@ -0,0 +1,25 @@
+Overview: 
+
+# Folder 'DeMeo_csv_files' - The folder containing the Bus-DeMeo asteroid spectra. These data are loaded into a dataframe in the working.py file and classified using the balanced training dataset included in this repo.
+
+specmap.py - The Spectral Meteorite Mapping Tool are several functions that create balanced training datasets and use a supervised machine learning approach (logistic regression) to classify data. The logistic regression is trained using the provided balanced dataset and indicated hyperparameter. These functions are used to classify asteroid or meteorite spectra into one of 10 groups based on their spectral similarity.
+
+Working.py - This file is an example of how the software is called and used in python. This file can be localized by other scientists to use and understand our software.
+
+# Pickle files:
+
+balanced_training_dataset2026_01_20_1620 - This file is the balanced training dataset that was selected based on its performance classifying the test dataset. We generated 30 balanced datasets then assessed their performance on the test data. This training dataset was selected based on how well it classified hydrated carbonaceous chondrites and iron meteorites. The C parameter for this training dataset is 64. This was obtained by a cross-validation analysis for this balanced dataset where we determined the best overall accuracy for this C parameter.
+
+Bus-DeMeo-results - This pickle includes the assigned classification and group probability for each asteroid in the dataset. Can be loaded independently of the working.py and used as needed.
+
+BD_asteroid_characteristics.xlsx - This file contains selected orbital and physical parameters of the 371 asteroids that we classify in the paper. These are all publicly available data that were gathered from the JPL Solar System Dynamics Small Bodies search query (ssd.jpl.nasa.gov). 
+
+Dyar_etal(2023)_dataset.xlsx - This file is the original 1422 meteorite spectra. It is a duplicate of the supplemental data file from Dyar et al., 2023. Citation: Dyar, M.D., Wallace, S.M., Burbine, T.H. and Sheldon, D.R., 2023. A machine learning classification of meteorite spectra applied to understanding asteroids. Icarus, 406, p.115718.
+
+new_RELAB_spectra_resampled.xlsx - This file includes the ~600 new spectra that were put on the RELAB archive after Dyar et al., (2023) created their dataset. These spectra augment the original dataset.
+
+trainTestSplitTesting.xlsx - This file contains 20% of the spectral dataset (including the original 1422 meteorite spectra from Dyar et al., 2023 as well as the new data). We stratified the dataset (ordering by unnormalized reflectance value at 1.25-µm) and selected every 5th spectrum. These spectra we reserved as a test cohort to validate the accuracy of our classifier.
+
+trainTestSplitTraining.xlsx - This file contains 80% of the dataset. These data are used to generate a balanced training dataset then to train the logistic regression classifier.
+
+Bus-DeMeo-results_vfinal.xlsx - This file contains the physical characteristics as well as the predicted class and probability of each class for all 371 asteroids in the study.