data_explore_clean.Rmd

---
title: "Code to clean and explore data"
author: "Caitlin Mothes"
---

# Workflow to clean and prep all data for the shiny app

```{r}
library(tidyverse)
library(readxl)
library(sf)
library(mapview)
```

## Read in datasets

```{r}
# Data from Kaela
sample <- read_xlsx("data/SAMPLE_geospatial_metadata_3.28.24.xlsx") %>% 
  janitor::clean_names()

genome <- read_xlsx("data/GENOME_geospatial_metadata_02.17.2026.xlsx")
  

mag_abundance <- read_xlsx("data/MAGs_average_relativeabund_by_basin.xlsx")

mag_cores <- readxl::read_xlsx("data/core_ALL.xlsx")
# Basins and Plays

# Downloaded from here: https://atlas.eia.gov/search?q=plays&source=u.s.%2520energy%2520information%2520administration
play_basin <- read_sf("data/TightOil_ShaleGas_Plays_Lower48_EIA.shp") %>% 
  st_transform(crs = 4326) %>% 
  # remove spaces around hyphens for joining and change Denver basin name
  mutate(Shale_play = str_replace(Shale_play, " - ", "-"),
         Basin = str_replace(Basin, " Basin", ""))

sediment_basin <- read_sf("data/SedimentaryBasins_US_EIA.shp") %>% 
  mutate(NAME = str_to_title(NAME)) %>% 
  # need to transform all of these to WGS84
  st_transform(crs = 4326)

```

## Explore International Basins

```{r}
bowland <- read_sf("data/BGS_DECC_Bowland_study/topBHU_limit.shp") %>% 
  mutate(shale_play = "uk", shale_basin = "Bowland Shale") %>% 
  st_transform(4326)
mapview(bowland)


west_canada <- read_sf("data/wcsb_shapefiles/WCSBAtlasShp/fg0301_py_ll.shp")
mapview(west_canada)
```

Quick map of basins

```{r}
mapview(play_basin) +
  mapview(sediment_basin, color = "red")
```

## Sample Explorer Data

Filter plays to join to sample df to get centroids:

```{r}
play_basin_filtered <- play_basin %>%
  st_make_valid() %>% 
  group_by(Shale_play, Basin) %>% 
  summarise(st_union(geometry)) %>% 
  ungroup() %>% 
  filter(Shale_play %in% sample$shale_play) %>% 
  rename(shale_play = Shale_play, shale_basin = Basin)
```

Create sample layers for app

```{r}
sample_sf <- sample %>%
  left_join(play_basin_filtered, by = c("shale_play", "shale_basin")) %>%
  st_as_sf() %>%
  # convert to points
  st_centroid() %>%
  # edit lithology typo
  mutate(
    # make days_since_frack numeric
    days_since_frack = na_if(days_since_frack, "N/A|unk") %>%
      as.numeric(.),
    # edit time series stage for filtering
    timeseries_stage = if_else(timeseries_stage == "N/A" | timeseries == "no", "none", timeseries_stage),
    #edit salinity values for filtering
    #salinity_classification = if_else(salinity_classification == "N/A", NA_character_, salinity_classification),
    salinity_conductivity_m_s_cm = if_else(salinity_conductivity_m_s_cm == "not measured", NA_character_, salinity_conductivity_m_s_cm) %>% 
      as.numeric(.),
    # edit 'N/A' values and make numeric
    across(contains("perc"), ~if_else(. == "N/A", NA_character_, .) %>% as.numeric(.)))

# create variables for map symbology
sample_app <- sample_sf %>%
  st_drop_geometry() %>%
  group_by(well_id) %>%
  mutate(
   # n_samples = n(),
    max_days_since_frack = max(days_since_frack, na.rm = TRUE),
    min_days_since_frack = min(days_since_frack, na.rm = TRUE),
    # change 'Inf' values to 0
    max_days_since_frack = if_else(max_days_since_frack == -Inf, 0, max_days_since_frack),
    min_days_since_frack = if_else(min_days_since_frack == Inf, 0, min_days_since_frack)
  ) %>%
  distinct(well_id, .keep_all = TRUE) %>%
  select(shale_play,
         well_id, 
         #n_samples,
         min_days_since_frack, 
         max_days_since_frack) %>%
  right_join(sample_sf, by = c("shale_play", "well_id")) %>%
  ungroup() %>%
  st_as_sf()

# separate samples from inputs
sample_inputs <- sample_app %>% 
  filter(sample_type == "input")

sample_app <- sample_app %>% 
  filter(sample_type != "input")
```

## Genome Explorer Data

### Relative abudance data

Tie abundance data to genome data

```{r}
genome_joined <- mag_abundance %>% 
  dplyr::select(MAG_ID = MAG, rel_abundance = aveRelAbund, Play, Basin) %>% 
  left_join(genome, by = c("MAG_ID")) %>% 
  select(-basin)
```

Clean genome data and join to basin shapefile

```{r}
genome_app <- genome_joined %>% 
  mutate(across(domain:species, ~str_remove(.x, ".*__")))
```

Visualize for app:

```{r}
# **will need to do this within the app using filtered data for 'genome_clean'

# basins 
genome_app %>% 
  filter(genus == "Thermotoga") %>% 
  group_by(Basin, species) %>% 
  summarise(max_rel_abundance = max(rel_abundance)) %>% 
  group_by(Basin) %>% 
  summarise(avg_rel_abundance = mean(max_rel_abundance)) %>% 
  right_join(sediment_basin, by = c("Basin" = "NAME")) %>%
  View()


# plays
genome_app %>% 
 group_by(Basin, Play, domain, phylum, class, order, family, genus, species) %>% 
  summarise(max_rel_abundance = max(rel_abundance)) %>% 
  group_by(Basin, Play) %>% 
  summarise(avg_rel_abundance = mean(max_rel_abundance)) %>% 
  inner_join(play_basin, by = c("Play" = "Shale_play")) %>%
  st_as_sf() %>% 
  View()


# test creating column for popup

# Old CODE
# basin_genome <- play_basin %>%
#   # trailing spaces causing errors in join
#   mutate(
#     Basin = str_remove_all(Basin, " ") %>%
#       str_remove_all(., "Basin"),
#     # fix typo
#     Lithology = if_else(Lithology == "MIxed Shale & Chalk", "Mixed Shale & Chalk", Lithology)
#   ) %>%
#    st_make_valid() %>% 
#   group_by(Basin) %>% 
#   summarise(geometry = st_union(geometry)) %>% 
#   ungroup() %>% 
#   st_make_valid() %>% 
#   # join to sample
#   inner_join(genome_clean, by = "Basin") %>%
#   group_by(Basin) %>% 
#   summarise(n_MAG_samples = n())
  
```

### Core Data

```{r}

# prep for app

## basins
core_app <- mag_cores %>% 
  filter(perc_samples_present_per_basin >= 0.50) %>% 
  group_by(basin) %>% 
  count() %>% 
  # join to shapefile
  inner_join(sediment_basin, by = c("basin" = "NAME")) %>% 
  st_as_sf()


mapview(core_app)
```

## Save all data for app

```{r}
#save updated data
save(sample_app, sample_inputs, play_basin, sediment_basin, genome_app, mag_cores, file = "data/app_data_Update.RData")

```

## Timeseries Plot

```{r}

fig_1 <- sample_app %>% 
plotly::plot_ly(colors = "Dark2") %>% 
  add_trace(x = sample_app$days_since_frack,
            y = sample_app$well_id,
            split = ~sample_app$well_id,
            color = ~sample_app$shale_basin,
            name = ~sample_app$shale_basin,
            legendgroup = ~sample_app$shale_basin,
            type = 'scatter',
            mode = 'lines+markers',
            connectgaps = TRUE) %>% 
  layout(showlegend = TRUE, yaxis = list(type = "category",
  categoryorder = "array", categoryarray = sort(unique(sample_app$shale_basin)))) %>% 

  # hacky way to get one trace per group in legend
  style(showlegend = FALSE, traces = c(1:4, 7:11, 13:16, 17, 21, 22:24, 26, 28))

fig_2 <-  sample_app %>% 
plotly::plot_ly(colors = "Dark2") %>% 
  add_trace(x = sample_app$days_since_frack,
            y = sample_app$well_id,
            split = ~sample_app$well_id,
            color = ~sample_app$shale_basin,
            name = ~sample_app$shale_basin,
            legendgroup = ~sample_app$shale_basin,
            type = 'scatter',
            mode = 'lines+markers',
            connectgaps = TRUE) %>% 
  layout(showlegend = TRUE, yaxis = list(type = "category",
  categoryorder = "array", categoryarray = sort(unique(sample_app$shale_basin)))) %>% 
  # hacky way to get one trace per group in legend
  style(showlegend = FALSE, traces = c(1:4, 7:11, 13:16, 17, 21, 22:24, 26, 28))

subplot(fig_1, fig_2, nrows = 1, shareY = TRUE, margin = 0.02) %>% 
  layout(xaxis = list(range = c(0, 2500)),
         xaxis2 = list(range = c(4500, 5000)))
  
```

## Well Area Plots

```{r}
# read in data
load("data/app_data_Update.RData")

# function to create and save all plots
area_plot <- function(x, data) {
  data %>%
    filter(well_id == x) %>%
    ggplot(aes(x = days_since_frack, y = value, fill = group)) +
    geom_area() +
    scale_fill_manual(
      values = c("#de3064", "#077969", "#1052ba"),
      name = "Type",
      labels = c("Acetate Producers", "Methanogens", "Sulfide Producers")
    ) +
    theme_minimal() +
    labs(y = "% Relative Abundance", x = "Day")+
    theme(text=element_text(size = 16))
  
  ggsave(
    paste0("www/", x, "_area_plot.png"),
    width = 6,
    height = 4,
    dpi = 300,
    units = "in"
  )
  
}


# arrange/clean well data
plot_data <- sample_app %>%
  st_drop_geometry() %>%
  select(well_id, days_since_frack, contains("perc")) %>%
  pivot_longer(cols = contains("perc"), names_to = "group")


# create and save plot for each well
walk(unique(plot_data$well_id), area_plot, plot_data)


area_plot("AND-1", plot_data)

ggsave("test_area_plot.png", width = 6, height = 4, dpi = 300, units = "in")
```