Merge branch 'main' of https://github.com/Gulf-IEA/Caribbean-ESR

indicator_processing/CalculateAllIndicators.R

@@ -0,0 +1,57 @@
+# 
+# script to update ALL indicators for report 
+# M. Karnauskas Mar 3 2024
+
+# load libraries ---------------------------------
+rm(list = ls())
+
+library(plotTimeSeries)
+
+#  find root directory for project ---------------
+
+directory <- rprojroot::find_rstudio_root_file()
+
+# first process automated downloads --------------
+
+setwd(directory)
+setwd("indicator_processing/automated_download/")
+dir()
+
+# run all scripts in folder ---------------------
+
+plot(1)
+
+source("ACE_index_Carib.R")       # hurricane energy index
+source("chl_caribbean.R")         # primary productivity
+source("DHW.R")                   # degree heating weeks
+source("earthquakes.R")           # earthquakes
+source("kd490.R")                 # turbidity from Kd490
+source("sst.R")                   # sea surface temperature
+
+
+###############################################################
+
+# Notes for standardizing scripts:  
+
+#  Top: 
+
+# specification file and libraries -----------------------------
+rm(list = ls())
+dev.off()
+
+library(maps)
+library(plotTimeSeries)
+
+load("spec_file.RData")
+
+# define years  --------------------------------
+styear <- 1961
+enyear <- terminal_year
+
+
+#  Bottom: 
+#  save all indicators as "ind" object
+
+save(ind, file = "../../indicator_objects/INDICATOR_NAME.RData")
+
+##

indicator_processing/automated_download/ACE_index_Carib.R

@@ -13,26 +13,29 @@
 
 # specification file and libraries -----------------------------
 rm(list = ls())
+dev.off()
+
 library(maps)
 library(plotTimeSeries)
 
-#load("spec_file.RData")
+load("spec_file.RData")
 
 # define years  --------------------------------
 styear <- 1961
-enyear <- 2023
+enyear <- terminal_year
 
 # download data directly from site -----------------------------
 options(download.file.method="libcurl")
 
-url <- "https://www.ncei.noaa.gov/data/international-best-track-archive-for-climate-stewardship-ibtracs/v04r00/access/csv/ibtracs.NA.list.v04r00.csv"
 url <- "https://www.ncei.noaa.gov/data/international-best-track-archive-for-climate-stewardship-ibtracs/v04r00/access/csv/ibtracs.NA.list.v04r00.csv"
 
-download.file(url = url, destfile = "C:/Users/mandy.karnauskas/Desktop/Caribbean-ESR/indicator_data/ibtracs.csv")
-dat <- read.csv("C:/Users/mandy.karnauskas/Desktop/Caribbean-ESR/indicator_data/ibtracs.csv", skip = 2, header = F)
-datn <- read.csv("C:/Users/mandy.karnauskas/Desktop/Caribbean-ESR/indicator_data/ibtracs.csv", skip = 0, header = T)
+download.file(url = url, destfile = "ibtracs.csv")
+dat <- read.csv("ibtracs.csv", skip = 2, header = F)
+datn <- read.csv("ibtracs.csv", skip = 0, header = T)
 names(dat) <- names(datn)
 
+file.remove("ibtracs.csv")
+
 # cut columns --------------------------------
 
 head(dat)
@@ -152,12 +155,10 @@ indnames <- data.frame(matrix(labs, nrow = 3, byrow = F))
 s <- list(labels = indnames, indicators = inddata, datelist = datdata) #, ulim = ulidata, llim = llidata)
 class(s) <- "indicatordata"
 
-setwd("C:/Users/mandy.karnauskas/Desktop/Caribbean-ESR/indicator_plots/")
-
-plotIndicatorTimeSeries(s, widadj = 0.5, outtype = "png")
+plotIndicatorTimeSeries(s, widadj = 0.5)
 
-inddata <- s
-save(inddata, file = "C:/Users/mandy.karnauskas/Desktop/Caribbean-ESR/indicator_objects/ACEindex.RData")
+ind <- s
+save(ind, file = "../../indicator_objects/ACEindex.RData")
 
 #################################################################################
 

indicator_processing/automated_download/DHW.R

@@ -1,10 +1,15 @@
-# DHW indicators
+#
+# Degree Heating Weeks indicators
 
 rm(list = ls())
+dev.off()
+
+load("spec_file.RData")
+
 options(download.file.method="libcurl")
 
 styear <- 2012
-enyear <- 2023
+enyear <- terminal_year
 
 # download data and calculate mean -------------------
 
@@ -13,10 +18,12 @@ for (i in 1:2) {
   if (i == 1) {  url <- "https://coralreefwatch.noaa.gov/product/vs/data/puerto_rico.txt"  }
   if (i == 2) {  url <- "https://coralreefwatch.noaa.gov/product/vs/data/usvi.txt"  }
 
-  download.file(url = url, destfile = "C:/Users/mandy.karnauskas/Downloads/dhw.txt")
-  d <- read.table("C:/Users/mandy.karnauskas/Downloads/dhw.txt", skip = 21, header = T)
+  download.file(url = url, destfile = "dhw.txt")
+  d <- read.table("dhw.txt", skip = 21, header = T)
+
+  file.remove("dhw.txt")
 
-  d <- d[which(d$YYYY >= styear & d$YYYY <= 2023), ]
+  d <- d[which(d$YYYY >= styear & d$YYYY <= enyear), ]
 
   head(d)
   d$yrmon <- paste0(d$YYYY, sprintf("%02.f", d$MM))
@@ -50,13 +57,16 @@ labs <- c(rep("Average monthly degree heating weeks", 2),
           "Puerto Rico", "USVI")
 indnames <- data.frame(matrix(labs, nrow = 3, byrow = T))
 
-inddata <- list(labels = indnames, indicators = inddats, datelist = datdata) #, ulim = ulidata, llim = llidata)
-class(inddata) <- "indicatordata"
+ind <- list(labels = indnames, indicators = inddats, datelist = datdata) #, ulim = ulidata, llim = llidata)
+class(ind) <- "indicatordata"
 
 # save and plot -----------------------------------
-save(inddata, file = "C:/Users/mandy.karnauskas/Desktop/Caribbean-ESR/indicator_objects/DegreeHeatingWeeks.RData")
 
-plotIndicatorTimeSeries(inddata, coltoplot = 1:2, plotrownum = 2, sublabel = T, dateformat = "%Y%m", yposadj = 0.7,
+save(ind, file = "../../indicator_objects/DegreeHeatingWeeks.RData")
+
+plotIndicatorTimeSeries(ind, coltoplot = 1:2, plotrownum = 2, sublabel = T, dateformat = "%Y%m", yposadj = 0.7,
                         sameYscale = T, type = "allLines")
 
+# END ---------------
+
 

indicator_processing/automated_download/chl_caribbean.R

@@ -1,20 +1,45 @@
+#
+# ocean productivity indicator 
+#
+# download Ocean Color data from ERDDAP 
+# monthly time step
+# 
+# best database as of Mar 1 2024: 
+#    https://coastwatch.pfeg.noaa.gov/erddap/griddap/pmlEsaCCI60OceanColorMonthly_Lon0360
+# can also search ERDDAP https://coastwatch.pfeg.noaa.gov/erddap/index.html
+#    for "ESA CCI Ocean Colour Product"
+#
+############################################################################################
 
 rm(list = ls())
+dev.off()
 
-# https://coastwatch.pfeg.noaa.gov/erddap/griddap/pmlEsaCCI50OceanColorMonthly_Lon0360.graph
+load("spec_file.RData")
 
 # download data ------------------------------------------------
-url <- "https://coastwatch.pfeg.noaa.gov/erddap/griddap/erdMH1chlamday.csv?chlorophyll[(2003-01-16T00:00:00Z):1:(2021-12-16T00:00:00Z)][(18.5):1:(17.5)][(-67.4):1:(-64.4)]"
+
+url <- paste0("https://coastwatch.pfeg.noaa.gov/erddap/griddap/pmlEsaCCI60OceanColorMonthly_Lon0360.csv?chlor_a", 
+              "[(1997-12-01T00:00:00Z):1:(", terminal_year, "-12-01T00:00:00Z)][(", 
+              max_lat, "):1:(", min_lat, ")][(", 360 + min_lon, "):1:(", 360 + max_lon, ")]")
+url
+
+# old data source: "https://coastwatch.pfeg.noaa.gov/erddap/griddap/erdMH1chlamday.csv?chlorophyll[(2003-01-16T00:00:00Z):1:(2021-12-16T00:00:00Z)][(18.5):1:(17.5)][(-67.4):1:(-64.4)]"
+
+options(download.file.method="libcurl")
+
 download.file(url = url, destfile = "chl.csv")
+
 labs <- read.table("chl.csv", sep = ",", header = T, skip = 0)
 chl <- read.table("chl.csv", sep = ",", header = T, skip = 1)
 names(chl) <- names(labs)
 head(chl)
 
+file.remove("chl.csv")
+
 # calculate monthly mean and quantiles ------------------------------
-ind <- tapply(chl$chlorophyll, chl$time, mean, na.rm = T)
-uli <- tapply(chl$chlorophyll, chl$time, function(x) quantile(x, probs = c(0.975), na.rm = T))
-lli <- tapply(chl$chlorophyll, chl$time, function(x) quantile(x, probs = c(0.025), na.rm = T))
+ind <- tapply(chl$chlor_a, chl$time, mean, na.rm = T)
+uli <- tapply(chl$chlor_a, chl$time, function(x) quantile(x, probs = c(0.975), na.rm = T))
+lli <- tapply(chl$chlor_a, chl$time, function(x) quantile(x, probs = c(0.025), na.rm = T))
 
 # look for monthly trend -------------------------------------------
 mon <- strftime(names(ind), format = "%m")
@@ -24,7 +49,7 @@ summary(lm(ind ~ mon))    # significant seasonal trend
 # format into indicator object --------------------------------------
 datdata <- strftime(names(ind), format = "%m-%Y")
 inddata <- data.frame(ind, row.names = datdata)
-#ulidata <- data.frame(uli, row.names = datdata)
+#ulidata <- data.frame(uli, row.names = datdata)   # data are highly variable -- CIs very wide and not useful
 #llidata <- data.frame(lli, row.names = datdata)
 labs <- c("U.S. Caribbean primary productivity", "Mean Chlorophyll a Concentration (mg per m^3)", "monthly anomalies")
 indnames <- data.frame(matrix(labs, nrow = 3, byrow = T))
@@ -33,32 +58,33 @@ s <- list(labels = indnames, indicators = inddata, datelist = datdata) #, ulim =
 class(s) <- "indicatordata"
 
 # save and plot ------------------------------------------------------
-inddata <- s
-save(inddata, file = "C:/Users/mandy.karnauskas/Desktop/Caribbean-ESR/indicator_objects/carib_Chl.RData")
 
-plotIndicatorTimeSeries(s, dateformat = "%m-%Y", sublabel = T, trendAnalysis = T, anom = "mon", 
+ind <- s
+save(ind, file = "../../indicator_objects/carib_Chl.RData")
+
+plotIndicatorTimeSeries(ind, dateformat = "%m-%Y", sublabel = T, trendAnalysis = T, anom = "mon", 
                         widadj = 0.7, hgtadj = 0.8, type = "allLines")
 
 
  # plot spatiotemporal patterns ----------------------------------------
 
-nc <- nc_open("C:/Users/mandy.karnauskas/Downloads/erdMH1chlamday_a61e_882b_bec3.nc")
-v1 <- nc$var[[1]]
-chl <- ncvar_get(nc, v1)
-lon <- v1$dim[[1]]$vals 
-lat <- v1$dim[[2]]$vals 
-tim <- strftime(as.POSIXct(v1$dim[[3]]$vals, origin = "1970-01-01"), format = "%b %Y")
-nc_close(nc)
-
-par(mfrow = c(6, 12), mar = c(2, 1, 2, 1), mex = 0.75) 
-cols <- viridis::viridis(9)
-for (i in 1: 72) {
-  image(lon, lat[length(lat):1], chl[, length(lat):1, i], breaks = c(seq(0, 2, 0.25), 7), col = cols, main = tim[i])
-  }
-
-for (i in 73: 144) {
-  image(lon, lat[length(lat):1], chl[, length(lat):1, i], breaks = c(seq(0, 2, 0.25), 7), col = cols, main = tim[i])
-}
+#nc <- nc_open("C:/Users/mandy.karnauskas/Downloads/erdMH1chlamday_a61e_882b_bec3.nc")
+#v1 <- nc$var[[1]]
+#chl <- ncvar_get(nc, v1)
+#lon <- v1$dim[[1]]$vals 
+#lat <- v1$dim[[2]]$vals 
+#tim <- strftime(as.POSIXct(v1$dim[[3]]$vals, origin = "1970-01-01"), format = "%b %Y")
+#nc_close(nc)
+
+#par(mfrow = c(6, 12), mar = c(2, 1, 2, 1), mex = 0.75) 
+#cols <- viridis::viridis(9)
+#for (i in 1: 72) {
+#  image(lon, lat[length(lat):1], chl[, length(lat):1, i], breaks = c(seq(0, 2, 0.25), 7), col = cols, main = tim[i])
+#  }
+
+#for (i in 73: 144) {
+#  image(lon, lat[length(lat):1], chl[, length(lat):1, i], breaks = c(seq(0, 2, 0.25), 7), col = cols, main = tim[i])
+#}