finish P:D and Lmax indicators

indicator_processing/CalculateAllIndicators.R

@@ -55,9 +55,7 @@ source("indicator_processing/fishery_dependent/INDICATOR_disturbance.R")      #
 source("indicator_processing/fishery_dependent/INDICATOR_fishery_indicators_PR.R")      # % revenue, pelagic:demersal ratio, and Lmax indicators for Puerto Rico
 source("indicator_processing/fishery_dependent/INDICATOR_fishery_indicators_STT.R")     # % revenue, pelagic:demersal ratio, and Lmax indicators for STT/STJ
 source("indicator_processing/fishery_dependent/INDICATOR_fishery_indicators_STX.R")     # % revenue, pelagic:demersal ratio, and Lmax indicators for STX 
-
-
-
+source("indicator_processing/fishery_dependent/compile_indicators.R")         # compile Lmax and P:D indicators for all islands
 
 
 ###############################################################

indicator_processing/fishery_dependent/INDICATOR_fishery_indicators_PR.R

@@ -8,6 +8,7 @@ rm(list = ls())
 library(maps)
 library(plotTimeSeries)
 library(pals)
+library(Hmisc)
 
 load("indicator_processing/spec_file.RData")
 
@@ -200,6 +201,7 @@ rownames(tab2)[nsp] <- "other"
 tabr <- apply(tab2, 2, function(x) x/sum(x))
 
 # plot of proportion revenue by year and species group ------------------------
+
 colgd <- read.csv("indicator_processing/fishery_dependent/cols.csv", header = F) 
 
 barplot(tabr, col = as.character(colgd$V2[match(rownames(tabr), colgd$V1)]), 

indicator_processing/fishery_dependent/INDICATOR_fishery_indicators_STT.R

@@ -8,6 +8,7 @@ rm(list = ls())
 library(maps)
 library(plotTimeSeries)
 library(pals)
+library(Hmisc)
 
 load("indicator_processing/spec_file.RData")
 

indicator_processing/fishery_dependent/INDICATOR_fishery_indicators_STX.R

@@ -8,6 +8,7 @@ rm(list = ls())
 library(maps)
 library(plotTimeSeries)
 library(pals)
+library(Hmisc)
 
 load("indicator_processing/spec_file.RData")
 
@@ -175,7 +176,7 @@ barplot(tabr, col = as.character(colgd$V2[match(rownames(tabr), colgd$V1)]),
 
 plot(tab3, tabr)
 
-png(filename="C:/Users/mandy.karnauskas/Desktop/Caribbean-ESR/indicator_plots/per_landings_STX.png", 
+png(filename="indicator_plots/per_landings_STX.png", 
     units="in", width=7, height=4.5, pointsize=12, res=72*10)
 
 barplot(tabr, col = as.character(colgd$V2[match(rownames(tabr), colgd$V1)]), 

indicator_processing/fishery_dependent/compile_indicators.R

@@ -0,0 +1,111 @@
+# M. Karnauskas 10/19/2023
+# compile pelagic:demersal ratio and Lmax indicators
+# uses logbook data for PR and USVI 
+# need to run first: 
+#   INDICATOR_fishery_indicators_PR.R
+#   INDICATOR_fishery_indicators_STT.R
+#   INDICATOR_fishery_indicators_STX.R
+# then compile outputs here
+
+
+# specification file and libraries -----------------------------
+rm(list = ls())
+
+library(maps)
+library(plotTimeSeries)
+
+dir("indicator_data/fish-dep-indicators/")
+
+load("indicator_data/fish-dep-indicators/Lmax_PR.RData")
+pr <- findat
+
+load("indicator_data/fish-dep-indicators/Lmax_STT.RData")
+st <- findat
+
+load("indicator_data/fish-dep-indicators/Lmax_STX.RData")
+sx <- findat
+
+yrs <- (min(c(pr$V1, st$yrs, sx$yrs))) : (max(c(pr$V1, st$yrs, sx$yrs)))
+
+mat <- data.frame(matrix(data = NA, nrow = length(yrs), ncol = 3))
+rownames(mat) <- yrs
+mat
+
+mat[match(pr$V1, yrs), 1] <- pr$lmax
+mat[match(st$yrs, yrs), 2] <- st$lmax
+mat[match(sx$yrs, yrs), 3] <- sx$lmax
+
+datdata <- yrs
+inddata <- data.frame(mat)
+labs <- c("Average maximum length of species in catch", "length (cm)", "Puerto Rico", 
+          "Average maximum length of species in catch", "length (cm)", "St. Thomas and St. John",
+          "Average maximum length of species in catch", "length (cm)", "St. Croix")
+indnames <- data.frame(matrix(labs, nrow = 3, byrow = F))
+ind <- list(labels = indnames, indicators = inddata, datelist = datdata) #, ulim = ulidata, llim = llidata)
+
+class(ind) <- "indicatordata"
+plotIndicatorTimeSeries(ind, coltoplot = 1:3, plotrownum = 3, sublabel = T, sameYscale = F)
+
+save(ind, file = "indicator_objects/mean_Lmax.RData")
+
+Lmax <- mat
+
+
+# P:D ratios -----------------------------
+
+rm(list = ls()[-match(c("Lmax"), ls())])
+
+dir("indicator_data/fish-dep-indicators/")
+
+load("indicator_data/fish-dep-indicators/PDRatioPR.RData")
+pr <- pdrat
+
+load("indicator_data/fish-dep-indicators/PDRatioSTT.RData")
+st <- pdrat
+
+load("indicator_data/fish-dep-indicators/PDRatioSTX.RData")
+sx <- pdrat
+
+yrspr <- as.numeric(names(pr))
+yrsst <- as.numeric(names(st))
+yrssx <- as.numeric(names(sx))
+
+yrs <- (min(c(yrspr, yrsst, yrssx))) : (max(c(yrspr, yrsst, yrssx)))
+
+mat <- data.frame(matrix(data = NA, nrow = length(yrs), ncol = 3))
+rownames(mat) <- yrs
+mat
+
+mat[match(yrspr, yrs), 1] <- pr
+mat[match(yrsst, yrs), 2] <- st
+mat[match(yrssx, yrs), 3] <- sx
+
+datdata <- yrs
+inddata <- data.frame(mat)
+labs <- c("Ratio of pelagic to demersal catch", "landings ratio", "Puerto Rico", 
+          "Ratio of pelagic to demersal catch", "landings ratio", "St. Thomas and St. John",
+          "Ratio of pelagic to demersal catch", "landings ratio", "St. Croix")
+indnames <- data.frame(matrix(labs, nrow = 3, byrow = F))
+ind <- list(labels = indnames, indicators = inddata, datelist = datdata) #, ulim = ulidata, llim = llidata)
+
+class(ind) <- "indicatordata"
+plotIndicatorTimeSeries(ind, coltoplot = 1:3, plotrownum = 3, sublabel = T, sameYscale = F)
+
+save(ind, file = "indicator_objects/PD_ratio.RData")
+
+
+# compare indicators --------------------------------
+
+plot(Lmax[, 1], mat[, 1])
+tecor.test(Lmax[, 1], mat[, 1], na.action=na.omit)
+
+plot(Lmax[, 2], mat[, 2])
+cor.test(Lmax[, 2], mat[, 2], na.action=na.omit)
+
+plot(Lmax[, 3], mat[, 3])
+cor.test(Lmax[, 3], mat[, 3], na.action=na.omit)
+
+# indicators are very highly correlated 
+# use P:D ratio only as it is more direct indicator of what is going on
+
+