LIBS_starter_Jul2024.Rmd

---
title: 'DAR Starter Notebook: M20 LIBS Analysis'
author: "Kristin Bennett and John Erickson"
date: "`r format(Sys.time(), '%d %B, %Y')`"
output:
  html_document:
    df_print: paged
---

# Introduction

This notebook demonstrates some basic analysis of data from the M20 Laser-Induced Breakdown Spectroscopy (LIBS) experiment, part of the data collected by the Perseverance Rover's SuperCam instrument. 

* More about SuperCam: https://an.rsl.wustl.edu/help/Content/About%20the%20mission/M20/Instruments/M20%20SuperCam.htm
* Raw data source: https://pds-geosciences.wustl.edu/m2020/urn-nasa-pds-mars2020_supercam/data_derived_spectra/supercam_libs_moc.csv

The data we'll be working with has had Mars latitude and longitude added to it. We load it as a binary Rds file to save


```{r include = FALSE}

# Package loads
library(tidyverse)
library(stringr)
library(BBmisc)
library(pheatmap)
library(ggplot2)
library(ggtern) # For ternary plots

```

# NOTES: 

* NASA Analysis Notebooks at: https://an.rsl.wustl.edu/m20/AN/an3.aspx

* About SuperCam: 

* LIBS source: https://pds-geosciences.wustl.edu/m2020/urn-nasa-pds-mars2020_supercam/data_derived_spectra/supercam_libs_moc.csv

# Load LIBS Data

```{r}
# Saved LIBS data with locations added
supercam_libs_moc_loc <- readRDS("/academics/MATP-4910-F24/Mars-F24/Data/supercam_libs_moc_loc.Rds")

# Add "Unknown" column
supercam_libs_moc_loc <- supercam_libs_moc_loc %>% 
  mutate(Unknown=100-Total)
```

# Create Ternary Plot

```{r}

# See e.g. https://serc.carleton.edu/mathyouneed/geomajors/ternary/index.html 

# Select just the essentials
# "Unknown" not used for ternary, but is for PCA
supercam_trim <- supercam_libs_moc_loc %>% select(c(target, SiO2, Al2O3, FeOT, MgO, CaO, Na2O, K2O, Unknown, sol))

libs_ternary <- supercam_trim %>%
  mutate(x=(SiO2+Al2O3)/100,y=(FeOT+MgO)/100,z=(CaO+Na2O+K2O)/100, value=sol) %>%
  select(-c(target,SiO2,Al2O3,FeOT,MgO,CaO,Na2O,K2O,Unknown,sol)) %>%
  drop_na() # Just in case...

libs_ternary.mat <- as.matrix(libs_ternary[,1:3])

tern.km <- kmeans(libs_ternary.mat, centers=4)

libs_ternary <- cbind(libs_ternary, cluster=as.factor(tern.km$cluster))

ggtern(libs_ternary, ggtern::aes(x=x,y=y,z=z,cluster=cluster)) +
  geom_point(aes(color=cluster)) + 
  theme_rgbw() +
  labs(title="Example of Mars 2020 LIBS ternary Plot",
       x="Si+Al",
       y="Fe+Mg",
       z="Ca+Na+K")
```

# PCA on LIBS Data

```{r}

# Include "Unknown" in matrix
libs_full.mat <- as.matrix(supercam_trim[2:9])

libs_full.mat.pca <- prcomp(libs_full.mat, scale=TRUE)

# Scree plot
eigenvalues <- libs_full.mat.pca$sdev^2
plot(eigenvalues/sum(eigenvalues), type = "b",
     xlab = "Principal Component",
     ylab = "Percentage of Variance Explained")

k <- 4
km <- kmeans(libs_full.mat,k)

# Prepare df for plotting
plot.df <- cbind.data.frame(PC1=libs_full.mat.pca$x[,1], PC2=libs_full.mat.pca$x[,2],
                            Cluster=as.factor(km$cluster))

expvar.pc1 <- summary(libs_full.mat.pca)$importance[2,1]
expvar.pc2 <- summary(libs_full.mat.pca)$importance[2,2]


libs_pca_plot <- ggplot() +
  ## In case we want ellipses around clusters
  # ggforce::geom_mark_ellipse(data=plot.df,
  #                            aes(x=PC1, y=PC2,
  #                                group=Cluster,
  #                                label=Cluster
  #                            ),
  #                            show.legend=FALSE) + 
  geom_point(data=plot.df,
             aes(x=PC1, y=PC2,
                 group=Cluster,
                 color=Cluster
             ),
             size=1) +
  # geom_text_repel(data=plot.df, 
  #                 aes(x=PC1, y=PC2)) +
  coord_fixed(ratio=1) +
  geom_hline(yintercept = 0, color = "gray70") +
  geom_vline(xintercept = 0, color = "gray70") +
  xlim(-5,5) + ylim(-6,6) +
  ggtitle('Mars LIBS PCA Biplot') +
  geom_segment(aes(x=0, y=0,
                   xend = 7*libs_full.mat.pca$rotation[,1],
                   yend = 7*libs_full.mat.pca$rotation[,2]),
               arrow = arrow(length = unit(1/2, 'picas'))) +
  geom_text(aes(x = 7.5*libs_full.mat.pca$rotation[,1],
                y = 7.5*libs_full.mat.pca$rotation[,2],
                label=rownames(libs_full.mat.pca$rotation)),
            size=4) +
  xlab(paste0("PC1: ",expvar.pc1*100,"% explained variance")) +
  ylab(paste0("PC2: ",expvar.pc2*100,"% explained variance")) + 
  ggpubr::theme_pubclean(base_size = 15)

libs_pca_plot
```