In raphidoc/lighthouse: Database management and processing workflow for ocean optics science

library(dplyr)
library(tidyr)
library(DT)
library(stringr)
library(lubridate)
library(leaflet)
library(crosstalk)
library(plotly)
library(purrr)
library(leafem)
#library(mapview)
library(leaflet.extras)
library(htmltools)
library(leafpop)


knitr::opts_knit$set(root.dir = "~/R/lighthouse/R/")
knitr::opts_chunk$set(echo = F, cache = F, message = F, warning = F)

load("/mnt/D/Documents/Maitrise/Paper/Models/Data/models_db.RData")

# models_db <- models_db %>% mutate(ROI= str_extract(PID,"[:alpha:]+[:digit:]{0,1}"))

# Create a DataLog to keep track of comments on data and outlier removal
if (!file.exists("/mnt/D/Documents/Maitrise/Paper/Models/DataLog.csv")) {
    DataLog <- models_db %>% select(ID,PID) %>%
        mutate(QC_SPM=1,
              QC_Ag=1,
              QC_Bbp=1,
              Comment="")

    write.csv(DataLog, "/mnt/D/Documents/Maitrise/Paper/Models/DataLog.csv", row.names=F)

} else {
    DataLog <- readr::read_csv("/mnt/D/Documents/Maitrise/Paper/Models/DataLog.csv")
}

# Check for bottom effect, linearly based on secchi depth
br_check <- function(Zsecchi,Zstation){

    if (is.na(Zsecchi) | is.na(Zstation)) {
        NA
    } else if (Zsecchi > Zstation){
        # Case when error of entry
        NA
    } else if (Zsecchi == 999) {
        # Case when Secchi touch the bottom
        T
    } else if (Zsecchi*2 > Zstation) {
        # assume a linear travel of light in water...
        T
    } else {
        F
    }
    #str_c(Zstation, Zsecchi)
}

models_db <- models_db %>% left_join(DataLog %>% select(ID,QC_SPM,QC_Ag,QC_Bbp,Comment), by="ID") %>%
    mutate(Region = ifelse(str_detect(Project,"CHONe|WISE|PMZA"),"EGSL","JB"),
          OptShallow = br_check(Zsecchi,Zstation),
          PIM_frac = ifelse((PIM/SPM)*100 <= 100,(PIM/SPM)*100,NA)) %>%
    relocate(matches("QC|Comment|ROI|Region"), .after = "DateTime") %>%
    relocate("OptShallow", .after = "Zsecchi")


# Simulated data for matchup and validation
wide_sensor <- models_db %>% 
    select(!c("data")) %>%
    unnest(cols = c(OLI,S2A), names_sep="_")

wide_sensor <- wide_sensor %>% rename_with(~str_replace(.,"S2A","MSI"), starts_with("S2A")) %>%
    mutate(across(matches("OLI|MSI"), ~na_if(., .<0)))

long_sensor <- wide_sensor %>%
    pivot_longer(cols = all_of(str_subset(names(wide_sensor),"OLI|MSI")),
               names_to = c("Sensor","Band","Cwl"),
               values_to = "Rrs",
               values_drop_na = T,
               names_sep = "_") %>%
    mutate(Cwl = as.numeric(Cwl))

# In-situ data for algo development
possibleApprox <- purrr::possibly(approx, otherwise=NA)

approx_op <- models_db %>%
    mutate(dataprox = purrr::map(.x = data, ~ data.frame(possibleApprox(.x$Lambda,.x$Rrs,c(380:800), na.rm = T)))) %>%
    mutate(dataprox = purrr::map(dataprox, ~ dplyr::rename(.,Lambda = x, Rrs = y))) %>%
    relocate(dataprox, .after=data)

# long_insitu_raw <- approx_op %>%
#   #mutate(data = (purrr::map(data, ~ dplyr::filter(., Lambda %in% c(440:450,540:560,650:670,725:750))))) %>%
#   select(!c("OLI","S2A","S3A","dataprox")) %>%
#   unnest(cols= c(data))
# 
# wide_insitu_raw <- long_insitu_raw %>%
#   pivot_wider(names_from = Lambda,
#             values_from = Rrs,
#             names_prefix = "Rrs_") %>%
#   relocate(matches("Rrs"), .after = Region)

long_insitu <- approx_op %>%
    mutate(dataprox = (purrr::map(dataprox, ~ dplyr::filter(., Lambda %in% c(380,395,440,560,665,710,740,800))))) %>%
    select(!c("OLI","S2A","data")) %>%
    unnest(cols= c(dataprox))

wide_insitu <- long_insitu %>%
    pivot_wider(names_from = Lambda,
              values_from = Rrs,
              names_prefix = "Rrs_") %>%
    relocate(matches("Rrs"), .after = Region) %>%
    ungroup()

# Check that aprox function dont mess up
# plot(wide_insitu_raw$Rrs_650, wide_insitu$Rrs_665)

iop_long <- wide_insitu %>%
    select(matches("^ID$|A_|Ap_|Aph_|Ag_|Anap_|Bbp_|Bb_")) %>%
    pivot_longer(cols = all_of(str_subset(names(wide_insitu), "A_|Ap_|Aph_|Ag_|Anap_|Bbp_|Bb_")),
               names_to = c(".value","Lambda"),
               names_pattern = "(.+)_(.+)",
               values_drop_na = T
               ) %>%
    mutate(Lambda = as.numeric(Lambda))



iop_nest <- iop_long %>% arrange(ID,Lambda) %>%
    group_by(ID) %>% nest()

possibleApprox <- purrr::possibly(approx, otherwise=NA)

Approxiop <- iop_nest %>%
    mutate(
        Lambda = purrr::map(.x = data, ~ data.frame(Lambda = .x$Lambda)),
        A = purrr::map(.x = data, ~ data.frame(possibleApprox(.x$Lambda,.x$A,.x$Lambda, na.rm = T))),
        Ap = purrr::map(.x = data, ~ data.frame(possibleApprox(.x$Lambda,.x$Ap,.x$Lambda, na.rm = T))),
        Aph = purrr::map(.x = data, ~ data.frame(possibleApprox(.x$Lambda,.x$Aph,.x$Lambda, na.rm = T))),
        Anap = purrr::map(.x = data, ~ data.frame(possibleApprox(.x$Lambda,.x$Anap,.x$Lambda, na.rm = T))),
        Ag = purrr::map(.x = data, ~ data.frame(possibleApprox(.x$Lambda,.x$Ag,.x$Lambda, na.rm = T))),
        Bbp = purrr::map(.x = data, ~ data.frame(possibleApprox(.x$Lambda,.x$Bbp,.x$Lambda, na.rm = T))),
        Bb = purrr::map(.x = data, ~ data.frame(possibleApprox(.x$Lambda,.x$Bb,.x$Lambda, na.rm = T))),
          ) %>%
    select(ID,Lambda,A,Ap,Aph,Anap,Ag,Bbp,Bb)

Approxiop <- Approxiop %>%
    dplyr::filter(purrr::map_lgl(A, ~ !is_empty(.)) | purrr::map_lgl(Ap, ~ !is_empty(.)) | purrr::map_lgl(Aph, ~ !is_empty(.)) | purrr::map_lgl(Anap, ~ !is_empty(.)) | purrr::map_lgl(Ag, ~ !is_empty(.)) | purrr::map_lgl(Bbp, ~ !is_empty(.)) | purrr::map_lgl(Bb, ~ !is_empty(.)))

Approxiop <- Approxiop %>% unnest(cols = c(Lambda, A, Ap, Aph, Anap, Ag, Bbp, Bb), names_sep="_") # %>% na.omit()

Approxiop <- Approxiop %>% select(A_x,!contains("x")) %>% mutate(
    Lambda = as.numeric(Lambda_Lambda),
    A = A_y,
    Ap = Ap_y,
    Aph = Aph_y,
    Anap = Anap_y,
    Ag = Ag_y,
    Bbp = Bbp_y,
    Bb = Bb_y
    ) %>%
    select(!contains(c("y","x","Lambda_Lambda"))) %>%
    filter(Lambda %in%
            c(305,320,330,340,380,394,395,412,420,443,465,470,490,510,532,555,560,589,620,625,665,683,694,700,710,740,780,875))

# Rrs at 9 wl give 255 731 row ... rather slow but working
# iopgate <- long_sensorBand %>% filter(Sensor=="MSI") %>%
#   select(matches("ID|ROI|Sensor|Band|Rrs|SPM|Ag_440")) %>%
#   right_join(Approxiop, by="ID") %>%
#   na.omit() %>% ungroup()


iopgate <- wide_insitu %>%
    select(ID,PID,Project,Station,Zstation,DateTime,Comment,SPM,QC_SPM,QC_Ag) %>%
    right_join(Approxiop, by="ID")

sd_iops <- SharedData$new(iopgate %>% group_by(ID), key = ~str_pad(ID,3,pad = "0"), group = "iops")

submap <- wide_insitu %>% mutate(DateTime = as.POSIXct(DateTime, tz='UTC', format = "%Y-%m-%dT%H:%M:%S")) %>%
    select(ID,PID,Project,Station,DateTime,SPM,Zstation,Zsecchi,Lat,Lon,Ag_440,SID) %>%
    group_by(ID)
sd_map <- SharedData$new(submap, key = ~str_pad(ID,3,pad = "0"), group = "iops")

subrrs <- long_insitu %>% #dplyr::filter(Sensor=="MSI") %>%
    select(ID,PID,Station,Zstation,Zsecchi,DateTime,Comment,Lambda,Rrs,SPM,Ag_440,QC_SPM,QC_Ag) %>%
    group_by(ID)
sd_rrs <- SharedData$new(subrrs, key = ~str_pad(ID,3,pad = "0"), group = "iops")

subtab <- wide_insitu %>%
    mutate(ID = str_pad(ID,3,pad = "0"),
          DateTime = as.POSIXct(DateTime, tz='UTC', format = "%Y-%m-%dT%H:%M:%S")) %>%
    select(ID,SID,PID,SPID,Station,Zstation,Zsecchi,DateTime,QC_SPM,QC_Ag,Comment) %>% group_by(ID)

sd_tab <- SharedData$new(subtab, key = ~str_pad(ID,3,pad = "0"), group = "iops")

Row {.sidebar}

filters

filter_checkbox(id = "Project",
              label = "Project",
              sharedData = sd_map,
            group = ~Project,
            inline = T
            )

filter_select(id = "ID",
              label = "ID",
              sharedData = sd_map,
            group = ~ID
            )

filter_select(id = "PID",
              label = "PID",
              sharedData = sd_map,
            group = ~PID
            )

filter_select(id = "Station",
              label = "Station",
              sharedData = sd_map,
            group = ~Station
            )
# not working ... some bug in display
# filter_slider(id = "DateTime",
#               label = "",
#               sharedData = sd_map,
#           column = ~DateTime, step = NULL, width = '100%', dragRange = TRUE
#           )

filter_select(id = "DateTime",
              label = "Date",
              sharedData = sd_map,
            group = ~lubridate::date(DateTime)
            )

filter_slider(id = "Depth",
              label = "Station Depth",
              sharedData = sd_map,
            column = ~Zstation
            )

filter_slider(id = "SPM",
              label = "SPM concentration",
              sharedData = sd_map,
            column = ~signif(SPM,2)
            )

filter_slider(id = "Ag_440",
              label = "CDOM concentration",
              sharedData = sd_map,
            column = ~signif(Ag_440,2)
            )

Column {data-width=450}

map

# Temp <- CTD_sd %>% plot_ly(x=~Temp, y=~Depth, showlegend=T, connectgaps = F) %>% add_paths() %>%
#   layout(yaxis = Noax, xaxis = list(title = "\u00B0C"))
# 
# PSU <- CTD_sd %>% plot_ly(x=~PSU, y=~Depth, showlegend=T, connectgaps = F) %>% add_paths() %>%
#   layout(yaxis = Noax, xaxis = list(title = "PSU")) 
# 
# CTD <- subplot(Temp, PSU, shareY= T, titleY = T, titleX = T) %>%
#   highlight(on = "plotly_click", off = "plotly_doubleclick", selectize = F, dynamic = F, persistent = F)
# 
# CTD <- bscols()

#leafpop::popupGraph should be the way to add a plot_ly graph as leaflet popup

link <- "/mnt/D/Data/WISEMan/L2/066/COPS_Kildir/PDF/WISE_CAST_006_190821_143528_URC.csv.pdf"

map <- leaflet(sd_map) %>%
    addScaleBar("bottomright") %>% leafem::addMouseCoordinates() %>%
    addProviderTiles(provider = providers$CartoDB.Positron, group = 'Positron') %>%
    addProviderTiles("Esri.WorldImagery", group = 'Aerial') %>%
    addProviderTiles("OpenTopoMap", group = 'Terrain') %>%
    addAwesomeMarkers(lat = ~Lat, lng = ~Lon, 
                  group = "Stations",
                  label = ~paste0("PID: ", PID),
                  popup = ~paste0('<h3>Station Details</h3><br>',
                                   '<b>ID</b>: ', ID, '<br>',
                            '<b>SID</b>: ', str_c(SID, collapse = ", "), '<br>',
                            '<b>PID</b>: ', PID, '<br>',
                            '<b>Station</b>: ', Station, '<br>',
                            '<b>DateTime</b>: ', DateTime, '<br>',
                            '<b>SPM</b>: ', SPM, ' [mg.l-1]<br>',
                            '<b>Zstation</b>: ', Zstation, ' [m-1]<br>',
                            '<b>Zsecchi</b>: ', Zsecchi, ' [m-1]<br>',
                            paste0('<a href="',link,'">COPS Report:</a>')
                            )
                  ) %>%
  addLayersControl(
        baseGroups = c("Positron", "Aerial", "Terrain"),
        overlayGroups = 'Stations',
        options = layersControlOptions(collapsed = TRUE)
      ) %>%
  addEasyButton(easyButton(
    icon="fa-globe", title="Zoom to Level 1",
    onClick=JS("function(btn, map){ map.setZoom(1); }"))) %>%
  addEasyButton(easyButton(
    icon="fa-crosshairs", title="Locate Me",
    onClick=JS("function(btn, map){ map.locate({setView: true}); }"))) %>% 
  addSearchFeatures(targetGroups = c("Stations"))

map

#bscols(widths = c(2, NA) ,Qfilter ,map)

# Cannot filter and cluster at the same time for now : https://github.com/rstudio/leaflet/issues/478
# clusterOptions = markerClusterOptions(disableClusteringAtZoom = 10)

data table

datatable(sd_tab,
  extensions = c("Buttons", "Scroller", "Select"),
  filter = "top",
  escape = TRUE, rownames = FALSE,
  style = "bootstrap",
  class = "compact",
  height = 100,
  options = list(
    dom = "Brtip",
    select = list(style = 'os', items = 'row'),
    buttons = list(I("colvis"),"selectNone","csv"),
    columnDefs = list(
        list(
          visible = FALSE,
          targets = c(1,3,5,6,8,9)
        )),
    deferRender = TRUE,
    scrollY = 10,
    pageLength = 10,
    scroller = TRUE
  ),
  selection = "none"
)

Column {.tabset}

Rrs

prrs <- plot_ly(sd_rrs, x = ~Lambda, y = ~Rrs,
                      text=~paste0(
                            '<b>ID</b>: ', ID, '<br>',
                            '<b>PID</b>: ', PID, '<br>',
                            '<b>Station</b>: ', Station, '<br>',
                            '<b>DateTime</b>: ', DateTime, '<br>',
                            '<b>SPM</b>: ', SPM, ' [mg.L<sup>-1</sup>]<br>',
                            '<b>a<sub>g</sub>(440)</b>: ', Ag_440, ' [m<sup>-1</sup>]<br>',
                            '<b>Zstation</b>: ', Zstation, ' [m-1]<br>',
                            '<b>Zsecchi</b>: ', Zsecchi, ' [m-1]<br>',
                            '<b>QC</b>: ', QC_SPM, '<br>',
                            '<b>Comment</b>: ', Comment, '<br>'
                            )) %>%
    add_lines(color = ~as.character(QC_SPM), colors = c("0" = "red", "1" = "steelblue3", "2"="orange"), showlegend= T) %>%
    layout(
        shapes = list(
            type = "rect", fillcolor = "transparent", xref = "paper", yref = "paper", x0 = 0, x1 = 1, y0 = 0.01, y1 = 1),
        xaxis = list(title = TeX("wavelength\\ [nm]"), type = ""),
        yaxis = list(title = TeX("R_{rs}\\ [sr^{-1}]"), type = "")
        ) %>%
    highlight(on = "plotly_click", off = "plotly_relayout", selectize = F, dynamic = F, persistent = F) %>%
    config(mathjax = "cdn", displayModeBar = T)

prrs

IOPs

piops <- plot_ly(sd_iops, x = ~Lambda, visible = "legendonly",
              # linetype = ~as.character(QC_SPM), linetypes = c("0" = "dash", "1" = "solid"), showlegend = T,

              text=~paste0(
                  '<b>PID</b>: ', PID, '<br>',
                  '<b>Station</b>: ', Station, '<br>',
                  '<b>DateTime</b>: ', DateTime, '<br>',
                  '<b>SPM</b>: ', SPM, ' [mg.l-1]<br>',
                  '<b>Zstation</b>: ', Zstation, ' [m-1]<br>',
                  '<b>QC</b>: ', QC_Ag, '<br>',
                  '<b>Comment</b>: ', Comment, '<br>'
                  )) %>%
    add_lines(y = ~A , color="A", showlegend = T) %>%
    add_lines(y = ~Ap , color="Ap", showlegend = T) %>%
    add_lines(y = ~Aph , color="Aph", showlegend = T) %>%
    add_lines(y = ~Anap , color="Anap", showlegend = T) %>%
    add_lines(y = ~Ag , color="Ag", showlegend = T) %>%
    add_lines(y = ~Bb , color="Bb", showlegend = T) %>%
    add_lines(y = ~Bbp , color="Bbp", showlegend = T)  %>%
     layout(
        shapes = list(
            type = "rect", fillcolor = "transparent", xref = "paper", yref = "paper", x0 = 0, x1 = 1, y0 = 0.01, y1 = 1),
        xaxis = list(title = TeX("wavelength\\ [nm]"), type = ""),
        yaxis = list(title = TeX("IOP\\ coefficient\\ [m^{-1}]"), type = "")
        ) %>%
    highlight(on = "plotly_click", off = "plotly_relayout", selectize = F, dynamic = F, persistent = F) %>%
    config(mathjax = "cdn", displayModeBar = T)

piops

raphidoc/lighthouse documentation built on June 13, 2022, 10:06 a.m.