From eb01cb989d1502ec25c47836be56e7764d2568da Mon Sep 17 00:00:00 2001 From: Brendan Turley Date: Thu, 6 Jun 2024 17:55:24 -0400 Subject: [PATCH 1/7] testing rerddap functionality as replacement --- indicator_processing/automated_download/sst.R | 79 +++++++++++++------ 1 file changed, 55 insertions(+), 24 deletions(-) diff --git a/indicator_processing/automated_download/sst.R b/indicator_processing/automated_download/sst.R index 603e877..bcffec8 100644 --- a/indicator_processing/automated_download/sst.R +++ b/indicator_processing/automated_download/sst.R @@ -5,8 +5,10 @@ rm(list = ls()) dev.off() +library(lubridate) library(maps) library(plotTimeSeries) +library(rerddap) load("indicator_processing/spec_file.RData") @@ -17,37 +19,66 @@ load("indicator_processing/spec_file.RData") styear <- 1982 enyear <- terminal_year +# get ERDDAP info -------------------------------- +sst <- info('ncdcOisst21Agg') +sst <- info('ncdcOisst21Agg_LonPM180') # this may work better + # empty data ------------------------------------------------- -dat <- data.frame(row.names = c("year", "mon", "PR_mean", "PR_min", "PR_max", "VI_mean", "VI_min", "VI_max")) +# dat <- data.frame(row.names = c("year", "mon", "PR_mean", "PR_min", "PR_max", "VI_mean", "VI_min", "VI_max")) + +dat <- setNames(data.frame(matrix(NA,length(styear:enyear)*12,5)), + c("year", "mon", 'mean', 'min', 'max')) +m <- 1 +n <- 0 # download by year to avoid timeout errors -------------------- for (yr in styear:enyear) { -# url from ERDDAP for OISST, download and read ---------------- - url <- paste0("https://coastwatch.pfeg.noaa.gov/erddap/griddap/ncdcOisst21Agg.csv?sst[(", yr, "-01-01T12:00:00Z):1:(", yr, "-12-31T12:00:00Z)][(0.0):1:(0.0)][(17):1:(19)][(292):1:(296)]") - - # url_vi <- paste0("https://coastwatch.pfeg.noaa.gov/erddap/griddap/ncdcOisst21Agg.csv?sst[(", yr, "-01-01T12:00:00Z):1:(", yr, "-12-31T12:00:00Z)][(0.0):1:(0.0)][(17):1:(19)][(294.75):1:(296)]") - # url_pr <- paste0("https://coastwatch.pfeg.noaa.gov/erddap/griddap/ncdcOisst21Agg.csv?sst[(", yr, "-01-01T12:00:00Z):1:(", yr, "-12-31T12:00:00Z)][(0.0):1:(0.0)][(17):1:(19)][(292):1:(294.75)]") - # download for entire Caribbean because USVI and PR highly correlated. - download.file(url = url, destfile = "st.csv") - - labs <- read.table("st.csv", sep = ",", header = T, skip = 0) - sst_vi <- read.table("st.csv", sep = ",", header = T, skip = 1) - names(sst_vi) <- names(labs) - -# extract month ----------------------------------------------- - sst_vi$mon <- strftime(sst_vi$time, format = "%m") - -# calculate mean, min, max and concatenate -------------------- - ind_vi <- as.numeric(tapply(sst_vi$sst, sst_vi$mon, mean, na.rm = T)) - min_vi <- as.numeric(tapply(sst_vi$sst, sst_vi$mon, min, na.rm = T)) - max_vi <- as.numeric(tapply(sst_vi$sst, sst_vi$mon, max, na.rm = T)) - - tempdat <- data.frame(yr, unique(sst_vi$mon), ind_vi, min_vi, max_vi, stringsAsFactors = F) - dat <- data.frame(rbind(dat, tempdat), stringsAsFactors = F) + ### BDT rERDDAP fix + sst_grab <- griddap(sst, fields = 'sst', + time = c(paste0(yr,'-01-01'), paste0(yr,'-12-31')), + longitude = c(360 + min_lon, 360 + max_lon), + latitude = c(min_lat, max_lat)) + + sst_grab <- griddap(sst, fields = 'sst', + time = c(paste0(yr,'-01-01'), paste0(yr,'-12-31')), + longitude = c(min_lon, max_lon), + latitude = c(min_lat, max_lat)) + + sst_agg <- aggregate(sst_grab$data$sst, + by = list(year(sst_grab$data$time), month(sst_grab$data$time)), + function(x) c(mean(x, na.rm = T), min(x, na.rm = T), max(x, na.rm = T))) + + n <- n + 12 + dat[m:n,] <- data.frame(sst_agg[,1:2], unlist(sst_agg[,3])) + m <- n + 1 + +# +# # url from ERDDAP for OISST, download and read ---------------- +# url <- paste0("https://coastwatch.pfeg.noaa.gov/erddap/griddap/ncdcOisst21Agg.csv?sst[(", yr, "-01-01T12:00:00Z):1:(", yr, "-12-31T12:00:00Z)][(0.0):1:(0.0)][(17):1:(19)][(292):1:(296)]") +# +# # url_vi <- paste0("https://coastwatch.pfeg.noaa.gov/erddap/griddap/ncdcOisst21Agg.csv?sst[(", yr, "-01-01T12:00:00Z):1:(", yr, "-12-31T12:00:00Z)][(0.0):1:(0.0)][(17):1:(19)][(294.75):1:(296)]") +# # url_pr <- paste0("https://coastwatch.pfeg.noaa.gov/erddap/griddap/ncdcOisst21Agg.csv?sst[(", yr, "-01-01T12:00:00Z):1:(", yr, "-12-31T12:00:00Z)][(0.0):1:(0.0)][(17):1:(19)][(292):1:(294.75)]") +# # download for entire Caribbean because USVI and PR highly correlated. +# download.file(url = url, destfile = "st.csv") +# +# labs <- read.table("st.csv", sep = ",", header = T, skip = 0) +# sst_vi <- read.table("st.csv", sep = ",", header = T, skip = 1) +# names(sst_vi) <- names(labs) +# +# # extract month ----------------------------------------------- +# sst_vi$mon <- strftime(sst_vi$time, format = "%m") +# +# # calculate mean, min, max and concatenate -------------------- +# ind_vi <- as.numeric(tapply(sst_vi$sst, sst_vi$mon, mean, na.rm = T)) +# min_vi <- as.numeric(tapply(sst_vi$sst, sst_vi$mon, min, na.rm = T)) +# max_vi <- as.numeric(tapply(sst_vi$sst, sst_vi$mon, max, na.rm = T)) +# +# tempdat <- data.frame(yr, unique(sst_vi$mon), ind_vi, min_vi, max_vi, stringsAsFactors = F) +# dat <- data.frame(rbind(dat, tempdat), stringsAsFactors = F) } -file.remove("st.csv") +# file.remove("st.csv") # add row names and yearmonth column -------------------------- names(dat) <- c("year", "mon", "mean", "min", "max") From 100a95a231402819fbf42bd7655b65360b86c382 Mon Sep 17 00:00:00 2001 From: Brendan Turley Date: Fri, 7 Jun 2024 10:48:10 -0400 Subject: [PATCH 2/7] fix ERDDAP timeout issue using rerddap package --- .Rhistory | 512 ------------------ .gitignore | 2 + indicator_objects/Carib_SST.RData | Bin 0 -> 8006 bytes indicator_processing/automated_download/sst.R | 11 +- 4 files changed, 4 insertions(+), 521 deletions(-) delete mode 100644 .Rhistory create mode 100644 indicator_objects/Carib_SST.RData diff --git a/.Rhistory b/.Rhistory deleted file mode 100644 index c53d1a8..0000000 --- a/.Rhistory +++ /dev/null @@ -1,512 +0,0 @@ -class(inddata) <- "indicatordata" -save(inddata, file = "indicator_objects/RVC_STX.RData") -load("../indicator_objects/RVC_PR.RData") -plotIndicatorTimeSeries(inddata, coltoplot = 1:6, trendAnalysis = T, sublabel = T) -# Puerto Rico -datdata <- 2001:2021 -inddata <- data.frame(cbind(RUVdensity_PRICO_BAL_VETU$species_data.density, RUVdensity_PRICO_EPI_GUTT$species_data.density, RUVdensity_PRICO_LUT_ANAL$species_data.density, RUVdensity_PRICO_OCY_CHRY$species_data.density, RUVdensity_PRICO_SPA_AURO$species_data.density, RUVdensity_PRICO_SPA_VIRI$species_data.density)) -labs <- c("PR" , "Density", "queen triggerfish", -"PR" , "Density", "red hind", -"PR" , "Density", "mutton snapper", -"PR" , "Density", "yellowtail snapper", -"PR" , "Density", "redband parrotfish", -"PR" , "Density", "stoplight parrotfish") -indnames <- data.frame(matrix(labs, nrow = 6, byrow = F)) -inddata <- list(labels = indnames, indicators = inddata, datelist = datdata) -class(inddata) <- "indicatordata" -save(inddata, file = "indicator_objects/RVC_PR.RData") -View(RUVdensity_PRICO_SPA_VIRI) -View(RUVdensity_PRICO_SPA_VIRI) -load("../indicator_objects/RVC_PR.RData") -plotIndicatorTimeSeries(inddata, coltoplot = 1:6, trendAnalysis = T, sublabel = T) -load("../indicator_objects/landings.RData") -load("../indicator_objects/RVC_PR.RData") -load("../indicator_objects/landings.RData") -load("../indicator_objects/fish_density.RData") -# Puerto Rico -datdata <- 2001:2021 -inddata <- data.frame(cbind(RUVdensity_PRICO_BAL_VETU$species_data.density, RUVdensity_PRICO_EPI_GUTT$species_data.density, RUVdensity_PRICO_LUT_ANAL$species_data.density, RUVdensity_PRICO_OCY_CHRY$species_data.density, RUVdensity_PRICO_SPA_AURO$species_data.density, RUVdensity_PRICO_SPA_VIRI$species_data.density)) -labs <- c("PR" , "Density", "queen triggerfish", -"PR" , "Density", "red hind", -"PR" , "Density", "mutton snapper", -"PR" , "Density", "yellowtail snapper", -"PR" , "Density", "redband parrotfish", -"PR" , "Density", "stoplight parrotfish") -indnames <- data.frame(matrix(labs, nrow = 3, byrow = F)) -inddata <- list(labels = indnames, indicators = inddata, datelist = datdata) -class(inddata) <- "indicatordata" -save(inddata, file = "indicator_objects/RVC_PR.RData") -# St. Thomas & St. John -datdata <- 2001:2021 -inddata <- data.frame(cbind(RUVdensity_STTSTJ_BAL_VETU$species_data.density, RUVdensity_STTSTJ_EPI_GUTT$species_data.density, RUVdensity_STTSTJ_LUT_ANAL$species_data.density, RUVdensity_STTSTJ_OCY_CHRY$species_data.density, RUVdensity_STTSTJ_SPA_AURO$species_data.density, RUVdensity_STTSTJ_SPA_VIRI$species_data.density)) -labs <- c("STTSTJ" , "Density", "queen triggerfish", -"STTSTJ" , "Density", "red hind", -"STTSTJ" , "Density", "mutton snapper", -"STTSTJ" , "Density", "yellowtail snapper", -"STTSTJ" , "Density", "redband parrotfish", -"STTSTJ" , "Density", "stoplight parrotfish") -indnames <- data.frame(matrix(labs, nrow = 3, byrow = F)) -inddata <- list(labels = indnames, indicators = inddata, datelist = datdata) -class(inddata) <- "indicatordata" -save(inddata, file = "indicator_objects/RVC_STSJ.RData") -# St. Croix -datdata <- 2001:2021 -inddata <- data.frame(cbind(RUVdensity_STX_BAL_VETU$species_data.density, RUVdensity_STX_EPI_GUTT$species_data.density, RUVdensity_STX_LUT_ANAL$species_data.density, RUVdensity_STX_OCY_CHRY$species_data.density, RUVdensity_STX_SPA_AURO$species_data.density, RUVdensity_STX_SPA_VIRI$species_data.density)) -labs <- c("STX" , "Density", "queen triggerfish", -"STX" , "Density", "red hind", -"STX" , "Density", "mutton snapper", -"STX" , "Density", "yellowtail snapper", -"STX" , "Density", "redband parrotfish", -"STX" , "Density", "stoplight parrotfish") -indnames <- data.frame(matrix(labs, nrow = 3, byrow = F)) -inddata <- list(labels = indnames, indicators = inddata, datelist = datdata) -class(inddata) <- "indicatordata" -save(inddata, file = "indicator_objects/RVC_STX.RData") -load("../indicator_objects/RVC_PR.RData") -plotIndicatorTimeSeries(inddata, coltoplot = 1:6, trendAnalysis = T, sublabel = T) -# St. Croix -datdata <- 2001:2021 -inddata <- data.frame(cbind(RUVdensity_STX_BAL_VETU$species_data.density, RUVdensity_STX_EPI_GUTT$species_data.density, RUVdensity_STX_LUT_ANAL$species_data.density, RUVdensity_STX_OCY_CHRY$species_data.density, RUVdensity_STX_SPA_AURO$species_data.density, RUVdensity_STX_SPA_VIRI$species_data.density)) -labs <- c("STX" , "Density", "queen triggerfish", -"STX" , "Density", "red hind", -"STX" , "Density", "mutton snapper", -"STX" , "Density", "yellowtail snapper", -"STX" , "Density", "redband parrotfish", -"STX" , "Density", "stoplight parrotfish") -indnames <- data.frame(matrix(labs, nrow = 3, byrow = F)) -inddata <- list(labels = indnames, indicators = inddata, datelist = datdata) -class(inddata) <- "indicatordata" -load("../indicator_objects/RVC_STX.RData") -View(RUVdensity_PRICO_BAL_VETU) -# Puerto Rico -datdata <- RUVdensity_PRICO_BAL_VETU$datdata -inddata <- data.frame(cbind(RUVdensity_PRICO_BAL_VETU$species_data.density, RUVdensity_PRICO_EPI_GUTT$species_data.density, RUVdensity_PRICO_LUT_ANAL$species_data.density, RUVdensity_PRICO_OCY_CHRY$species_data.density, RUVdensity_PRICO_SPA_AURO$species_data.density, RUVdensity_PRICO_SPA_VIRI$species_data.density)) -labs <- c("PR" , "Density", "queen triggerfish", -"PR" , "Density", "red hind", -"PR" , "Density", "mutton snapper", -"PR" , "Density", "yellowtail snapper", -"PR" , "Density", "redband parrotfish", -"PR" , "Density", "stoplight parrotfish") -indnames <- data.frame(matrix(labs, nrow = 3, byrow = F)) -inddata <- list(labels = indnames, indicators = inddata, datelist = datdata) -class(inddata) <- "indicatordata" -save(inddata, file = "indicator_objects/RVC_PR.RData") -# St. Thomas & St. John -datdata <- RUVdensity_PRICO_BAL_VETU$datdata -inddata <- data.frame(cbind(RUVdensity_STTSTJ_BAL_VETU$species_data.density, RUVdensity_STTSTJ_EPI_GUTT$species_data.density, RUVdensity_STTSTJ_LUT_ANAL$species_data.density, RUVdensity_STTSTJ_OCY_CHRY$species_data.density, RUVdensity_STTSTJ_SPA_AURO$species_data.density, RUVdensity_STTSTJ_SPA_VIRI$species_data.density)) -labs <- c("STTSTJ" , "Density", "queen triggerfish", -"STTSTJ" , "Density", "red hind", -"STTSTJ" , "Density", "mutton snapper", -"STTSTJ" , "Density", "yellowtail snapper", -"STTSTJ" , "Density", "redband parrotfish", -"STTSTJ" , "Density", "stoplight parrotfish") -indnames <- data.frame(matrix(labs, nrow = 3, byrow = F)) -inddata <- list(labels = indnames, indicators = inddata, datelist = datdata) -class(inddata) <- "indicatordata" -save(inddata, file = "indicator_objects/RVC_STSJ.RData") -# St. Croix -datdata <- RUVdensity_PRICO_BAL_VETU$datdata -inddata <- data.frame(cbind(RUVdensity_STX_BAL_VETU$species_data.density, RUVdensity_STX_EPI_GUTT$species_data.density, RUVdensity_STX_LUT_ANAL$species_data.density, RUVdensity_STX_OCY_CHRY$species_data.density, RUVdensity_STX_SPA_AURO$species_data.density, RUVdensity_STX_SPA_VIRI$species_data.density)) -labs <- c("STX" , "Density", "queen triggerfish", -"STX" , "Density", "red hind", -"STX" , "Density", "mutton snapper", -"STX" , "Density", "yellowtail snapper", -"STX" , "Density", "redband parrotfish", -"STX" , "Density", "stoplight parrotfish") -indnames <- data.frame(matrix(labs, nrow = 3, byrow = F)) -inddata <- list(labels = indnames, indicators = inddata, datelist = datdata) -class(inddata) <- "indicatordata" -save(inddata, file = "indicator_objects/RVC_STX.RData") -load("../indicator_objects/RVC_PR.RData") -plotIndicatorTimeSeries(inddata, coltoplot = 1:6, trendAnalysis = T, sublabel = T) -plotIndicatorTimeSeries(inddata, coltoplot = 1:6, trendAnalysis = T, sublabel = T, widadj = 10) -load("../indicator_objects/landings.RData") -plotIndicatorTimeSeries(inddata, coltoplot = 1:6, trendAnalysis = T, sublabel = T) -# Puerto Rico -datdata <- as.integer(RUVdensity_PRICO_BAL_VETU$datdata) -inddata <- data.frame(cbind(RUVdensity_PRICO_BAL_VETU$species_data.density, RUVdensity_PRICO_EPI_GUTT$species_data.density, RUVdensity_PRICO_LUT_ANAL$species_data.density, RUVdensity_PRICO_OCY_CHRY$species_data.density, RUVdensity_PRICO_SPA_AURO$species_data.density, RUVdensity_PRICO_SPA_VIRI$species_data.density)) -labs <- c("PR" , "Density", "queen triggerfish", -"PR" , "Density", "red hind", -"PR" , "Density", "mutton snapper", -"PR" , "Density", "yellowtail snapper", -"PR" , "Density", "redband parrotfish", -"PR" , "Density", "stoplight parrotfish") -indnames <- data.frame(matrix(labs, nrow = 3, byrow = F)) -inddata <- list(labels = indnames, indicators = inddata, datelist = datdata) -class(inddata) <- "indicatordata" -save(inddata, file = "indicator_objects/RVC_PR.RData") -load("../indicator_objects/RVC_PR.RData") -plotIndicatorTimeSeries(inddata, coltoplot = 1:6, trendAnalysis = T, sublabel = T, widadj = 10) -plotIndicatorTimeSeries(inddata, coltoplot = 1:6, trendAnalysis = T, sublabel = T) -# Puerto Rico -datdata <- as.integer(RUVdensity_PRICO_BAL_VETU$datdata) -inddata <- data.frame(cbind(RUVdensity_PRICO_BAL_VETU$species_data.density, RUVdensity_PRICO_EPI_GUTT$species_data.density, RUVdensity_PRICO_LUT_ANAL$species_data.density, RUVdensity_PRICO_OCY_CHRY$species_data.density, RUVdensity_PRICO_SPA_AURO$species_data.density, RUVdensity_PRICO_SPA_VIRI$species_data.density)) -labs <- c("Puerto Rico" , "Density", "queen triggerfish", -"Puerto Rico" , "Density", "red hind", -"Puerto Rico" , "Density", "mutton snapper", -"Puerto Rico" , "Density", "yellowtail snapper", -"Puerto Rico" , "Density", "redband parrotfish", -"Puerto Rico" , "Density", "stoplight parrotfish") -indnames <- data.frame(matrix(labs, nrow = 3, byrow = F)) -inddata <- list(labels = indnames, indicators = inddata, datelist = datdata) -class(inddata) <- "indicatordata" -save(inddata, file = "indicator_objects/RVC_PR.RData") -# St. Thomas & St. John -datdata <- as.integer(RUVdensity_PRICO_BAL_VETU$datdata) -inddata <- data.frame(cbind(RUVdensity_STTSTJ_BAL_VETU$species_data.density, RUVdensity_STTSTJ_EPI_GUTT$species_data.density, RUVdensity_STTSTJ_LUT_ANAL$species_data.density, RUVdensity_STTSTJ_OCY_CHRY$species_data.density, RUVdensity_STTSTJ_SPA_AURO$species_data.density, RUVdensity_STTSTJ_SPA_VIRI$species_data.density)) -labs <- c("St. Thomas & St. John" , "Density", "queen triggerfish", -"St. Thomas & St. John" , "Density", "red hind", -"St. Thomas & St. John" , "Density", "mutton snapper", -"St. Thomas & St. John" , "Density", "yellowtail snapper", -"St. Thomas & St. John" , "Density", "redband parrotfish", -"St. Thomas & St. John" , "Density", "stoplight parrotfish") -indnames <- data.frame(matrix(labs, nrow = 3, byrow = F)) -inddata <- list(labels = indnames, indicators = inddata, datelist = datdata) -class(inddata) <- "indicatordata" -save(inddata, file = "indicator_objects/RVC_STSJ.RData") -# St. Croix -datdata <- as.integer(RUVdensity_PRICO_BAL_VETU$datdata) -inddata <- data.frame(cbind(RUVdensity_STX_BAL_VETU$species_data.density, RUVdensity_STX_EPI_GUTT$species_data.density, RUVdensity_STX_LUT_ANAL$species_data.density, RUVdensity_STX_OCY_CHRY$species_data.density, RUVdensity_STX_SPA_AURO$species_data.density, RUVdensity_STX_SPA_VIRI$species_data.density)) -labs <- c("St. Croix" , "Density", "queen triggerfish", -"St. Croix" , "Density", "red hind", -"St. Croix" , "Density", "mutton snapper", -"St. Croix" , "Density", "yellowtail snapper", -"St. Croix" , "Density", "redband parrotfish", -"St. Croix" , "Density", "stoplight parrotfish") -indnames <- data.frame(matrix(labs, nrow = 3, byrow = F)) -inddata <- list(labels = indnames, indicators = inddata, datelist = datdata) -class(inddata) <- "indicatordata" -save(inddata, file = "indicator_objects/RVC_STX.RData") -load("../indicator_objects/RVC_PR.RData") -plotIndicatorTimeSeries(inddata, coltoplot = 1:6, trendAnalysis = T, sublabel = T) -load("../indicator_objects/RVC_STSJ.RData") -plotIndicatorTimeSeries(inddata, coltoplot = 1:6, trendAnalysis = T, sublabel = T) -load("../indicator_objects/RVC_STX.RData") -plotIndicatorTimeSeries(inddata, coltoplot = 1:6, trendAnalysis = T, sublabel = T) -# The data portal has more information on the data: https://grunt.sefsc.noaa.gov/rvc_analysis20/ -# Species of interest: -# OCY CHRY --> yellowtail snapper -# LUT ANAL --> mutton snapper* -# BAL VETU --> queen triggerfish -# EPI GUTT --> red hind -# SPA AURO --> redband parrotfish -# SPA VIRI --> stoplight parrotfish -# Regions of interest: -# PRICO --> Puerto Rico -# STTSTJ --> St. Thomas & St. John -# STX --> St. Croix -# NOTE: as of 3/1/2024 calibrations have not been completed for all species. This means the data on the portal only go back to 2016. To get the full time series (2001 onward) we needed to get the calibrated data from Jeremiah Blondeau (jeremiah.blondeau@noaa.gov). Species specific calibrations were done. The calibrated_species.csv file is a list of all the species that have thus far been calibrated. All species of interest have been calibrated except for mutton snapper. For mutton snapper we can only use data from 2017 onward. -calibrated = read.csv("indicator_data/RVC/calibrated_species.csv") -prico = readRDS("indicator_data/RVC/prico_2001_2021_calibrated.rds") -sttstj = readRDS("indicator_data/RVC/sttstj_2001_2021_calibrated.rds") -stx = readRDS("indicator_data/RVC/stx_2001_2021_calibrated.rds") -## Make a list of species -## You can use full scientific names, common names, or -## species codes (first 3 letters of genus, and first 4 of species) -## Only scientific names are case-sensitive -spcs = c("OCY CHRY", "LUT ANAL", "BAL VETU", "EPI GUTT", "SPA AURO", "SPA VIRI") -## Code to bring in data from the server if needed for updating -## Download desired Regions data from server -#carib = getRvcData(years = 2016:2021, regions = c("PRICO","STTSTJ","STX")) -# Extract the time series for each species and each area -regions = c("PRICO","STTSTJ","STX") -for(j in regions) { -# call data for the current region -dat = switch(j, -"PRICO" = prico, -"STTSTJ" = sttstj, -"STX" = stx) -## Calculate statistics for entire sampling domain -ddens = getDomainDensity(dat, species = spcs) -# Convert SPECIES_CD to factor to ensure correct ordering in plots -ddens$SPECIES_CD <- factor(ddens$SPECIES_CD) -species <- unique(ddens$SPECIES_CD) -for (i in species) { -# Subset data for the current species -species_data <- filter(ddens, SPECIES_CD == i) -# format indicator object ----------------------------- -datdata <- species_data$YEAR -inddata <- data.frame(species_data$density) -s <- cbind(datdata, inddata) -# save ----------------------------------------- -save(s, file = paste("indicator_data/RVC/RUVdensity_", j, "_", i, ".RData", sep = "")) -} -} -# Set the directory path where the .RData files are located -folder_path <- "indicator_data/RVC" -# List all the files in the directory -files <- list.files(path = folder_path, pattern = "\\.RData$", full.names = TRUE) -# Load all .RData files -for (file in files) { -# Replace spaces with underscores in object names -object_names <- gsub(" ", "_", gsub(".RData", "", basename(file))) -# Load the file -loaded_objects <- load(file) -# Assign loaded objects to global environment with modified names -for (i in seq_along(loaded_objects)) { -assign(object_names[i], get(loaded_objects[i]), envir = .GlobalEnv) -} -} -# format indicator object ----------------------------- -# OCY CHRY --> yellowtail snapper -# LUT ANAL --> mutton snapper* -# BAL VETU --> queen triggerfish -# EPI GUTT --> red hind -# SPA AURO --> redband parrotfish -# SPA VIRI --> stoplight parrotfish -# Puerto Rico -datdata <- as.integer(RUVdensity_PRICO_BAL_VETU$datdata) -inddata <- data.frame(cbind(RUVdensity_PRICO_BAL_VETU$species_data.density, RUVdensity_PRICO_EPI_GUTT$species_data.density, RUVdensity_PRICO_LUT_ANAL$species_data.density, RUVdensity_PRICO_OCY_CHRY$species_data.density, RUVdensity_PRICO_SPA_AURO$species_data.density, RUVdensity_PRICO_SPA_VIRI$species_data.density)) -labs <- c("Puerto Rico" , "Density", "queen triggerfish", -"Puerto Rico" , "Density", "red hind", -"Puerto Rico" , "Density", "mutton snapper*", -"Puerto Rico" , "Density", "yellowtail snapper", -"Puerto Rico" , "Density", "redband parrotfish", -"Puerto Rico" , "Density", "stoplight parrotfish") -indnames <- data.frame(matrix(labs, nrow = 3, byrow = F)) -inddata <- list(labels = indnames, indicators = inddata, datelist = datdata) -class(inddata) <- "indicatordata" -save(inddata, file = "indicator_objects/RVC_PR.RData") -# St. Thomas & St. John -datdata <- as.integer(RUVdensity_PRICO_BAL_VETU$datdata) -inddata <- data.frame(cbind(RUVdensity_STTSTJ_BAL_VETU$species_data.density, RUVdensity_STTSTJ_EPI_GUTT$species_data.density, RUVdensity_STTSTJ_LUT_ANAL$species_data.density, RUVdensity_STTSTJ_OCY_CHRY$species_data.density, RUVdensity_STTSTJ_SPA_AURO$species_data.density, RUVdensity_STTSTJ_SPA_VIRI$species_data.density)) -labs <- c("St. Thomas & St. John" , "Density", "queen triggerfish", -"St. Thomas & St. John" , "Density", "red hind", -"St. Thomas & St. John" , "Density", "mutton snapper*", -"St. Thomas & St. John" , "Density", "yellowtail snapper", -"St. Thomas & St. John" , "Density", "redband parrotfish", -"St. Thomas & St. John" , "Density", "stoplight parrotfish") -indnames <- data.frame(matrix(labs, nrow = 3, byrow = F)) -inddata <- list(labels = indnames, indicators = inddata, datelist = datdata) -class(inddata) <- "indicatordata" -save(inddata, file = "indicator_objects/RVC_STSJ.RData") -# St. Croix -datdata <- as.integer(RUVdensity_PRICO_BAL_VETU$datdata) -inddata <- data.frame(cbind(RUVdensity_STX_BAL_VETU$species_data.density, RUVdensity_STX_EPI_GUTT$species_data.density, RUVdensity_STX_LUT_ANAL$species_data.density, RUVdensity_STX_OCY_CHRY$species_data.density, RUVdensity_STX_SPA_AURO$species_data.density, RUVdensity_STX_SPA_VIRI$species_data.density)) -labs <- c("St. Croix" , "Density", "queen triggerfish", -"St. Croix" , "Density", "red hind", -"St. Croix" , "Density", "mutton snapper*", -"St. Croix" , "Density", "yellowtail snapper", -"St. Croix" , "Density", "redband parrotfish", -"St. Croix" , "Density", "stoplight parrotfish") -indnames <- data.frame(matrix(labs, nrow = 3, byrow = F)) -inddata <- list(labels = indnames, indicators = inddata, datelist = datdata) -class(inddata) <- "indicatordata" -save(inddata, file = "indicator_objects/RVC_STX.RData") -### Left to troubleshoot: -# Need to take out mutton snapper and do that one separately because the data are not calibrated. Need to chop it so only use data from 2017 onward. -load("../indicator_objects/RVC_PR.RData") -plotIndicatorTimeSeries(inddata, coltoplot = 1:6, trendAnalysis = T, sublabel = T) -dat = read.csv("indicator_data/pollution_sites.csv") -head(dat) -USVI = dat[,c(1,4,7,10,13,16)] -head(USVI) -PR = dat[,c(1,3,6,9,12,15)] -head(PR) -USVI = dat[-1,c(1,4,7,10,13,16)] -head(USVI) -PR = dat[-1,c(1,3,6,9,12,15)] -head(PR) -USVI$sum = rowSums(USVI[,-1], na.rm=T) -View(PR) -str(USVI) -# make sure columns are numeric -dat[] = lapply(dat, as.numeric) -USVI = dat[-1,c(1,4,7,10,13,16)] -head(USVI) -str(USVI) -PR = dat[-1,c(1,3,6,9,12,15)] -head(PR) -USVI$sum = rowSums(USVI[,-1], na.rm=T) -head(USVI) -PR$sum = rowSums(PR[,-1], na.rm=T) -head(dat) -dat = read.csv("indicator_data/pollution_sites.csv") -head(dat) -# make sure columns are numeric -dat[] = lapply(dat, as.numeric) -USVI = dat[-1,c(1,4,7,10,13,16)] -head(USVI) -str(USVI) -PR = dat[-1,c(1,3,6,9,12,15)] -head(PR) -USVI$sum = rowSums(USVI[,-1], na.rm=T) -PR$sum = rowSums(PR[,-1], na.rm=T) -# save as indicator object ---------------------- -datdata <- (min(dat$year):max(dat$year)) -inddata <- data.frame(PR, USVI) -# save as indicator object ---------------------- -datdata <- (min(dat$year, na.rm=T):max(dat$year, na.rm=T)) -# save as indicator object ---------------------- -datdata <- dat$year -# save as indicator object ---------------------- -datdata <- dat$year[-1,] -# save as indicator object ---------------------- -datdata <- USVI$year -inddata <- data.frame(PR, USVI) -labs <- c("Pollution sites" , "Number added", "Puerto Rico", -"Pollution sites" , "Number added", "USVI") -indnames <- data.frame(matrix(labs, nrow = 3, byrow = F)) -inddata <- list(labels = indnames, indicators = inddata, datelist = datdata) -class(inddata) <- "indicatordata" -# plot and save ---------------------------------- -save(inddata, file = "indicator_objects/pollution.RData") -plotIndicatorTimeSeries(inddata, coltoplot = 1:2, sublabel = TRUE) -library(plotTimeSeries) -plotIndicatorTimeSeries(inddata, coltoplot = 1:2, sublabel = TRUE) -PR = PR$sum -USVI = USVI$sum -# save as indicator object ---------------------- -datdata <- USVI$year -inddata <- data.frame(PR, USVI) -labs <- c("Pollution sites" , "Number added", "Puerto Rico", -"Pollution sites" , "Number added", "USVI") -indnames <- data.frame(matrix(labs, nrow = 3, byrow = F)) -inddata <- list(labels = indnames, indicators = inddata, datelist = datdata) -class(inddata) <- "indicatordata" -# plot and save ---------------------------------- -save(inddata, file = "indicator_objects/pollution.RData") -plotIndicatorTimeSeries(inddata, coltoplot = 1:2, sublabel = TRUE) -# save as indicator object ---------------------- -datdata <- as.integer(USVI$year) -dat = read.csv("indicator_data/pollution_sites.csv") -head(dat) -# make sure columns are numeric -dat[] = lapply(dat, as.numeric) -USVI = dat[-1,c(1,4,7,10,13,16)] -head(USVI) -str(USVI) -PR = dat[-1,c(1,3,6,9,12,15)] -head(PR) -USVI$sum = rowSums(USVI[,-1], na.rm=T) -PR$sum = rowSums(PR[,-1], na.rm=T) -PR_ind = PR$sum -USVI_ind = USVI$sum -# save as indicator object ---------------------- -datdata <- as.integer(USVI$year) -inddata <- data.frame(PR_ind, USVI_ind) -labs <- c("Pollution sites" , "Number added", "Puerto Rico", -"Pollution sites" , "Number added", "USVI") -indnames <- data.frame(matrix(labs, nrow = 3, byrow = F)) -inddata <- list(labels = indnames, indicators = inddata, datelist = datdata) -class(inddata) <- "indicatordata" -# plot and save ---------------------------------- -save(inddata, file = "indicator_objects/pollution.RData") -plotIndicatorTimeSeries(inddata, coltoplot = 1:2, sublabel = TRUE) -load("../indicator_objects/DegreeHeatingWeeks.RData") -plotIndicatorTimeSeries(ind, coltoplot = 1:2, trendAnalysis = T, sublabel = T, dateformat = "%Y%m", plotrownum = 2, plotcolnum = 1) -library(plotTimeSeries) -plotIndicatorTimeSeries(ind, coltoplot = 1:2, trendAnalysis = T, sublabel = T, dateformat = "%Y%m", plotrownum = 2, plotcolnum = 1) -# -# Degree Heating Weeks indicators -rm(list = ls()) -dev.off() -rm(list = ls()) -load("spec_file.RData") -dev.off() -load("indicator_data/spec_file.RData") -load("~/My Github Projects/Caribbean-ESR/indicator_objects/pollution.RData") -plotIndicatorTimeSeries(inddata, coltoplot = 1:3, trendAnalysis = F, sublabel = T) -plotIndicatorTimeSeries(inddata, coltoplot = 1:2, trendAnalysis = F, sublabel = T) -devtools::install_github("https://github.com/CarissaGervasi-NOAA/plotTimeSeries.git") -devtools::install_github("https://github.com/CarissaGervasi-NOAA/plotTimeSeries.git") -devtools::install_github("https://github.com/MandyKarnauskas-NOAA/plotTimeSeries.git") -devtools::install_github("https://github.com/MandyKarnauskas-NOAA/plotTimeSeries.git") -devtools::install_github("https://github.com/MandyKarnauskas-NOAA/plotTimeSeries.git") -library(plotTimeSeries) -library(spam) -library(plotTimeSeries) -library(spam) -load("../indicator_objects/DegreeHeatingWeeks.RData") -plotIndicatorTimeSeries(ind, coltoplot = 1:2, trendAnalysis = T, sublabel = T, dateformat = "%Y%m", plotrownum = 2, plotcolnum = 1) -load("../indicator_objects/OA.RData") -plotIndicatorTimeSeries(ind, trendAnalysis = T) -load("../indicator_objects/ACEindex.RData") -plotIndicatorTimeSeries(ind, trendAnalysis = T) -load("../indicator_objects/turbidity.RData") -plotIndicatorTimeSeries(ind, coltoplot = 1:3, plotrownum = 2, trendAnalysis = T, dateformat = "%m-%Y", sublabel = T) -load("../indicator_objects/Carib_SST.RData") -plotIndicatorTimeSeries(ind, coltoplot = 1:3, trendAnalysis = T, sublabel = T, dateformat = "%m-%Y") -load("../indicator_objects/marine_debris.RData") -plotIndicatorTimeSeries(ind, trendAnalysis = F) -load("../indicator_objects/pollution.RData") -plotIndicatorTimeSeries(inddata, coltoplot = 1:2, trendAnalysis = F, sublabel = T) -load("../indicator_objects/carib_Chl.RData") -plotIndicatorTimeSeries(ind, trendAnalysis = T, dateformat = "%m-%Y") -![](/indicator_plots/Land_Use_Land_Cover_2024.jpg) -[](/indicator_plots/Land_Use_Land_Cover_2024.jpg) -load("../indicator_objects/earthquakes.RData") -plotIndicatorTimeSeries(ind, trendAnalysis = T) -load("../indicator_objects/disturbance.RData") -plotIndicatorTimeSeries(inddata, coltoplot = 1:3, trendAnalysis = T, sublabel = T) -load("../indicator_objects/disturbance.RData") -load("../indicator_objects/Sargassum.RData") -plotIndicatorTimeSeries(inddata, trendAnalysis = T) -load("../indicator_objects/sargassum_innundation_monthly_mean_hu.RData") -plotIndicatorTimeSeries(inddata, coltoplot = 1:3, sublabel = T, trendAnalysis = T, dateformat = "%Y%b") -load("../indicator_objects/hotel_occupancy_rates_USVI_and_PR.RData") -plotIndicatorTimeSeries(inddata, coltoplot = 1:2, sublabel = T, trendAnalysis = T, dateformat = "%Y%b") -load("../indicator_objects/hotel_occupancy.RData") -plotIndicatorTimeSeries(inddata, coltoplot = 1:2, sublabel = T, trendAnalysis = T, dateformat = "%Y%b") -load("../indicator_objects/hotel_occupancy_rates_USVI_and_PR.RData") -load("../indicator_objects/RVC_PR.RData") -plotIndicatorTimeSeries(inddata, coltoplot = 1:6, trendAnalysis = T, sublabel = T) -load("../indicator_objects/RVC_PR.RData") -load("../indicator_objects/RVC_PR.RData") -plotIndicatorTimeSeries(ind, coltoplot = 1:6, trendAnalysis = T, sublabel = T) -load("../indicator_objects/RVC_STSJ.RData") -plotIndicatorTimeSeries(ind, coltoplot = 1:6, trendAnalysis = T, sublabel = T) -load("../indicator_objects/RVC_STX.RData") -plotIndicatorTimeSeries(ind, coltoplot = 1:6, trendAnalysis = T, sublabel = T) -load("../indicator_objects/landings.RData") -load("../indicator_objects/fish_density.RData") -plotIndicatorTimeSeries(inddata, coltoplot = 1:3, plotrownum = 2, trendAnalysis = T, sublabel = T) -load("../indicator_objects/fish_density.RData") -plotIndicatorTimeSeries(ind, coltoplot = 1:3, plotrownum = 2, trendAnalysis = T, sublabel = T) -load("../indicator_objects/total_landings.RData") -load("../indicator_objects/OceanNAICS.RData") -load("../indicator_objects/GDP.RData") -plotIndicatorTimeSeries(inddata, coltoplot = 1:2, trendAnalysis = T, sublabel = T) -load("../indicator_objects/unemployment.RData") -plotIndicatorTimeSeries(inddata, coltoplot = 1:2, trendAnalysis = T, sublabel = T, dateformat = "%Y%b") -load("../indicator_objects/gini.RData") -load("../indicator_objects/total_rec_catch.RData") -load("../indicator_objects/NCRMP_coral_cover_richness.RData") -load("../indicator_objects/coral_spprichness_cover.RData") -load("../indicator_objects/coral_spprichness_cover.RData") -plotIndicatorTimeSeries(inddata, coltoplot = 1:4, sublabel = T, trendAnalysis = T) -plotIndicatorTimeSeries(ind, coltoplot = 1:4, sublabel = T, trendAnalysis = T) -load("../indicator_objects/Carib_SST.RData") -load("../indicator_objects/pollution.RData") -plotIndicatorTimeSeries(ind, coltoplot = 1:2, trendAnalysis = F, sublabel = T) -library(plotTimeSeries) -library(spam) -load("../indicator_objects/marine_debris.RData") -plotIndicatorTimeSeries(ind, trendAnalysis = F) -load("../indicator_objects/pollution.RData") -plotIndicatorTimeSeries(ind, coltoplot = 1:2, trendAnalysis = F, sublabel = T) -load("../indicator_objects/carib_Chl.RData") -plotIndicatorTimeSeries(ind, trendAnalysis = T, dateformat = "%m-%Y") -load("../indicator_objects/earthquakes.RData") -plotIndicatorTimeSeries(ind, trendAnalysis = T) -load("../indicator_objects/Sargassum.RData") -plotIndicatorTimeSeries(inddata, trendAnalysis = T) -plotIndicatorTimeSeries(ind, trendAnalysis = T) -load("../indicator_objects/RVC_PR.RData") -plotIndicatorTimeSeries(ind, coltoplot = 1:6, trendAnalysis = T, sublabel = T) -load("../indicator_objects/RVC_STSJ.RData") -plotIndicatorTimeSeries(ind, coltoplot = 1:6, trendAnalysis = T, sublabel = T) -load("../indicator_objects/RVC_STX.RData") -plotIndicatorTimeSeries(ind, coltoplot = 1:6, trendAnalysis = T, sublabel = T) -load("../indicator_objects/fish_density.RData") -plotIndicatorTimeSeries(ind, coltoplot = 1:3, plotrownum = 2, trendAnalysis = T, sublabel = T) -load("../indicator_objects/GDP.RData") -plotIndicatorTimeSeries(inddata, coltoplot = 1:2, trendAnalysis = T, sublabel = T) -load("../indicator_objects/GDP.RData") -plotIndicatorTimeSeries(ind, coltoplot = 1:2, trendAnalysis = T, sublabel = T) -load("../indicator_objects/unemployment.RData") -plotIndicatorTimeSeries(inddata, coltoplot = 1:2, trendAnalysis = T, sublabel = T, dateformat = "%Y%b") -load("../indicator_objects/unemployment.RData") -plotIndicatorTimeSeries(ind, coltoplot = 1:2, trendAnalysis = T, sublabel = T, dateformat = "%Y%b") -load("../indicator_objects/coral_spprichness_cover.RData") -plotIndicatorTimeSeries(ind, coltoplot = 1:4, sublabel = T, trendAnalysis = T) diff --git a/.gitignore b/.gitignore index 7fe81de..4771689 100644 --- a/.gitignore +++ b/.gitignore @@ -2,3 +2,5 @@ .Rproj.user /.quarto/ + +.Rhistory \ No newline at end of file diff --git a/indicator_objects/Carib_SST.RData b/indicator_objects/Carib_SST.RData new file mode 100644 index 0000000000000000000000000000000000000000..eaff9e6e72cae2a35cfdbe235e1d514c81bdc0e9 GIT binary patch literal 8006 zcmd^@S5y;P+n}X?qJSNhl4AiW3Mf*fB%lWoDUpt}1Q0>05F*8d2pmOf1VlOkfdhmN z(n~-(geE0O?dq^>4GjE`8;D8Tk6i+r^`S*T-vp zL9fs)#$1~wOod_>vxjB_$pJ{h?fmsv5-3+v^A4m-ta}4w--2M+0_>YAV;MB1O{meo zjvJ|8^Pn$rO$}>vXl+S)3eb&Vu3&@SusEbjLw5(qX_H5{paYno`mP#6tJ($irsd}B ztnd4tE%pmS&+Hn&nIIHLLH6w6+DQzb#nZW@J2k;4*7nP~F%2ds@iJlT*ym)PW z3i|VBU$`(*2yegtZ5AB$u4zF{WNKNkC(OMYq;ag-Q<9f?`*OzVmm90|1IqALgE`1T znjx(9DkOI)5uwx;G;^CdeY~z+-SR-Wx>fBNiB25d&5=0`815^q_+gC^a%&0{Z?=52GA^r@UbR*y7*y zsJuq04rNT^Fu$%!qYDkb!_u4+q({;3LMv^O(#82Ed_d!?WcFWdX$6=!lcef+`xOXud&9hgwjHs=;|ZDi>0&m3Ef-lmwTFI-9+ zOIRQF7jE@B+w59pDYHK0iJ9XtXDAd@}F-DLBB{*I&*ex zv|AF=%swv+mz_v*IL6GUEW=i4=+gG{g ztpF39arZ19DW;D&8tjG!!n`IwYg{gdnC$YF74l+rYP~xD;!GUQnXVY3hsn?fOS6E1 zTFz{V>*JpcufkR{IEh26@~FTcVq;(PqCib)r?EjM+0#6=p|rKV z!;&{cC#qKGs;5&sLKssNe^$MpJr8^~&VKb|G6OcCXV)NHOb9NGXC1jfkw}(jnJDK( zpZlEmyevQV;ZL_aX0p%7{nM7N*5}7`(c!oQXKAa(%8S9?Y)ysk*p|_JN7+yP69P=e zl<%h^46QVH@5U#gsI>zYUP3`g_`!pqPlKc3=FDP|V(QD8$}oTV5}vJYFZdKVm>0Vh zW9a>Jh<3NqVsk=C>yDz)&gr{atKR2l=;k!-WMtggRa7~2q3U~tvaJa^r!w2y#uT`a8rn%M~uWt)cFV%J* zo#RIX0w1kghEMYLU)oYjt88#A3Wl~oa>!)>VJ-%=e-}25`@wjg>)Q+ixIxtRoNb)0 zT6z07dG5aApHPu#{HTyorvG6t*@_0wV6UFX6o1J|RoaQmnS1_>)X>gb4B?)@$$s!W z@)5Cm>GrH<-?Mdlw75Gt8`7E<|U6yh5$+j-EoF~K#Qpd zvNo~ItayvRy!OR8tN&*7M$iT~5hUwi(~-g~$MEzXhj26JoN^Bay8 zz*s+uYu;1CWelj&m6D10oE#nh(E5c#Z@K(3PFVMAs|zUw)O^2C@}7*<`A7jE=(fBS zIMueoC6jY}F+WBVw!)RTK2{$fljOG+nVx0T#>k*xk4^Z<9@M zg0+SEmbDb=FUUCcWQlcaB$jdu(%J zT>@I%^QfVHDfW#s(X-H*iV%uM&0Wc%WVFo!z4`Ig1yUC#zd_Kvgb+<*<@hz?}HChqz(4 zOn!PR_M#flvFoZrzL{9OKSTdhI2}b;B?lQQf;`+52@ z>??j@cd{ZAi@#LF5~mhk?F_3laky9B+RxC9rK%4h=FRr|{j~7VD{sl8OXQ5|dHRmY zSLd?|JCl(`jz3y$>3q7_sK@?dBv2+e)3-piz*OKZ&EEQ7 zW8W9{*7ykTl)mmI&1bT?@Y820`5{rQ>A1TG9yZOTB4BaXxWADPzRl!vj_$rw(iNYY zW|@}1Bv?W8X1QepJ;OX#J_Gx<6|1tL zwRq9LO*Hu4<`C|i5Ir)jIsnR@tnRgMvV zbbH*faIkQ)6J7+bS;UT$PeIfh9lo}Zy2{+-q9ibKd75E2iy8O&Sw6{o&~GCa#yV)ZG~vtu=bi^e$`` zv-7b;Lt2M^vKiaia=W$nyV#6qoMS;+I3mb=4yXTRm3$Xeo<@Lno+3 z!)mP5eIx4HVr7Y*8an7P!OExZQZgadsE3>yi&|WGdd<8V7mcfePiPIqnX@8a>>0kv zbdQ;yMF_qr3|1Ee_MdcXPfuA@J3S0ktK6x-H}lFqBv-2~NH*KK3QKS9W<1R=Ofzx9@R(5?O=r@g)eU$@Af0P5~U{Q7fL` zUNo+4Mv1{wGaCdf4&oP%+4p&>^tZ>NgS+n6>f8=2lJC*KW*XWB;DIf6L=+<7qV52U z>T>x$rVV+j1w=I*m_9fu`(=hBQKP}trctioR+5My?x)~()2+kkW2{&KxKiiE#6V=N z?nIFOA!|xl+jA7eRkh8;e4n1+Gd+I0m9!$(fRlks9$q8 zZ}Nw4uYZZ9_ZBI&Sf6~exYEWY=*GVN<>)#&O+4$z6ms~0J_XGicmfx)IrWk_SGp*m z#HS&K8rTCrVfl4v&0;+2S7WD}g^J^Z9N&aD=N+y~_OCr1Mvi^&R>>dYg;(OcuK*pB zo@^ZzWS52~pI365A8yi{X$-cfK5E!y9~}5pBtSFkXDhXUJBrP~lw|6Xn}TI}TJ28V znyqU}U1^#^4R?d>q?I(<#<0ckTLd~i`7nEbRlaO! z?)3??#45KPUNQQ|;ryXfo+GsCfXVny`BtyGK=&U6SgXAmiX8gOFiQ@OKpaikW()~9 zsd;Q<(K&WvjaY?Ezu$hAZ>`DD<0A$ol5YQwwzv!`8av`Oiv?lcOKi}@7VD>g<6Oj< zl08YEAk8RJM#^Lb_vA^_KPiCuuvzwSF1-lNt^NXaXlL9YkGzqO+lW`Rre|l@&g+k{ zTp6BW=kr7BWnb-wd0*XXLRM^b!NroQPjX@K!$_a1PD}c{JBLi3bB*#)U&sy-h3i6Q zu&*&zfB4y1c88m9BA(J8DR!3>&sm1R*SFKJB$>8RG*Ow)aii>?!mwekHR8iS`M)*C z5uz&oBL#%E&_0Z?&Ey%g8h#-~ie91q-7BmN>$K1DBlDBj9NwS*y<^#5Yk zz{h4$nZ-!g3w03hl0eF|Ctc*{ewgM$?XX~*Q(I=KNy~ih;oP(ezi6SK4ynJAfZy&{ z(=?i{3%qW{lR1O0cz(E!3t!wK5r&a%%J!DZcSh5UU03&;mq#d@F<wT{-!I6p~TKKZd^zP{ORd7qeaS(Iz~<=1Zh+EOsom7k+UU2)I-aS z7VMsX?<4+=CQO{=besCZP~&UXH}%9;%=$x&6N)e{YGv~}x2$21eGs>pQ>lI!@N@e~ zr1(41V$0o>?i7`K=(rqSL|kA;#RH|+nVY`aYvjO%w~FI}2vjjF1-~p$RRkqz4)Qq9 zk%p4#irsG(5c3KBDTMwxAw|+H{{>a)*UcSU1H@%6dcR0ODNS$TfyF6=5y3O$&|YY! zJDj;ue@{y4_4b)Mh=rJM`w#sTdf|nGneWtruMR0e2vrays0dHJi3I-fIzzt1?aU}r zm=sRE2jE_Mak1logV@G6T+5Mb0#e=_TvkSZ$~CY9HlcHpmcTyk@G2&>Y2&D#qxsyj zB&<}ned(9of*NzqLs8>3JEUOx^vwL9&6n)Gy%wl`yJ1@r^R0Dbr}5n#q=9|`#dY^Z z`r!IyNYqk|-0O^v9C?}78A7D>c4~T^M?0eOW`JAi+h&N5F+Ty-tK2gg@`#Ciqe96>Ew{D zo0(u;wISQH$@vlM#Zjz&yd(D$G?RY-TZ6bjg`I^KC?E1F!E;%p}hx{3L& znjn90j*Bn~GdP;-I7rnU&L;znW|!b*LPjhHyEs&@BO*9=`fJc*Oyvd)tBBX90IGQC zFhFTiRq1*FLb|l$$|-`5vHe+>9(u2^0vo}qJh;WMylk(_EZ;4u+|WH$(c&t;;+`yY z4wn=FxP@_LRoII{v@|Wcq?SL5Rtygru3#eque}0xEJmqz!Q)Ff8*bX^!vA8d9^BGQ zqMX?Ju=hYWQgPM!saAe2yj}I2O9s5PsagRf3V7IJB6dbaVvUa+GBKbL&G;Dj;ySwr z^(4H5_RH{eIzW$}*U+^U!Ios)ELCSzKAg8)y=%W``@^v10q?R_;KP=%mOp)zRUFbk zHS!b2uVbb2;n-db_=<~Pm>Pug!qFDJKF-Xk!BMI?Q<^$IJGyDyrYymL-?p%j9MZ}ix zJyil@nmpg3M+&+znjYx|FTTd*lg5&2r-ydYJVQbLMiI70Rr3P;EGw z_>S_Cf}x~OBv3vQKh{P0R27gq>PCWO@cudseWVi2@IHX({|tRzsMzup+)5Sv9&Sj> zE5OwCZw-?Ho>-g1u$bTMlyGCj6WNrOG{PRz`sP~s?AX>9cX$SjS4=avt8?t#;~KD%7LrNcU%xxTm9 zA-#8gsWt`~_MKK6C;$8B1hgj5IkIWvih9zN7c8<)T9GYL8Q_~Qn(Vn|?V`iFPP4m! zi3Ww)lO0x~zB8nbj@jdPnAtJG28_JQUts|&P3UXKo(tsKxI=zw8yX0`QFdW1M@M(* z8jWW~%RYBo@BCP4$cCmMG=1avGdejucu#^N0@CDLUdD6L;?d_Z zVsj?Vj&qX-wo8#Qx}eSOmQRTGuk<_L*sVud!SCxHYJ7_Wt)PvL|FA3&|0pzg{Elfr zw5d4Wd>dH>*y@ZpKiH3^E$XuO64=i`5sZx{p=_37DQ+2*>>UeoI!ZqrMSKEB?rV>^ zm}D4Z4306OKDwHNO$xy)FlB}e#c)w0v%r-**!VwJt-YhcX-#EwHJ7G{V!955LnmQa?b%;1E%EFNu&);rEy^z{Q%)&w+UPku8SSi=q)Q?cLGD! z9vy1Gr>%8}Q6mv}G-K!B+yI)XiiHHB9c9HXaCAx%V%a0d1k)H^7qwq9Lek`Ce@m7 zQvLR?bl4HwrRP2;-52_OWGLkeZUimeY{q1<&NAOGdyj!Wk#d5qEkUd1moiO))k&V^ z#~G~aj0@!VlnBaSYYA}NpIiDQaxbP(=qpQXKas!NDA)w`k(Gnl=jxE(bD+rbTt-n` ztN}|(J78B-_}1>qDcbpOMji5fU54|_OT(Vl@k7OoGyhr)oCu?v1pGT{NHPXH{N zP*Gly`pKSw$KuMAO5>X%jL9c0YJJjf)a->va<{1Z!U;r26FUQJlpg$Up`7G+XKq8X znQ#(u$;8g~-%zp8OLBZ5w;{z$C4uNv1*MEVkyMJ; zi+N#k<9`k1l;ZVcR(?}6KmyUp1ZjXp>(TBO0wu@0bD7CzfF$Cw3DWjAOfrggto0tt_PR=3 zDY==SxbiJN{XgJ;o4(e1MzX!H{%!s@{J$+7lAUPWxG(!H2)^QYLq#NTvZw!XoX3jO z4F*x_k&FmGaK-l??!rH7_CL#Li(N3`N%~#JFg>OJebQH|GYixIz3~5P%k{ZRTrT0x zNrZko7i|$b@>Ie4nqU~yg7+& zuKx^QTVF}H;LZOJ`VaU2oUq802ffFSs9@jvL70I4?~aGG=$*5P!=OxH@3^AOwFd9) h*1y;k#a_hV#Vern)BHyV&sZ literal 0 HcmV?d00001 diff --git a/indicator_processing/automated_download/sst.R b/indicator_processing/automated_download/sst.R index bcffec8..bbb7794 100644 --- a/indicator_processing/automated_download/sst.R +++ b/indicator_processing/automated_download/sst.R @@ -3,7 +3,7 @@ # https://stackoverflow.com/questions/68666187/internetopenurl-failed-a-connection-with-the-server-could-not-be-established rm(list = ls()) -dev.off() +# dev.off() library(lubridate) library(maps) @@ -20,7 +20,6 @@ styear <- 1982 enyear <- terminal_year # get ERDDAP info -------------------------------- -sst <- info('ncdcOisst21Agg') sst <- info('ncdcOisst21Agg_LonPM180') # this may work better # empty data ------------------------------------------------- @@ -35,11 +34,6 @@ n <- 0 for (yr in styear:enyear) { ### BDT rERDDAP fix - sst_grab <- griddap(sst, fields = 'sst', - time = c(paste0(yr,'-01-01'), paste0(yr,'-12-31')), - longitude = c(360 + min_lon, 360 + max_lon), - latitude = c(min_lat, max_lat)) - sst_grab <- griddap(sst, fields = 'sst', time = c(paste0(yr,'-01-01'), paste0(yr,'-12-31')), longitude = c(min_lon, max_lon), @@ -53,7 +47,6 @@ for (yr in styear:enyear) { dat[m:n,] <- data.frame(sst_agg[,1:2], unlist(sst_agg[,3])) m <- n + 1 -# # # url from ERDDAP for OISST, download and read ---------------- # url <- paste0("https://coastwatch.pfeg.noaa.gov/erddap/griddap/ncdcOisst21Agg.csv?sst[(", yr, "-01-01T12:00:00Z):1:(", yr, "-12-31T12:00:00Z)][(0.0):1:(0.0)][(17):1:(19)][(292):1:(296)]") # @@ -81,7 +74,7 @@ for (yr in styear:enyear) { # file.remove("st.csv") # add row names and yearmonth column -------------------------- -names(dat) <- c("year", "mon", "mean", "min", "max") +# names(dat) <- c("year", "mon", "mean", "min", "max") dat$yrmon <- paste0(dat$mon, "-", dat$year) dat head(dat) From 1f2b8192bdf60617248d0693d97bba9176dd86d6 Mon Sep 17 00:00:00 2001 From: Brendan Turley Date: Fri, 7 Jun 2024 10:50:14 -0400 Subject: [PATCH 3/7] uncomment dev.off() which caused issues initially rerunning code --- indicator_processing/automated_download/sst.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/indicator_processing/automated_download/sst.R b/indicator_processing/automated_download/sst.R index bbb7794..153a01d 100644 --- a/indicator_processing/automated_download/sst.R +++ b/indicator_processing/automated_download/sst.R @@ -3,7 +3,7 @@ # https://stackoverflow.com/questions/68666187/internetopenurl-failed-a-connection-with-the-server-could-not-be-established rm(list = ls()) -# dev.off() +dev.off() library(lubridate) library(maps) From f416c4fcf90ce687d96a73e7c840380d3813971d Mon Sep 17 00:00:00 2001 From: Brendan Turley Date: Fri, 7 Jun 2024 12:14:27 -0400 Subject: [PATCH 4/7] update kd490 code to use rerddap package to fix timeout issues --- .../automated_download/kd490.R | 204 +++++++++++++----- 1 file changed, 145 insertions(+), 59 deletions(-) diff --git a/indicator_processing/automated_download/kd490.R b/indicator_processing/automated_download/kd490.R index 663d169..836e5c7 100644 --- a/indicator_processing/automated_download/kd490.R +++ b/indicator_processing/automated_download/kd490.R @@ -7,74 +7,81 @@ # specification file and libraries ----------------------------- rm(list = ls()) -dev.off() +# dev.off() +library(lubridate) library(maps) library(plotTimeSeries) +library(rerddap) load("indicator_processing/spec_file.RData") -# define urls -------------------------------- - -beg <- "https://coastwatch.pfeg.noaa.gov/erddap/griddap/nesdisVHNSQkd490Monthly.csv?kd_490[(2012-01-02T12:00:00Z):1:(" - -url_pr <- paste0(beg, terminal_year, "-12-01T12:00:00Z)][(0.0):1:(0.0)][(18.7):1:(17.7)][(-67.5):1:(-65.2)]") -url_st <- paste0(beg, terminal_year, "-12-01T12:00:00Z)][(0.0):1:(0.0)][(18.41):1:(18.20)][(-65.12):1:(-64.66)]") -url_sc <- paste0(beg, terminal_year, "-12-01T12:00:00Z)][(0.0):1:(0.0)][(17.86):1:(17.59)][(-64.94):1:(-64.41)]") - -# previously separated St. Thomas and St. John but highly correlated -#url_sj <- paste0(beg, terminal_year, "-12-01T12:00:00Z)][(0.0):1:(0.0)][(18.37):1:(18.20)][(-64.8):1:(-64.66)]") - -# automated downloads ------------------ - -download.file(url = url_pr, destfile = "test_pr.csv") -download.file(url = url_st, destfile = "test_st.csv") -#download.file(url = url_sj, destfile = "test_sj.csv") -download.file(url = url_sc, destfile = "test_sc.csv") - -# read files ---------------------------- -kdlabs <- read.table("test_pr.csv", sep = ",", header = T, skip = 0) -kd_pr <- read.table("test_pr.csv", sep = ",", header = T, skip = 1) -kd_st <- read.table("test_st.csv", sep = ",", header = T, skip = 1) -#kd_sj <- read.table("test_sj.csv", sep = ",", header = T, skip = 1) -kd_sc <- read.table("test_sc.csv", sep = ",", header = T, skip = 1) - -# remove .csv files --------------------- -file.remove("test_pr.csv") -file.remove("test_st.csv") -#file.remove("test_sj.csv") -file.remove("test_sc.csv") - -# assign labels to data frames ---------- -names(kd_pr) <- names(kdlabs) -names(kd_st) <- names(kdlabs) -#names(kd_sj) <- names(kdlabs) -names(kd_sc) <- names(kdlabs) - -head(kd_pr) -head(kd_st) -#head(kd_sj) -head(kd_sc) - -# calculate means ----------------------- -ind_pr <- tapply(kd_pr$kd_490, kd_pr$time, mean, na.rm = T) -ind_st <- tapply(kd_st$kd_490, kd_st$time, mean, na.rm = T) -#ind_sj <- tapply(kd_sj$kd_490, kd_sj$time, mean, na.rm = T) -ind_sc <- tapply(kd_sc$kd_490, kd_sc$time, mean, na.rm = T) - -# check names are the same -------------- -summary(names(ind_pr) == names(ind_st)) -#summary(names(ind_pr) == names(ind_sj)) -summary(names(ind_pr) == names(ind_sc)) - -#uli <- tapply(kd$kd_490, kd$time, function(x) quantile(x, probs = c(0.95), na.rm = T)) -#lli <- tapply(kd$kd_490, kd$time, function(x) quantile(x, probs = c(0.05), na.rm = T)) +# define years -------------------------------- +styear <- 2012 +enyear <- terminal_year + +# define geographic scope -------------------------------- +pr_coord <- matrix(c(-67.5, -65.2, 17.7, 18.7), 2, 2) +sc_coord <- matrix(c(-64.94, -64.41, 17.59, 17.86), 2, 2) +st_coord <- matrix(c(-65.12, -64.66, 18.20, 18.41), 2, 2) + +# get ERDDAP info -------------------------------- +kd490 <- info('nesdisVHNSQkd490Monthly', + url = 'https://coastwatch.pfeg.noaa.gov/erddap') +# alternative with longer timeseries +# kd490 <- info('pmlEsaCCI60OceanColorDaily', +# url = 'https://coastwatch.pfeg.noaa.gov/erddap') + + +# empty data ------------------------------------------------- +pr_dat <- sc_dat <- st_dat <- setNames(data.frame(matrix(NA,length(styear:enyear)*12,3)), + c("year", "mon", 'mean')) +m <- 1 +n <- 0 + +# download by year to avoid timeout errors -------------------- +for (yr in styear:enyear) { + + ### BDT rERDDAP fix + n <- n + 12 + ### Puerto Rico + pr_kd490_grab <- griddap(kd490, fields = 'kd_490', + time = c(paste0(yr,'-01-15'), paste0(yr,'-12-15')), + longitude = pr_coord[ ,1], + latitude = pr_coord[ ,2]) + + pr_dat[m:n,] <- aggregate(pr_kd490_grab$data$kd_490, + by = list(year(pr_kd490_grab$data$time), month(pr_kd490_grab$data$time)), + mean, na.rm = T) + + ### St Croix + sc_kd490_grab <- griddap(kd490, fields = 'kd_490', + time = c(paste0(yr,'-01-15'), paste0(yr,'-12-15')), + longitude = sc_coord[ ,1], + latitude = sc_coord[ ,2]) + + sc_dat[m:n,] <- aggregate(sc_kd490_grab$data$kd_490, + by = list(year(sc_kd490_grab$data$time), month(sc_kd490_grab$data$time)), + mean, na.rm = T) + + ### St Thomas + st_kd490_grab <- griddap(kd490, fields = 'kd_490', + time = c(paste0(yr,'-01-15'), paste0(yr,'-12-15')), + longitude = st_coord[ ,1], + latitude = st_coord[ ,2]) + + st_dat[m:n,] <- aggregate(st_kd490_grab$data$kd_490, + by = list(year(st_kd490_grab$data$time), month(st_kd490_grab$data$time)), + mean, na.rm = T) + + m <- n + 1 +} # create indicator object -------------- -datdata <- strftime(names(ind_pr), format = "%m-%Y") -inddata <- data.frame(cbind(ind_pr, ind_st, ind_sc), row.names = datdata) -labs <- c(rep("Turbidity from ocean color data", 3), rep("Diffuse attenuation coefficient at 490 nm (kd^2)", 3), +datdata <- paste(sprintf('%02d',pr_dat$mon),pr_dat$year,sep='-') +inddata <- data.frame(cbind(pr_dat$mean, st_dat$mean, sc_dat$mean), row.names = datdata) +labs <- c(rep("Turbidity from ocean color data", 3), rep("Diffuse attenuation coefficient at 490 nm (m^-1)", 3), "Puerto Rico", "St. Thomas", "St. Croix") indnames <- data.frame(matrix(labs, nrow = 3, byrow = T)) @@ -90,3 +97,82 @@ plotIndicatorTimeSeries(ind, coltoplot = 1:3, plotrownum = 3, dateformat = "%m-% + + +# terminal_year <- 2013 +# # define urls -------------------------------- +# +# beg <- "https://coastwatch.pfeg.noaa.gov/erddap/griddap/nesdisVHNSQkd490Monthly.csv?kd_490[(2012-01-02T12:00:00Z):1:(" +# +# url_pr <- paste0(beg, terminal_year, "-12-01T12:00:00Z)][(0.0):1:(0.0)][(18.7):1:(17.7)][(-67.5):1:(-65.2)]") +# url_st <- paste0(beg, terminal_year, "-12-01T12:00:00Z)][(0.0):1:(0.0)][(18.41):1:(18.20)][(-65.12):1:(-64.66)]") +# url_sc <- paste0(beg, terminal_year, "-12-01T12:00:00Z)][(0.0):1:(0.0)][(17.86):1:(17.59)][(-64.94):1:(-64.41)]") +# +# # previously separated St. Thomas and St. John but highly correlated +# #url_sj <- paste0(beg, terminal_year, "-12-01T12:00:00Z)][(0.0):1:(0.0)][(18.37):1:(18.20)][(-64.8):1:(-64.66)]") +# +# # automated downloads ------------------ +# +# download.file(url = url_pr, destfile = "test_pr.csv") +# download.file(url = url_st, destfile = "test_st.csv") +# #download.file(url = url_sj, destfile = "test_sj.csv") +# download.file(url = url_sc, destfile = "test_sc.csv") +# +# # read files ---------------------------- +# kdlabs <- read.table("test_pr.csv", sep = ",", header = T, skip = 0) +# kd_pr <- read.table("test_pr.csv", sep = ",", header = T, skip = 1) +# kd_st <- read.table("test_st.csv", sep = ",", header = T, skip = 1) +# #kd_sj <- read.table("test_sj.csv", sep = ",", header = T, skip = 1) +# kd_sc <- read.table("test_sc.csv", sep = ",", header = T, skip = 1) +# +# # remove .csv files --------------------- +# file.remove("test_pr.csv") +# file.remove("test_st.csv") +# #file.remove("test_sj.csv") +# file.remove("test_sc.csv") +# +# # assign labels to data frames ---------- +# names(kd_pr) <- names(kdlabs) +# names(kd_st) <- names(kdlabs) +# #names(kd_sj) <- names(kdlabs) +# names(kd_sc) <- names(kdlabs) +# +# head(kd_pr) +# head(kd_st) +# #head(kd_sj) +# head(kd_sc) +# +# # calculate means ----------------------- +# ind_pr <- tapply(kd_pr$kd_490, kd_pr$time, mean, na.rm = T) +# ind_st <- tapply(kd_st$kd_490, kd_st$time, mean, na.rm = T) +# #ind_sj <- tapply(kd_sj$kd_490, kd_sj$time, mean, na.rm = T) +# ind_sc <- tapply(kd_sc$kd_490, kd_sc$time, mean, na.rm = T) +# +# # check names are the same -------------- +# summary(names(ind_pr) == names(ind_st)) +# #summary(names(ind_pr) == names(ind_sj)) +# summary(names(ind_pr) == names(ind_sc)) +# +# #uli <- tapply(kd$kd_490, kd$time, function(x) quantile(x, probs = c(0.95), na.rm = T)) +# #lli <- tapply(kd$kd_490, kd$time, function(x) quantile(x, probs = c(0.05), na.rm = T)) +# +# # create indicator object -------------- +# +# datdata <- strftime(names(ind_pr), format = "%m-%Y") +# inddata <- data.frame(cbind(ind_pr, ind_st, ind_sc), row.names = datdata) +# labs <- c(rep("Turbidity from ocean color data", 3), rep("Diffuse attenuation coefficient at 490 nm (kd^2)", 3), +# "Puerto Rico", "St. Thomas", "St. Croix") +# indnames <- data.frame(matrix(labs, nrow = 3, byrow = T)) +# +# ind <- list(labels = indnames, indicators = inddata, datelist = datdata) #, ulim = ulidata, llim = llidata) +# class(ind) <- "indicatordata" +# +# # save and plot --------------------------------------- +# +# save(ind, file = "indicator_objects/turbidity.RData") +# +# plotIndicatorTimeSeries(ind, coltoplot = 1:3, plotrownum = 3, dateformat = "%m-%Y", sublabel = T, trendAnalysis = F, widadj = 1.5) +# plotIndicatorTimeSeries(ind, coltoplot = 1:3, plotrownum = 3, dateformat = "%m-%Y", sublabel = T, trendAnalysis = F, widadj = 1.5, anom = "mon") +# + + From ebb075c46fcbad6f7f7ac6aa28821754a33fc1f9 Mon Sep 17 00:00:00 2001 From: Brendan Turley Date: Fri, 7 Jun 2024 12:25:10 -0400 Subject: [PATCH 5/7] create turbindity indicator object; uncomment dev.off() casuing issues --- indicator_objects/turbidity.RData | Bin 0 -> 3976 bytes indicator_processing/automated_download/kd490.R | 2 +- 2 files changed, 1 insertion(+), 1 deletion(-) create mode 100644 indicator_objects/turbidity.RData diff --git a/indicator_objects/turbidity.RData b/indicator_objects/turbidity.RData new file mode 100644 index 0000000000000000000000000000000000000000..4132814583ead7d9bfd1a659326e82f8afc26f45 GIT binary patch literal 3976 zcmbu?`#%$ozX$N4wsK7oxovz>C>CzUQ1D&iSdI&ijYwDU-q=+rEiXC#2>zF=q|ua=KnrSuO|elNtA895<5| z`Uj4QyOi4u7loEZ7eAj>ON@?jE{`#0I5|akJ2RYgxs6|!e0Ma%NJyy1Y5ETG&gI9& zWAB}_QNg-04%<$i)@zx`a^tqjFbStVLq{UAEx004?#Jk{)A_h#Sq@@K(uLkgC4q<@ zMt5{)pSz56L`2`oQ3{Rmim4Pm<{-oHH!uD)2h<;HMn7zKxUXBVnKKQPQ+QhHfC>wP zjl~!Fx^Pg>bW zwY_*0U`ok(bkyW6BK&Tru57TPgrOx^z0?|@X96|q@J}jvA3pkhvb7v8sa#mEqyo;o zTOX{8|2N9`UAYYHgC>F)bf}twW}aM%I^-iSFemk3$u{zA5^q|8&jqhXJI#!IS=$6c zTDKo7D1S#R@p}I%3%=Pem;bicT|M??2GQ;JzBlHg;eB2udORWu^c^^Y`3tdU6=f7g!ije4t+w_BZZ!xEfzKofZo3#N911c7GAP zc_N%t`6#x??j!3`gtyW{AObM;bhn-mA)n^_nR3S<1yEoAWFdPD8s+jB4>bleSBU+h z)MZA%g`3nedgaf&p4v#4X~=>!9=iV|=bx{hjWe4M-j050U!0y;t&3>NIL_LKYNzWw zo9YnYP5Zh~dF@@@8EwaLMU~de37QXLR$OzKJ+{=jmc-$MjEwz9`Y}>et$rFc;70<# z|6UV;Wy+#aM|Fe)k@w8@EAIm9D@#3(3>AGXSNYmbpXZ>eM}PjH$LxXCcWVg((Vt2h z8dXAR43&~Oxi~8oftWwt3u_s|`*4qR12d{@xSXVYwtMFbAA8zK6@QLG1eyrI%5~ZOx7#}vwxO6n^7_jL z6Anm(LG13lGMmd*0A>MD{p@;Sdl3h;=?bELXpXv5#G%*g77aF)inDczXDsFPfd=QH*Qm2KVJS6Vnr~3SV01?@l;EZbz0oWN@KBZF~;whNCDvE; zWvEe7EZUAqNCS__U!qk_-2LYnY^#?bW+F)0_a0m!YuE`=-!St|VZ4#HTqVXOB%1un zM+`#prYpNGBhxgVDcm*;0D9!U!5d!#%@!@^;ZlMjSvcIZ>j%tDm-kP;Y2T1YGT_Lq z%PvxqZceo@4`P1Ms!!SpVt0gb3G+UxXTHVgYsZrMFb8>-1}7tTO4?w;?BUM(?{viV z^~AU(Ec2FUG-o)>!8*Nd0?raix3?*GLZjt;JoD#GP;@Yn-d7za3LzFgyK8ix%lx15 zUY(UI-|tAznNYR6Qi<5^gL9{!B%L$bm`ulI!8KvP1Ree4_TyBRLsLS^py;2r$quu$ zWp7|_rP~_|a?Yu_DNcm?=QGO5kwxNq$MlnFUi#|V`69(0Cbl6$kyrS0`}l%YPiE?_ zPyC&*Fi;xk8mHkVuK3uhbu4isxkOn5I*3|Fr8x|=xCo7}CtkaycxkEnc!uh{3n{i~ z{*kq7iCBy~yKh+x8ZR@s9ihc*kUeazNu)nS*>RjF@6GVlWyfU5*tefK_!ze_75R~6 zOJl1y!;>j3dCCi|slsX)j71jjK8d7=iU~y4GNbnE+h_fg&x@@jx=5=xp@+p?2%KQOG@{j1}8-Ka=>QwF}mK=FT?TO zS!!Zmwka_OQJzu?Kl!?Zxo|v{oh-v;RasQN#5zCKiBs2m*(|w!nV)=}8*NvSr0U4) z6I@BK@h^`yP-o?aW`0%o5~I4dr`~nk61d8B`$mVuEav@$X!Oo=F868vcQvInhxSdf zMR2q^6(wa!;^@DQ%<54KkD5X3rD4ID1B`q2LvTbWl2quLt}tL)V{3eo3n{XIUZ)z`@%{=E61Lsd@22YT@0TVZ2o&SeX*yQE>WoQc zA36v~2VI!0clGr?jV^nUKmX3Y16WacTIV@V9Mo-fq>G1}4N1a=^O`?-RG1l>$ z$H=0yj-lt@4u{fbmbmR-X25IFQ`ikZ&Z~PFBF5-lz|w&8h#alOeR=spZ>>mLp;0?> zSdpXuhV?*}hmKYdeuJd;I13^@kXKU_2|9Y=G}kj2JA_4v=Sa7AOKk!6sK}p=G=2C{ z<_JHd8_c?cd!@8`!WZwI=zPQ^@9*p2t4GAgKIs`jIA{;}{Y@a!{dS1#q_kIzLb^Xr##*(S+HaOBqH7^S=xlpsEsGZqTJ zn0YjYG=_0&<;|1uBFyh+a+7!00iK~jaudXg*}hYY{l>L~$VVHlX^x}Sd%-Paw=dS_ zc>j}@Mh50c-zr(Z`I>}VEOGx66!60-=WFOG#Ifdju9}L<7)2F%F#!Bi$q-x!b&vh0 zJi8%xomaX>Xo^;K9C0-OJZ>Gcn@iKJ1=~cx=1(~_j>EYfqk;+8-77{>`I?v7Fao10 zjIQa95`K-a#TUgJQKK|ii&j3fm~1$Fr+|PNVvW#~Lg8;t+{3qcG9g&Ve7YjOKtGgx zM!^q+Il{qa|MCSAnjhDRx4E}F%#EL`j6W2bFPyTqa|bO>1`vyGmCm|{oL$ZG z!AQmj_Bwoai+90Kwu2%e&8=48t;ZYK^+MJQIs4#>pgyTP)hnkpqOG}BOx6^Y-|iBn zWjB;OLoti1QsIvpuJ;#Q{5%|z{ms|ZX>f?+Qt|6!yG2F4&BL#z@AY=c#-|A#J^st~ zGR)^6RHHdr*ow>kH!yj^h(;FAAUNUw#$kOz_~c)?2`4a5 z{vBO{_T Date: Tue, 11 Jun 2024 15:46:40 -0400 Subject: [PATCH 6/7] clean up code --- indicator_processing/automated_download/sst.R | 40 ++----------------- 1 file changed, 4 insertions(+), 36 deletions(-) diff --git a/indicator_processing/automated_download/sst.R b/indicator_processing/automated_download/sst.R index 153a01d..facffa9 100644 --- a/indicator_processing/automated_download/sst.R +++ b/indicator_processing/automated_download/sst.R @@ -1,4 +1,3 @@ - # fix for dealing with internet timeout error # https://stackoverflow.com/questions/68666187/internetopenurl-failed-a-connection-with-the-server-could-not-be-established @@ -12,9 +11,6 @@ library(rerddap) load("indicator_processing/spec_file.RData") -# devtools::install_github("mdsumner/ncdf4") -# library(ncdf4) - # define years -------------------------------- styear <- 1982 enyear <- terminal_year @@ -23,8 +19,6 @@ enyear <- terminal_year sst <- info('ncdcOisst21Agg_LonPM180') # this may work better # empty data ------------------------------------------------- -# dat <- data.frame(row.names = c("year", "mon", "PR_mean", "PR_min", "PR_max", "VI_mean", "VI_min", "VI_max")) - dat <- setNames(data.frame(matrix(NA,length(styear:enyear)*12,5)), c("year", "mon", 'mean', 'min', 'max')) m <- 1 @@ -46,35 +40,10 @@ for (yr in styear:enyear) { n <- n + 12 dat[m:n,] <- data.frame(sst_agg[,1:2], unlist(sst_agg[,3])) m <- n + 1 - -# # url from ERDDAP for OISST, download and read ---------------- -# url <- paste0("https://coastwatch.pfeg.noaa.gov/erddap/griddap/ncdcOisst21Agg.csv?sst[(", yr, "-01-01T12:00:00Z):1:(", yr, "-12-31T12:00:00Z)][(0.0):1:(0.0)][(17):1:(19)][(292):1:(296)]") -# -# # url_vi <- paste0("https://coastwatch.pfeg.noaa.gov/erddap/griddap/ncdcOisst21Agg.csv?sst[(", yr, "-01-01T12:00:00Z):1:(", yr, "-12-31T12:00:00Z)][(0.0):1:(0.0)][(17):1:(19)][(294.75):1:(296)]") -# # url_pr <- paste0("https://coastwatch.pfeg.noaa.gov/erddap/griddap/ncdcOisst21Agg.csv?sst[(", yr, "-01-01T12:00:00Z):1:(", yr, "-12-31T12:00:00Z)][(0.0):1:(0.0)][(17):1:(19)][(292):1:(294.75)]") -# # download for entire Caribbean because USVI and PR highly correlated. -# download.file(url = url, destfile = "st.csv") -# -# labs <- read.table("st.csv", sep = ",", header = T, skip = 0) -# sst_vi <- read.table("st.csv", sep = ",", header = T, skip = 1) -# names(sst_vi) <- names(labs) -# -# # extract month ----------------------------------------------- -# sst_vi$mon <- strftime(sst_vi$time, format = "%m") -# -# # calculate mean, min, max and concatenate -------------------- -# ind_vi <- as.numeric(tapply(sst_vi$sst, sst_vi$mon, mean, na.rm = T)) -# min_vi <- as.numeric(tapply(sst_vi$sst, sst_vi$mon, min, na.rm = T)) -# max_vi <- as.numeric(tapply(sst_vi$sst, sst_vi$mon, max, na.rm = T)) -# -# tempdat <- data.frame(yr, unique(sst_vi$mon), ind_vi, min_vi, max_vi, stringsAsFactors = F) -# dat <- data.frame(rbind(dat, tempdat), stringsAsFactors = F) -} -# file.remove("st.csv") +} -# add row names and yearmonth column -------------------------- -# names(dat) <- c("year", "mon", "mean", "min", "max") +# add yearmonth column -------------------------- dat$yrmon <- paste0(dat$mon, "-", dat$year) dat head(dat) @@ -83,7 +52,7 @@ tail(dat) # check outputs and look at correlations --------------------- table(dat$year) table(dat$mon) -matplot(dat[3:5], type = "l") #, col = c(4, 4, 4, 2, 2, 2), lty = c(1, 2, 3, 1, 2, 3)) +matplot(dat[3:5], type = "l") cor(dat[3:5]) plot(dat[3:5]) @@ -91,13 +60,12 @@ plot(dat[3:5]) labs <- c(rep("U.S. Caribbean sea surface temperature", 3), rep("degrees Celsius", 3), "monthly mean", "monthly minimum", "monthly maximum") - # "Puerto Rico - mean", "USVI - mean", "Puerto Rico - monthly minimum", "USVI - monthly minimum", "Puerto Rico - monthly maximum", "USVI - monthly maximum") indnames <- data.frame(matrix(labs, nrow = 3, byrow = T)) inddata <- data.frame(dat[c(3:5)]) datdata <- dat$yrmon -s <- list(labels = indnames, indicators = inddata, datelist = datdata) #, ulim = ulidata, llim = llidata) +s <- list(labels = indnames, indicators = inddata, datelist = datdata) class(s) <- "indicatordata" s From cc5edd13825e01d6727ba35954f2429f35ac4554 Mon Sep 17 00:00:00 2001 From: Brendan Turley Date: Tue, 11 Jun 2024 15:47:32 -0400 Subject: [PATCH 7/7] remove depreciated code --- .../automated_download/kd490.R | 82 ------------------- 1 file changed, 82 deletions(-) diff --git a/indicator_processing/automated_download/kd490.R b/indicator_processing/automated_download/kd490.R index 1ccb269..cce46b2 100644 --- a/indicator_processing/automated_download/kd490.R +++ b/indicator_processing/automated_download/kd490.R @@ -1,4 +1,3 @@ -# # Turbidity indicator # # direct download from ERDDAP @@ -95,84 +94,3 @@ save(ind, file = "indicator_objects/turbidity.RData") plotIndicatorTimeSeries(ind, coltoplot = 1:3, plotrownum = 3, dateformat = "%m-%Y", sublabel = T, trendAnalysis = F, widadj = 1.5) plotIndicatorTimeSeries(ind, coltoplot = 1:3, plotrownum = 3, dateformat = "%m-%Y", sublabel = T, trendAnalysis = F, widadj = 1.5, anom = "mon") - - - - -# terminal_year <- 2013 -# # define urls -------------------------------- -# -# beg <- "https://coastwatch.pfeg.noaa.gov/erddap/griddap/nesdisVHNSQkd490Monthly.csv?kd_490[(2012-01-02T12:00:00Z):1:(" -# -# url_pr <- paste0(beg, terminal_year, "-12-01T12:00:00Z)][(0.0):1:(0.0)][(18.7):1:(17.7)][(-67.5):1:(-65.2)]") -# url_st <- paste0(beg, terminal_year, "-12-01T12:00:00Z)][(0.0):1:(0.0)][(18.41):1:(18.20)][(-65.12):1:(-64.66)]") -# url_sc <- paste0(beg, terminal_year, "-12-01T12:00:00Z)][(0.0):1:(0.0)][(17.86):1:(17.59)][(-64.94):1:(-64.41)]") -# -# # previously separated St. Thomas and St. John but highly correlated -# #url_sj <- paste0(beg, terminal_year, "-12-01T12:00:00Z)][(0.0):1:(0.0)][(18.37):1:(18.20)][(-64.8):1:(-64.66)]") -# -# # automated downloads ------------------ -# -# download.file(url = url_pr, destfile = "test_pr.csv") -# download.file(url = url_st, destfile = "test_st.csv") -# #download.file(url = url_sj, destfile = "test_sj.csv") -# download.file(url = url_sc, destfile = "test_sc.csv") -# -# # read files ---------------------------- -# kdlabs <- read.table("test_pr.csv", sep = ",", header = T, skip = 0) -# kd_pr <- read.table("test_pr.csv", sep = ",", header = T, skip = 1) -# kd_st <- read.table("test_st.csv", sep = ",", header = T, skip = 1) -# #kd_sj <- read.table("test_sj.csv", sep = ",", header = T, skip = 1) -# kd_sc <- read.table("test_sc.csv", sep = ",", header = T, skip = 1) -# -# # remove .csv files --------------------- -# file.remove("test_pr.csv") -# file.remove("test_st.csv") -# #file.remove("test_sj.csv") -# file.remove("test_sc.csv") -# -# # assign labels to data frames ---------- -# names(kd_pr) <- names(kdlabs) -# names(kd_st) <- names(kdlabs) -# #names(kd_sj) <- names(kdlabs) -# names(kd_sc) <- names(kdlabs) -# -# head(kd_pr) -# head(kd_st) -# #head(kd_sj) -# head(kd_sc) -# -# # calculate means ----------------------- -# ind_pr <- tapply(kd_pr$kd_490, kd_pr$time, mean, na.rm = T) -# ind_st <- tapply(kd_st$kd_490, kd_st$time, mean, na.rm = T) -# #ind_sj <- tapply(kd_sj$kd_490, kd_sj$time, mean, na.rm = T) -# ind_sc <- tapply(kd_sc$kd_490, kd_sc$time, mean, na.rm = T) -# -# # check names are the same -------------- -# summary(names(ind_pr) == names(ind_st)) -# #summary(names(ind_pr) == names(ind_sj)) -# summary(names(ind_pr) == names(ind_sc)) -# -# #uli <- tapply(kd$kd_490, kd$time, function(x) quantile(x, probs = c(0.95), na.rm = T)) -# #lli <- tapply(kd$kd_490, kd$time, function(x) quantile(x, probs = c(0.05), na.rm = T)) -# -# # create indicator object -------------- -# -# datdata <- strftime(names(ind_pr), format = "%m-%Y") -# inddata <- data.frame(cbind(ind_pr, ind_st, ind_sc), row.names = datdata) -# labs <- c(rep("Turbidity from ocean color data", 3), rep("Diffuse attenuation coefficient at 490 nm (kd^2)", 3), -# "Puerto Rico", "St. Thomas", "St. Croix") -# indnames <- data.frame(matrix(labs, nrow = 3, byrow = T)) -# -# ind <- list(labels = indnames, indicators = inddata, datelist = datdata) #, ulim = ulidata, llim = llidata) -# class(ind) <- "indicatordata" -# -# # save and plot --------------------------------------- -# -# save(ind, file = "indicator_objects/turbidity.RData") -# -# plotIndicatorTimeSeries(ind, coltoplot = 1:3, plotrownum = 3, dateformat = "%m-%Y", sublabel = T, trendAnalysis = F, widadj = 1.5) -# plotIndicatorTimeSeries(ind, coltoplot = 1:3, plotrownum = 3, dateformat = "%m-%Y", sublabel = T, trendAnalysis = F, widadj = 1.5, anom = "mon") -# - -