Merge pull request #51 from noaa-afsc/revisions

Revisions in Response to Peer Review

.devcontainer/Dockerfile

@@ -0,0 +1,13 @@
+FROM rocker/geospatial
+
+RUN R -e "install.packages('renv', repos = c(CRAN = 'https://cloud.r-project.org'))"
+
+WORKDIR /project
+COPY renv.lock renv.lock
+
+RUN mkdir -p renv
+COPY .Rprofile .Rprofile
+COPY renv/activate.R renv/activate.R
+COPY renv/settings.json renv/settings.json
+
+RUN R -e "renv::restore()"

.devcontainer/devcontainer.json

@@ -0,0 +1,6 @@
+{
+  "name": "${localWorkspaceFolderBasename}",
+  "image": "ghcr.io/rocker-org/devcontainer/geospatial:4",
+
+}
+

.gitignore

@@ -6,3 +6,8 @@
 _targets
 renv
 *.Rproj
+.secrets
+data_raw/era5_ssrd
+/Users/joshlondon/dev/berchukseals-haulout/.secrets
+
+/.luarc.json

R/add_ssrd.R

@@ -0,0 +1,35 @@
+add_ssrd <- function(d) {
+  years <- 2005:2021
+
+  raster_merge <- vector(mode = "list", length = length(years))
+
+  for (i in 1:length(years)) {
+    left_file <- here::here('data_raw','era5_ssrd',paste0('era5_ssrd_',years[i],'_left.nc'))
+    right_file <- here::here('data_raw','era5_ssrd',paste0('era5_ssrd_',years[i],'_right.nc'))
+    r1 <- terra::rast(left_file)
+    r2 <- terra::rast(right_file)
+    raster_merge[[i]] <- terra::merge(r1,r2)
+  }
+
+  ssrd_dataset <- terra::rast(raster_merge)
+
+  pts_vec <- project(terra::vect(d), 'EPSG:4326')
+
+  # list to store for-loop results
+  pts_ssrd_list <- vector(mode = "list", length = length(years))
+
+  for (i in 1:length(years)) {
+    r <- subset(ssrd_dataset, lubridate::year(time(ssrd_dataset)) == years[i])
+    v <- subset(pts_vec, lubridate::year(pts_vec$haulout_dt) == years[i])
+    tm <- match(v$haulout_dt,time(r))
+    pts_ssrd_list[[i]] <- terra::extract(r,v,layer=tm, bind=TRUE)
+  }
+
+  ssrd_data <- do.call(rbind,pts_ssrd_list) %>%
+    sf::st_as_sf() %>% 
+    dplyr::rename(era5_ssrd = value) %>% 
+    dplyr::mutate(era5_ssrd_watts = era5_ssrd/3600) %>%
+    dplyr::select(-layer)
+
+  return(ssrd_data)
+}

R/adfg_clean.R

@@ -18,7 +18,7 @@ adfg_clean_deploy <- function(file1,file2) {
                                                    'America/Anchorage')) %>%
     dplyr::mutate(DateTagged = lubridate::as_datetime(DateTagged,
                                                       tz = 'America/Anchorage')) %>%
-    dplyr::filter( ! (DeployIDs %in% c("EB16SMK-158429")) )
+    dplyr::filter( ! (DeployIDs %in% c("EB16SMK-158429")) ) #included in adfg file2
 
   cols <- readr::cols(
     DeployIDs = readr::col_character(),
@@ -48,7 +48,7 @@ adfg_clean_deploy <- function(file1,file2) {
     dplyr::mutate(species = tolower(species),
                   age = tolower(age),
                   age = case_when(
-                    age %in% c("1","3") ~ "SUBADULT",
+                    age %in% c("1","3") ~ "subadult",
                     TRUE ~ age
                   ),
                   sex = tolower(sex)) %>%
@@ -82,7 +82,7 @@ adfg_clean_locs <- function(file1, file2, adfg_deployments) {
     dplyr::mutate(Date = as.Date(Date, origin = "1899-12-30", tz ="UTC") %>%
                     lubridate::as_datetime(),
                   Ptt = as.integer(Ptt)) %>%
-    dplyr::filter( ! (DeployID %in% c("EB16SMK-158429")) )
+    dplyr::filter( ! (DeployID %in% c("EB16SMK-158429")) ) # included in adfg file2
 
   cols <- readr::cols(
     DeployID = readr::col_character(),
@@ -162,6 +162,11 @@ adfg_clean_tl <- function(file1,file2, adfg_deployments) {
                   hist_type = HistType,
                   timeline_start_dt = GMTDate) %>%
     dplyr::rename_all(tolower) %>%
+    dplyr::mutate(tag_family = case_when(
+      ptt %in% c(110597, 158422) ~ "SPLASH",
+      ptt %in% c(167946, 137526) ~ "SPOT",
+      TRUE ~ tag_family
+    )) %>% 
     dplyr::left_join(adfg_deployments, by = c("speno","tag_family")) %>%
     mutate(species_code = case_when(
       species == "bearded" ~ "EB",

R/create_ar1_id.R

@@ -10,3 +10,15 @@ create_ar1_id <- function(h_dt) {
 
   return(ar1_id)
 }
+
+create_ar1_start <- function(ar1_id) {
+  rle_ar1 = ar1_id %>% rle()
+  end = cumsum(rle_ar1$lengths)
+  start = c(1, lag(end)[-1] + 1)
+
+  ar1_start <- vector(length = length(ar1_id))
+  ar1_start[seq_along(ar1_start)] <- FALSE
+  ar1_start[start] <- TRUE
+
+  return(ar1_start)
+}

R/create_data_sf.R

@@ -34,9 +34,11 @@ create_data_sf <- function(locs_sf, source_data) {
     mutate(age = forcats::fct_relevel(age, c("ADULT", "SUBADULT", "YOUNG OF YEAR"))) %>%
     dplyr::select(-c(bday,age_days)) %>%
     relocate(deploy_dt, .after = age) %>%
-    filter(!(species == "Bearded seal" & age == "YOUNG OF YEAR"))
+    filter(!(species == "Bearded seal" & age == "YOUNG OF YEAR"),
+           !(speno %in% c("EB2009_7010","EB2009_3002")) #Sea of Okhotsk bearded seals
+           ) 
 
-  # read in world land polygons and buffer in by 100m; want to ensure only locations well
+  # read in world land polygons and buffer in by 100m; want to ensure only locations well inland filtered out
   world <- ne_countries(scale = "large", returnclass = "sf") %>%
     filter(continent != "Antarctica") %>%
     st_wrap_dateline() %>% st_transform(st_crs(dat.sf)) %>%

R/create_db_input_data.R

@@ -31,37 +31,51 @@ create_source_data <- function(locs_sf, timeline_data) {
                      y = weighted.mean(y,1/error_radius))
 
   tbl_percent_locs <- timeline_data %>%
+    ungroup() %>%
+    # let's make sure we only have hourly summarized timeline data
+    group_by(speno,species,sex,age,unique_day,n_tags,
+             timeline_hour = lubridate::hour(timeline_start_dt)) %>%
+    reframe( #using reframe here b/c summarize creates a warning about more (or less) than 1 row per `summarise()`
+      timeline_start_dt = lubridate::floor_date(timeline_start_dt, "hours"),
+      percent_dry = mean(percent_dry)
+    ) %>%
+    ungroup() %>%
+    dplyr::select(-timeline_hour) %>%
+    # now join with daily locations
     full_join(locs_daily,
                      by = c("speno","unique_day","age",
                             "sex","species")) %>%
     arrange(speno,unique_day,timeline_start_dt) %>%
-    group_by(speno) %>% nest() %>%
-    mutate(start_idx = map_int(data,~ which.max(!is.na(.x$x))),
-           data = map2(data, start_idx, ~ slice(.x, .y:nrow(.x)))) %>%
+    group_by(speno) %>% 
+    tidyr::nest() %>%
+    mutate(start_idx = purrr::map_int(data,~ which.max(!is.na(.x$x))),
+           data = purrr::map2(data, start_idx, ~ slice(.x, .y:nrow(.x)))) %>%
     dplyr::select(-start_idx) %>%
-    unnest(cols = c(data)) %>%
+    tidyr::unnest(cols = c(data)) %>%
     mutate(fill_xy = ifelse(is.na(x), TRUE, FALSE)) %>%
     group_by(speno) %>%
-    fill(x,y) %>%
+    tidyr::fill(x,y) %>%
     ungroup() %>%
     dplyr::filter(!is.na(percent_dry)) %>%
-    dplyr::filter(!speno %in% c("EB2005_5995","PL2006_5984","PL2006_5987")) %>%
+    dplyr::filter(!speno %in% c("EB2005_5995")) %>% #peard bay capture; only 1 day of data
     dplyr::filter(!is.na(x)) %>%
     st_as_sf(coords = c("x","y")) %>%
     st_set_crs(3571) %>%
     rename(haulout_dt = timeline_start_dt) %>%
+
     dplyr::select(speno,species,age,sex,haulout_dt,percent_dry,n_tags,fill_xy)
 
-  stopifnot(
-    "PEP Postgres Database Not Available; did you start VPN? ;)" =
-      is_up("161.55.120.122", "5432")
-  )
 
-  con <- dbConnect(RPostgres::Postgres(),dbname = 'pep',
-                   host = '161.55.120.122',
-                   port = 5432,
-                   user = keyringr::get_kc_account("pgpep_sa"),
-                   password = keyringr::decrypt_kc_pw("pgpep_sa"))
+  tryCatch({
+    con <- dbConnect(RPostgres::Postgres(),
+                     dbname = 'pep', 
+                     host = Sys.getenv('PEP_PG_IP'),
+                     user = keyringr::get_kc_account("pgpep_sa"),
+                     password = keyringr::decrypt_kc_pw("pgpep_sa"))
+  },
+  error = function(cond) {
+    print("Unable to connect to Database.")
+  })
   on.exit(dbDisconnect(con))
 
   st_write(obj = tbl_percent_locs,
@@ -71,6 +85,7 @@ create_source_data <- function(locs_sf, timeline_data) {
   )
 
   dbExecute(con, "ALTER TABLE telem.res_iceseal_haulout RENAME COLUMN geometry TO geom")
+  dbExecute(con, "ALTER TABLE IF EXISTS telem.res_iceseal_haulout OWNER TO pep_manage_telem")
   dbExecute(con, "SELECT telem.fxn_iceseal_pred_idx();")
   dbExecute(con, "SELECT telem.fxn_iceseal_haulout_cov();")
 
@@ -80,9 +95,7 @@ create_source_data <- function(locs_sf, timeline_data) {
   WHERE
   EXTRACT(MONTH FROM haulout_dt) IN (3,4,5,6,7) AND
   rast_vwnd IS NOT NULL AND
-  species != 'Ph' AND
-  speno != 'EB2009_7010' AND
-  speno != 'EB2009_3002' "
+  species != 'Ph'"
   }
 
   sf::st_read(con, query = qry)

R/create_deploy_details.R

@@ -0,0 +1,29 @@
+create_deploy_details <- function(deploy_details_file) {
+  deploy_tbl <- readr::read_csv(deploy_details_file)
+
+  tryCatch({
+    con <- dbConnect(RPostgres::Postgres(),
+                      dbname = 'pep', 
+                      host = Sys.getenv('PEP_PG_IP'),
+                      user = keyringr::get_kc_account("pgpep_londonj"),
+                      password = keyringr::decrypt_kc_pw("pgpep_londonj"))
+  },
+  error = function(cond) {
+    print("Unable to connect to Database.")
+  })
+  on.exit(dbDisconnect(con))
+
+  deploy_db <- tbl(con, in_schema("telem","tbl_tag_deployments")) %>% 
+    dplyr::select(speno, deploy_lat, deploy_long, capture_lat, capture_long) %>% 
+    dplyr::mutate(deploy_lat = if_else(is.na(deploy_lat), capture_lat, deploy_lat),
+                  deploy_long = if_else(is.na(deploy_long), capture_long, deploy_long)) %>% 
+    collect()
+
+  dbDisconnect(con)
+
+  deploy_tbl <- deploy_tbl %>% 
+    # dplyr::select(speno, species, location, region, latitude, longitude) %>% 
+    left_join(deploy_db, by = "speno") %>% 
+    mutate(latitude = coalesce(latitude,deploy_lat),
+           longitude = coalesce(longitude, deploy_long))
+}

R/create_grid.R

@@ -1,15 +1,17 @@
 create_grid_sf <- function(data_sf) {
-  stopifnot(
-    "PEP Postgres Database Not Available; did you start VPN? ;)" =
-      is_up("161.55.120.122", "5432")
-  )
 
-  con <- dbConnect(RPostgres::Postgres(),dbname = 'pep',
-                   host = '161.55.120.122',
-                   port = 5432,
-                   user = keyringr::get_kc_account("pgpep_londonj"),
-                   password = keyringr::decrypt_kc_pw("pgpep_londonj"))
-  on.exit(dbDisconnect(con))
+  tryCatch({
+    con <- dbConnect(RPostgres::Postgres(),
+                     dbname = 'pep', 
+                     host = Sys.getenv('PEP_PG_IP'),
+                     user = keyringr::get_kc_account("pgpep_londonj"),
+                     password = keyringr::decrypt_kc_pw("pgpep_londonj"))
+  },
+  error = function(cond) {
+    print("Unable to connect to Database.")
+  })
+  on.exit(odbc::dbDisconnect(con))
+
   sf::st_read(con,
               query = "SELECT * FROM base.geo_analysis_grid") %>%
     st_filter(st_bbox(data_sf) %>% st_as_sfc())

R/create_model_input.R

@@ -1,5 +1,5 @@
 create_model_input <- function(dat.sf) {
-  dat.sf  %>%
+  temp_dat <- dat.sf  %>%
   mutate(dry = round(percent_dry/100)) %>%
   mutate(hour_utc = lubridate::hour(haulout_dt),
          yday = lubridate::yday(haulout_dt),
@@ -14,12 +14,13 @@ create_model_input <- function(dat.sf) {
          sin3 = sin(pi*solar_hour/4),
          cos3 = sin(pi*solar_hour/4),
          day = (yday-120)/10,
-         day2 = day^2,
-         day3 = day^3
+         day_2 = day^2,
+         day_3 = day^3,
+         day_4 = day^4
   ) %>%
   mutate(precip = rast_acpcp,
-         temp2 = (rast_air2m-270)/27,
-         temp0 = (rast_airsfc-270)/27,
+         temp2m = (rast_air2m-273.15)/27,
+         temp0 = (rast_airsfc-273.15)/27,
          pressure = (rast_prmsl-100000)/10000,
          wind = sqrt(rast_uwnd^2 + rast_vwnd^2)/10
   ) %>%
@@ -33,11 +34,13 @@ create_model_input <- function(dat.sf) {
          solar_hour = forcats::as_factor(solar_hour)) %>%
   arrange(speno,haulout_dt) %>%
   group_by(speno) %>%
-  mutate(ar1_id = paste(speno,create_ar1_id(haulout_dt),sep="_")) %>%
+  mutate(ar1_id = paste(speno,create_ar1_id(haulout_dt),sep="_"),
+         ar1_start = create_ar1_start(ar1_id)) %>%
   ungroup() %>%
   mutate(ar1_id = forcats::as_factor(ar1_id),
          speno = forcats::as_factor(speno)) %>%
   relocate(ar1_id, .after = speno) %>%
+  relocate(ar1_start, .after = ar1_id) %>%
   dplyr::select(-c(coords_x, coords_y)) %>%
   group_by(ar1_id) %>% mutate(time_vec = row_number(ar1_id) - 1) %>%
   ungroup() %>%
@@ -56,4 +59,6 @@ create_model_input <- function(dat.sf) {
     age_sex_inter = forcats::as_factor(age_sex_inter) %>%
       forcats::fct_relevel(c("ADULT.F","ADULT.M","SUBADULT"))
   )
+
+  return(temp_dat)
 }