tweaks to some indicators and plots

indicator_data/oceanEconomy_ENOW.csv

@@ -0,0 +1,61 @@
+geoid,geoName,year,sector,establishments,employment,wages,gdp
+72000,Puerto Rico,2019,Government,NA,ND,ND,ND
+72000,Puerto Rico,2020,Government,NA,ND,ND,ND
+72000,Puerto Rico,2021,Government,NA,ND,ND,ND
+72000,Puerto Rico,2019,Living Resources,10,134,5255000,14000000
+72000,Puerto Rico,2020,Living Resources,9,136,3605000,9700000
+72000,Puerto Rico,2021,Living Resources,10,147,5828000,15600000
+72000,Puerto Rico,2019,Marine Construction,14,152,3461000,5400000
+72000,Puerto Rico,2020,Marine Construction,15,174,4488000,6300000
+72000,Puerto Rico,2021,Marine Construction,15,355,7352000,9600000
+72000,Puerto Rico,2019,Marine Transportation,89,4594,155527000,290800000
+72000,Puerto Rico,2020,Marine Transportation,84,4278,132108000,320200000
+72000,Puerto Rico,2021,Marine Transportation,95,3841,144916000,352500000
+72000,Puerto Rico,2019,Ocean Economy,3072,53699,895344000,1639900000
+72000,Puerto Rico,2020,Ocean Economy,3115,54694,746315000,1437700000
+72000,Puerto Rico,2021,Ocean Economy,3607,53736,1014456000,1914300000
+72000,Puerto Rico,2019,Offshore Mineral Resources,17,198,4482000,8300000
+72000,Puerto Rico,2020,Offshore Mineral Resources,15,209,4594000,9000000
+72000,Puerto Rico,2021,Offshore Mineral Resources,20,296,7521000,14900000
+72000,Puerto Rico,2019,Research and Education,ND,ND,ND,ND
+72000,Puerto Rico,2020,Research and Education,ND,ND,ND,ND
+72000,Puerto Rico,2021,Research and Education,ND,ND,ND,ND
+72000,Puerto Rico,2019,Ship and Boat Building,NA,NA,NA,NA
+72000,Puerto Rico,2020,Ship and Boat Building,NA,NA,NA,NA
+72000,Puerto Rico,2021,Ship and Boat Building,NA,NA,NA,NA
+72000,Puerto Rico,2019,Tourism and Recreation,2941,48620,726575000,1321400000
+72000,Puerto Rico,2020,Tourism and Recreation,2990,49895,601483000,1092500000
+72000,Puerto Rico,2021,Tourism and Recreation,3464,49096,848788000,1521700000
+72000,Puerto Rico,2019,Utilities,9,233,18349776,95631717
+72000,Puerto Rico,2020,Utilities,9,258,18740176,88717101
+72000,Puerto Rico,2021,Utilities,14,252,19083604,106510188
+78000,U.S. Virgin Islands,2019,Government,NA,128,5896526,8425316
+78000,U.S. Virgin Islands,2020,Government,NA,148,6033590,8667984
+78000,U.S. Virgin Islands,2021,Government,NA,104,6243744,9293810
+78000,U.S. Virgin Islands,2019,Living Resources,NA,NA,NA,NA
+78000,U.S. Virgin Islands,2020,Living Resources,NA,NA,NA,NA
+78000,U.S. Virgin Islands,2021,Living Resources,NA,NA,NA,NA
+78000,U.S. Virgin Islands,2019,Marine Construction,3,62,2494000,4117341
+78000,U.S. Virgin Islands,2020,Marine Construction,4,19,587000,958924
+78000,U.S. Virgin Islands,2021,Marine Construction,NA,NA,NA,NA
+78000,U.S. Virgin Islands,2019,Marine Transportation,14,410,13077000,28298204
+78000,U.S. Virgin Islands,2020,Marine Transportation,16,474,16246000,28739403
+78000,U.S. Virgin Islands,2021,Marine Transportation,17,427,14954000,29612398
+78000,U.S. Virgin Islands,2019,Ocean Economy,300,5882,174698000,443380724
+78000,U.S. Virgin Islands,2020,Ocean Economy,309,6286,142127000,311563380
+78000,U.S. Virgin Islands,2021,Ocean Economy,303,5372,166226000,467645267
+78000,U.S. Virgin Islands,2019,Offshore Mineral Resources,NA,NA,NA,NA
+78000,U.S. Virgin Islands,2020,Offshore Mineral Resources,NA,NA,NA,NA
+78000,U.S. Virgin Islands,2021,Offshore Mineral Resources,NA,NA,NA,NA
+78000,U.S. Virgin Islands,2019,Research and Education,NA,62,1953146,2809067
+78000,U.S. Virgin Islands,2020,Research and Education,NA,178,5350517,6840404
+78000,U.S. Virgin Islands,2021,Research and Education,ND,ND,ND,ND
+78000,U.S. Virgin Islands,2019,Ship and Boat Building,NA,NA,NA,NA
+78000,U.S. Virgin Islands,2020,Ship and Boat Building,NA,NA,NA,NA
+78000,U.S. Virgin Islands,2021,Ship and Boat Building,NA,NA,NA,NA
+78000,U.S. Virgin Islands,2019,Tourism and Recreation,279,5315,152024000,296169289
+78000,U.S. Virgin Islands,2020,Tourism and Recreation,285,5713,121293000,216980238
+78000,U.S. Virgin Islands,2021,Tourism and Recreation,280,4836,145723000,310810246
+78000,U.S. Virgin Islands,2019,Utilities,ND,ND,ND,ND
+78000,U.S. Virgin Islands,2020,Utilities,ND,ND,ND,ND
+78000,U.S. Virgin Islands,2021,Utilities,ND,ND,ND,ND

indicator_processing/PlotAllIndicators.R

@@ -401,3 +401,33 @@ png(filename = "indicator_plots/total_landings_plot_TAP.png", width = 9, height
 load("indicator_objects/total_landings.RData")
 plotIndicatorTimeSeries(ind, coltoplot = 1:9, plotrownum = 3, plotcolnum = 3, trendAnalysis = TRUE, dateformat = "%b%Y", sublabel = TRUE, widadj = 0.9, hgtadj = 0.5, anom = "none", yposadj = 1, cex.axis = 1.2, cex.lab = 1.4, cex.main = 1.3)
 dev.off()
+
+
+
+png(filename = "indicator_plots/RVC_PR_plot_TAP.png", width = 7, height = 6, units = "in", res = 300)
+load("indicator_objects/RVC_PR.RData") 
+plotIndicatorTimeSeries(ind, coltoplot = 1:6, plotrownum = 3, plotcolnum = 2, trendAnalysis = TRUE, dateformat = "%b%Y", sublabel = TRUE, widadj = 1, hgtadj = 0.5, anom = "none", yposadj = 1, cex.axis = 1.2, cex.lab = 1.4, cex.main = 1.3)
+dev.off()
+
+png(filename = "indicator_plots/RVC_STSJ_plot_TAP.png", width = 7, height = 6, units = "in", res = 300)
+load("indicator_objects/RVC_STSJ.RData") 
+plotIndicatorTimeSeries(ind, coltoplot = 1:6, plotrownum = 3, plotcolnum = 2, trendAnalysis = TRUE, dateformat = "%b%Y", sublabel = TRUE, widadj = 1, hgtadj = 0.4, anom = "none", yposadj = 1, cex.axis = 1.2, cex.lab = 1.4, cex.main = 1.3)
+dev.off()
+
+png(filename = "indicator_plots/RVC_STX_plot_TAP.png", width = 7, height = 6, units = "in", res = 300)
+load("indicator_objects/RVC_STX.RData") 
+plotIndicatorTimeSeries(ind, coltoplot = 1:6, plotrownum = 3, plotcolnum = 2, trendAnalysis = TRUE, dateformat = "%b%Y", sublabel = TRUE, widadj = 1, hgtadj = 0.5, anom = "none", yposadj = 1, cex.axis = 1.2, cex.lab = 1.4, cex.main = 1.3)
+dev.off()
+
+
+png(filename = "indicator_plots/FRsection_plot_TAP.png", width = 10, height = 4, units = "in", res = 300)
+load("indicator_objects/FRsection.RData") 
+plotIndicatorTimeSeries(ind,coltoplot = 1, plotrownum = 1, plotcolnum = 1, trendAnalysis = TRUE, dateformat = "%b%Y", sublabel = F, widadj = 1, hgtadj = 1, anom = "none", yposadj = 1, cex.lab = 1, type = "allLines")
+dev.off()
+
+
+
+png(filename = "indicator_plots/oceanNAICS_plot_TAP.png", width = 6, height = 7, units = "in", res = 300)
+load("indicator_objects/oceanNAICS.RData") 
+plotIndicatorTimeSeries(ind, coltoplot = 1:8, plotrownum = 4, plotcolnum = 2, trendAnalysis = TRUE, dateformat = "%b%Y", sublabel = TRUE, widadj = 5, hgtadj = 1, anom = "none", yposadj = 1.2, cex.axis = 1, cex.lab = 1.1, cex.main = 1.3)
+dev.off()

indicator_processing/non_automated/RVC.R

@@ -41,6 +41,11 @@ enyear = 2023 #to update the plots when new data are added to the online portal,
 
 # Define the regions and corresponding file paths
 regions <- c("prico", "sttstj", "stx")
+
+
+# Run this code if you need to pull more data from the server
+##########
+
 data_2001_2021_paths <- c("indicator_data/RVC/prico_2001_2021_calibrated.rds", 
                           "indicator_data/RVC/sttstj_2001_2021_calibrated.rds", 
                           "indicator_data/RVC/stx_2001_2021_calibrated.rds")
@@ -88,7 +93,8 @@ for (i in seq_along(regions)) {
   saveRDS(combined_data, paste0("indicator_data/RVC/combined_", region, "_2001_2023.rds"))
 }
 
-
+#########
+# Start here if you don't need to pull new data from the server
 
 prico = readRDS("indicator_data/RVC/combined_prico_2001_2023.rds")
 sttstj = readRDS("indicator_data/RVC/combined_sttstj_2001_2023.rds")
@@ -119,7 +125,7 @@ for(j in regions) {
 
   # 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
@@ -133,14 +139,15 @@ for(j in regions) {
     # format indicator object -----------------------------
 
     datdata <- species_data$YEAR
-    inddata <- data.frame(species_data$density)
-    vardata <- data.frame(species_data$var)
-    s1 <- cbind(datdata, inddata)
-    s2 <- cbind(s1, vardata)
+    inddata <- data.frame(density = species_data$density)
+    vardata <- data.frame(var = species_data$var)
+    ulidata <- data.frame(uli = species_data$density + sqrt(species_data$var))
+    llidata <- data.frame(lli = species_data$density - sqrt(species_data$var))
+    s <- cbind(datdata, inddata, vardata, ulidata, llidata)
 
     # save -----------------------------------------
 
-    save(s2, file = paste("indicator_data/RVC/RUVdensity_", j, "_", i, ".RData", sep = ""))
+    save(s, file = paste("indicator_data/RVC/RUVdensity_", j, "_", i, ".RData", sep = ""))
   }
 }
 
@@ -177,46 +184,97 @@ for (file in files) {
 
 # 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))
+inddata <- data.frame(cbind(RUVdensity_PRICO_BAL_VETU$density,
+                            RUVdensity_PRICO_EPI_GUTT$density,
+                            RUVdensity_PRICO_LUT_ANAL$density,
+                            RUVdensity_PRICO_OCY_CHRY$density,
+                            RUVdensity_PRICO_SPA_AURO$density,
+                            RUVdensity_PRICO_SPA_VIRI$density))
+ulidata <- data.frame(cbind(RUVdensity_PRICO_BAL_VETU$uli,
+                            RUVdensity_PRICO_EPI_GUTT$uli,
+                            RUVdensity_PRICO_LUT_ANAL$uli,
+                            RUVdensity_PRICO_OCY_CHRY$uli,
+                            RUVdensity_PRICO_SPA_AURO$uli,
+                            RUVdensity_PRICO_SPA_VIRI$uli))
+llidata <- data.frame(cbind(RUVdensity_PRICO_BAL_VETU$lli,
+                            RUVdensity_PRICO_EPI_GUTT$lli,
+                            RUVdensity_PRICO_LUT_ANAL$lli,
+                            RUVdensity_PRICO_OCY_CHRY$lli,
+                            RUVdensity_PRICO_SPA_AURO$lli,
+                            RUVdensity_PRICO_SPA_VIRI$lli))
 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)
+inddata <- list(labels = indnames, indicators = inddata, datelist = datdata, ulim = ulidata, llim = llidata)
 class(inddata) <- "indicatordata"
 ind <- inddata
 save(ind, 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))
+inddata <- data.frame(cbind(RUVdensity_STTSTJ_BAL_VETU$density,
+                            RUVdensity_STTSTJ_EPI_GUTT$density,
+                            RUVdensity_STTSTJ_LUT_ANAL$density,
+                            RUVdensity_STTSTJ_OCY_CHRY$density,
+                            RUVdensity_STTSTJ_SPA_AURO$density,
+                            RUVdensity_STTSTJ_SPA_VIRI$density))
+ulidata <- data.frame(cbind(RUVdensity_STTSTJ_BAL_VETU$uli,
+                            RUVdensity_STTSTJ_EPI_GUTT$uli,
+                            RUVdensity_STTSTJ_LUT_ANAL$uli,
+                            RUVdensity_STTSTJ_OCY_CHRY$uli,
+                            RUVdensity_STTSTJ_SPA_AURO$uli,
+                            RUVdensity_STTSTJ_SPA_VIRI$uli))
+llidata <- data.frame(cbind(RUVdensity_STTSTJ_BAL_VETU$lli,
+                            RUVdensity_STTSTJ_EPI_GUTT$lli,
+                            RUVdensity_STTSTJ_LUT_ANAL$lli,
+                            RUVdensity_STTSTJ_OCY_CHRY$lli,
+                            RUVdensity_STTSTJ_SPA_AURO$lli,
+                            RUVdensity_STTSTJ_SPA_VIRI$lli))
 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)
+inddata <- list(labels = indnames, indicators = inddata, datelist = datdata, ulim = ulidata, llim = llidata)
 class(inddata) <- "indicatordata"
 ind <- inddata
 save(ind, 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))
+inddata <- data.frame(cbind(RUVdensity_STX_BAL_VETU$density,
+                            RUVdensity_STX_EPI_GUTT$density,
+                            RUVdensity_STX_LUT_ANAL$density,
+                            RUVdensity_STX_OCY_CHRY$density,
+                            RUVdensity_STX_SPA_AURO$density,
+                            RUVdensity_STX_SPA_VIRI$density))
+ulidata <- data.frame(cbind(RUVdensity_STX_BAL_VETU$uli,
+                            RUVdensity_STX_EPI_GUTT$uli,
+                            RUVdensity_STX_LUT_ANAL$uli,
+                            RUVdensity_STX_OCY_CHRY$uli,
+                            RUVdensity_STX_SPA_AURO$uli,
+                            RUVdensity_STX_SPA_VIRI$uli))
+llidata <- data.frame(cbind(RUVdensity_STX_BAL_VETU$lli,
+                            RUVdensity_STX_EPI_GUTT$lli,
+                            RUVdensity_STX_LUT_ANAL$lli,
+                            RUVdensity_STX_OCY_CHRY$lli,
+                            RUVdensity_STX_SPA_AURO$lli,
+                            RUVdensity_STX_SPA_VIRI$lli))
 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)
+inddata <- list(labels = indnames, indicators = inddata, datelist = datdata, ulim = ulidata, llim = llidata)
 class(inddata) <- "indicatordata"
 ind <- inddata
 save(ind, file = "indicator_objects/RVC_STX.RData")
@@ -233,7 +291,4 @@ plotIndicatorTimeSeries(ind, coltoplot = 1:6, trendAnalysis = T, sublabel = T)
 
 
 
-### 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.
 

indicator_processing/non_automated/oceanNAICS.R

@@ -1,4 +1,76 @@
-## Ocean economy indicators from NAICS
+## Ocean economy indicators from NAICS. There is an API, but it is not working. Instead, the data were downloaded manually from here: https://coast.noaa.gov/quickreport/#/index.html
+
+
+rm(list = ls())
+dev.off()
+
+
+library(plotTimeSeries)
+
+
+dat = read.csv("indicator_data/oceanEconomy_ENOW.csv")
+head(dat)
+
+yrs = as.integer(dat$year)
+
+styear = min(yrs)
+enyear = max(yrs)
+
+PR = subset(dat, dat$geoName == "Puerto Rico" & dat$sector == "Ocean Economy")
+USVI = subset(dat, dat$geoName == "U.S. Virgin Islands" & dat$sector == "Ocean Economy")
+
+GDP.PR = as.numeric(PR$gdp)/1000000000
+est.PR = as.numeric(PR$establishments)
+emp.PR = as.numeric(PR$employment)/1000
+wag.PR = as.numeric(PR$wages)/1000000000
+
+GDP.USVI = as.numeric(USVI$gdp)/1000000000
+est.USVI = as.numeric(USVI$establishments)
+emp.USVI = as.numeric(USVI$employment)/1000
+wag.USVI = as.numeric(USVI$wages)/1000000000
+
+
+
+# save as indicator object ----------------------
+datdata <- styear:enyear
+inddata <- data.frame(cbind(GDP.PR,GDP.USVI,est.PR,est.USVI,emp.PR,emp.USVI,wag.PR,wag.USVI))
+labs <- c("Gross domestic product" , "Billions of dollars", "Puerto Rico",
+          "Gross domestic product" , "Billions of dollars", "USVI",
+          "Ocean economy establishments" , "Number of places of work", "Puerto Rico",
+          "Ocean economy establishments" , "Number of places of work", "USVI",
+          "Ocean economy employees" , "People (thousands)" , "Puerto Rico",
+          "Ocean economy employees" , "People (thousands)" , "USVI",
+          "Ocean economy wages" , "Billions of dollars", "Puerto Rico",
+          "Ocean economy wages" , "Billions of dollars", "USVI")
+indnames <- data.frame(matrix(labs, nrow = 3, byrow = F))
+inddata <- list(labels = indnames, indicators = inddata, datelist = datdata)
+class(inddata) <- "indicatordata"
+
+# plot and save ----------------------------------
+ind <- inddata
+plotIndicatorTimeSeries(ind, plotrownum = 4, coltoplot = 1:8, sublabel = TRUE, trendAnalysis = T)
+
+save(ind, file = "indicator_objects/oceanNAICS.RData")
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+########################################################################################
+
+##### Below code is old, do not use ####################################################
+
+########################################################################################
 
 
 #### ************ As of 7/5/24, can't get the API for this working. For now, will just plot the data Seann uploaded to Github. At the end of this script is some code to work with the API but it isn't yet working.