diff --git a/indicator_processing/fishery_dependent/INDICATOR_disturbance.R b/indicator_processing/fishery_dependent/INDICATOR_disturbance.R
index e7511f0..a7e9d13 100644
--- a/indicator_processing/fishery_dependent/INDICATOR_disturbance.R
+++ b/indicator_processing/fishery_dependent/INDICATOR_disturbance.R
@@ -373,7 +373,7 @@ class(s) <- "indicatordata"
 
 # plot and save -------------------------------------
 plotIndicatorTimeSeries(s, sublabel = T, coltoplot = 1:3, plotrownum = 3, yposadj = 1.2, type = "allLines", 
-                        sameYscale = T, trendAnalysis = F)
+                        sameYscale = F, trendAnalysis = F)
 
 inddata <- s
 save(inddata, file = "C:/Users/mandy.karnauskas/Desktop/Caribbean-ESR/indicator_objects/disturbance.RData")
diff --git a/indicator_processing/fishery_dependent/INDICATOR_fishery_indicators.R b/indicator_processing/fishery_dependent/INDICATOR_fishery_indicators.R
index e33d54d..cc73577 100644
--- a/indicator_processing/fishery_dependent/INDICATOR_fishery_indicators.R
+++ b/indicator_processing/fishery_dependent/INDICATOR_fishery_indicators.R
@@ -1,112 +1,96 @@
+# M. Karnauskas 10/19/2023
+# code for calculating pelagic:demersal ratio and Lmax indicators
+# uses logbook data for PR and USVI 
 
 rm(list = ls())
 
 library(pals)
+
+# define start and end years ---------------------------
+styear <- 1990
+enyear <- 2020
+
+# input data for Puerto Rico ---------------------------
 setwd("C:/Users/mandy.karnauskas/Desktop/Caribbean-ESR/indicator_processing/fishery_dependent/")
+dat <- read.csv("C:/Users/mandy.karnauskas/Desktop/CONFIDENTIAL/CaribbeanData/Jun2022/PR_landings_83_20_wSC_2005cor.csv")
+
+table(dat$YEAR_LANDED)
 
-d <- read.csv("C:/Users/mandy.karnauskas/Desktop/CONFIDENTIAL/CaribbeanData/Jun2022/PR_landings_83_20.csv")
-
-# take a look at data fields ----------------------------
-table(d$VESSEL, useNA = "always")
-table(d$YEAR_LANDED, useNA = "always")
-table(d$MONTH_LANDED, useNA = "always")
-table(d$DAY_LANDED, useNA = "always")
-table(d$FISHING_CENTER_ED, useNA = "always")
-table(d$FISHING_CENTER_NAME, useNA = "always")
-table(d$MUNICIPALITY, useNA = "always")
-table(d$AREA_FISHED1, useNA = "always")
-table(d$AREA_FISHED2, useNA = "always")
-table(d$AREA_FISHED3, useNA = "always")
-table(d$AREA_FISHED4, useNA = "always")
-table(d$FIN_GEAR_CODE, useNA = "always")
-table(d$FIN_GEAR_NAME, useNA = "always")
-table(d$PR_ID_CODE_ED, useNA = "always")
-table(d$NUMBER_OF_TRIPS_ED, useNA = "always")
-table(d$GEAR_QTY_ED, useNA = "always")
-table(d$GEAR_HOURS_ED, useNA = "always")
-hist(d$MAXIMUM_DEPTH_ED)
-hist(d$MINIMUM_DEPTH_ED)
-table(d$SPECIES_ITIS, useNA = "always")
-table(d$ITIS_COMMON_NAME, useNA = "always")
-table(d$ITIS_SCIENTIFIC_NAME, useNA = "always")
-hist(d$POUNDS_LANDED)
-hist(d$ADJUSTED_POUNDS)
-hist(d$VALUE_IN_DOLLARS)
-hist(d$PRICE_PER_LB)
-table(d$DISTANCE, useNA = "always")
-table(d$DISTANCE_DESCRIBE, useNA = "always")
-#table(d$TRIP_TICKET_NUMBER_ED, useNA = "always")
-
-# look at gear trends ---------------------------------
-tapply(d$ADJUSTED_POUNDS, d$FIN_GEAR_NAME, sum, na.rm = T)
-par(mar = c(5, 15, 2, 2))
-barplot(tapply(d$ADJUSTED_POUNDS, d$FIN_GEAR_NAME, sum, na.rm = T), horiz = T, las = 2)
-
-d$gear <- "other" 
-d$gear[grep("POTS", d$FIN_GEAR_NAME)] <- d$FIN_GEAR_NAME[grep("POTS", d$FIN_GEAR_NAME)]
-
-par(mar = c(4, 4, 1, 1))
-tab <- tapply(d$ADJUSTED_POUNDS, list(d$gear, d$YEAR_LANDED), sum, na.rm = T)
-barplot(tab, col = c(2, gray(0.1), gray(0.5), gray(0.8)), args.legend = c(x = "topright"), legend.text = rownames(tab))
-
-tab2 <- tab
-for (i in 1:ncol(tab)) {   tab2[, i] <- tab[, i] / sum(tab[, i])   }
-
-barplot(tab2, col = c(2, gray(0.2), gray(0.4), gray(0.6)))
-
-numtrap <- colSums(tab[2:4, ]) / (10^6)
-pertrap <- 1 - tab2[1, ]
-
-# format into indicator object --------------------------------------
-#datdata <- 1983:2020
-#inddata <- data.frame(cbind(numtrap, pertrap))
-#labs <- c("Number of traps", "millions of traps" , "Puerto Rico", "Proportion of effort from traps", "proportion", "Puerto Rico")
-#indnames <- data.frame(matrix(labs, nrow = 3, byrow = F))
-
-#s <- list(labels = indnames, indicators = inddata, datelist = datdata) #, ulim = ulidata, llim = llidata)
-#class(s) <- "indicatordata"
-
-#plotIndicatorTimeSeries(s, coltoplot = 1:2, plotrownum = 2, sublabel = T, outtype = "png")
+# subset years------------------------------------------
+
+d <- dat[which(dat$YEAR_LANDED >= styear & dat$YEAR_LANDED <= enyear), ]
+table(d$YEAR_LANDED)
 
 # look at main species landed --------------------------------
 tab <- sort(tapply(d$ADJUSTED_POUNDS, d$ITIS_COMMON_NAME, sum, na.rm = T), decreasing = T)
 par(mar = c(15, 5, 2, 2))
 barplot(tab, las = 2)
-barplot(tab[1:100], las = 2)
 barplot(tab[1:50], las = 2)
+barplot(tab[1:25], las = 2)
 
 tab <- tapply(d$ADJUSTED_POUNDS, list(d$ITIS_COMMON_NAME, d$YEAR_LANDED), sum, na.rm = T)
 tab <- tab[order(rowSums(tab, na.rm = T), decreasing = T), ]
 
 par(mar = c(4, 4, 1, 1)+1)
-matplot(1983:2020, t(tab[1:10, ]), type = "l", col = 1:10, lty = c(1, 1, 1, 1, 1, 2, 2, 2, 2, 2), lwd = 2)
+matplot(styear:enyear, t(tab[1:10, ]), type = "l", col = 1:10, lty = c(1, 1, 1, 1, 1, 2, 2, 2, 2, 2), lwd = 2)
 legend("topright", rownames(tab)[1:10], col = 1:10, lty = c(1, 1, 1, 1, 1, 2, 2, 2, 2, 2), lwd = 2)
 
 tab2 <- tab[1:10, ]
 for (i in 1:ncol(tab)) {  tab2[, i] <- tab[1:10, i] / sum(tab[1:10, i], na.rm = T)   }
 
-barplot(tab2, col = glasbey(10), xlim = c(0, 56), legend.text = rownames(tab2), args.legend = c(x = "right"))
+barplot(tab2, col = glasbey(10), xlim = c(0, 56), legend.text = rownames(tab2), args.legend = c(x = "right"), las = 2)
 
-rownames(tab)[4]
-plot(1983:2020, tab[4, ]/colSums(tab, na.rm = T), main = "proportion of unidentified landings by year", type = "b")
+j <- grep("FISHES,BONY,UNSPEC", rownames(tab)); rownames(tab)[j]
+plot(styear:enyear, tab[j, ]/colSums(tab, na.rm = T), main = "proportion of unidentified landings by year", type = "b")
 abline(h = 0.02)
 abline(v = 2005) # 2000 and on is 'data rich' period
 abline(h = 0.1)
 abline(v = 2000)
 
+# redo starting year cutoff based on data avail --------------------------
+styear <- 2000
+d <- d[which(d$YEAR_LANDED >= styear), ]
+
+# add missing scientific names
+
+d$ITIS_SCIENTIFIC_NAME <- as.character(d$ITIS_SCIENTIFIC_NAME)
+lis <- as.data.frame(unique(d$ITIS_SCIENTIFIC_NAME))
+names(lis) = "sci"
+lis$com <- NA
+head(lis)
+lis <- lis[1:(nrow(lis)-1), ]
+lis
+
+for (i in 1:nrow(lis))  { 
+  n <- sort(table(d$ITIS_COMMON_NAME[which(d$ITIS_SCIENTIFIC_NAME == lis$sci[i])]), decreasing = T)
+  lis$com[i] <- names(n[1])    }
+
+table(d$YEAR_LANDED[which(is.na(d$ITIS_SCIENTIFIC_NAME))])
+
+mis <- which(is.na(d$ITIS_SCIENTIFIC_NAME))
+missci <- as.character(lis$sci[match(d$ITIS_COMMON_NAME[mis], lis$com)])
+d$ITIS_SCIENTIFIC_NAME[which(is.na(d$ITIS_SCIENTIFIC_NAME))] <- missci
+dim(table(d$ITIS_SCIENTIFIC_NAME))
+
+mis <- which(is.na(d$ITIS_SCIENTIFIC_NAME))
+d[mis, ]   # still some missing but not going to fix until new data pull
+
 # pull in reference file and merge ------------------------------------
 
 ref <- read.csv("spp_ref_manualEdit.csv")
 head(ref)
 head(d)
+tail(d)
 
 hist(ref$Lmax)
 table(cut(ref$Lmax, breaks = c(0, 40, 60, 100, 200, 2000)))
 ref$Lmax_cat <- cut(ref$Lmax, breaks = c(0, 40, 60, 100, 200, 2000))
 
 db <- merge(d, ref, by.x = "ITIS_SCIENTIFIC_NAME", by.y = "SCIname", all.x = TRUE)
+dim(ref)
 dim(d)
 dim(db)
+head(db)
 
 tab <- tapply(db$ADJUSTED_POUNDS, list(db$famcode, db$YEAR_LANDED), sum, na.rm = T)
 tab <- tab[order(rowSums(tab, na.rm = T), decreasing = T), ]
@@ -115,13 +99,16 @@ head(tab, 15)
 nsp <- 15
 cols <- glasbey(nsp)
 
+# plot top species groups --------------------------------
+
 par(mar = c(4, 4, 1, 1)+1)
-matplot(1983:2020, t(tab[1:nsp, ]), type = "l", col = cols, lty = rep(1:3, (nsp/3)), lwd = 2)
+matplot(styear:enyear, t(tab[1:nsp, ]), type = "l", col = cols, lty = rep(1:3, (nsp/3)), lwd = 2)
 legend("topright", rownames(tab)[1:nsp], col = cols, lty = rep(1:3, (nsp/3)), lwd = 2)
 
 tab2 <- tab[1:(nsp-1), ]
 tab2 <- rbind(tab2, colSums(tab[nsp:nrow(tab), ], na.rm = T))
 rownames(tab2)[nsp] <- "other"
+tab2
 
 tab3 <- tab2
 for (i in 1:ncol(tab2)) {  tab3[, i] <- tab2[, i] / sum(tab2[, i], na.rm = T)   }
@@ -129,15 +116,16 @@ for (i in 1:ncol(tab2)) {  tab3[, i] <- tab2[, i] / sum(tab2[, i], na.rm = T)
 cols[which(rownames(tab3) == "UNID")] <- "white"
 barplot(tab3, col = cols, xlim = c(0, 56), legend.text = rownames(tab2), args.legend = c(x = "right"), las = 2)
 
-#png(filename="PR_proportion_landings.png", units="in", width=8, height=5, pointsize=12, res=72*10)
-
-barplot(tab3[, 23:38], col = cols, xlim = c(0, 23), legend.text = rownames(tab2), args.legend = c(x = "right", bty = "n"), las = 2, space = 0, border = NA, 
+barplot(tab3[, which(colnames(tab3) == styear):ncol(tab3)], col = cols, xlim = c(0, 23), legend.text = rownames(tab2), args.legend = c(x = "right", bty = "n"), las = 2, space = 0, border = NA, 
         xlab = "", ylab = "proportion of catch", main = "Puerto Rico landings composition by year                                ")
 abline(h = 0)
 
+#png(filename="PR_proportion_landings.png", units="in", width=8, height=5, pointsize=12, res=72*10)
+
 #dev.off()
 
-table(db$ITIS_COMMON_NAME == db$COMname)
+
+table(as.character(db$ITIS_COMMON_NAME) == as.character(db$COMname))
 length(which(db$ITIS_COMMON_NAME == "FISHES,BONY,UNSPECIFIED"))
 table(db$PD, useNA = "always")
 table(db$Lmax, useNA = "always")
@@ -156,16 +144,16 @@ dbf$PD2[grep("pelagic", dbf$PD)] <- "pelagic"
 table(dbf$PD, dbf$PD2)
 
 pd <- tapply(dbf$ADJUSTED_POUNDS, list(dbf$YEAR_LANDED, dbf$PD2), sum, na.rm = T)
-pd2 <- pd[18:38, ]
+pd
+matplot(pd)
 pdrat <- pd[, 2] / pd[, 1]
-plot(1983:2020, pdrat, type = "b")
+plot(styear:enyear, pdrat, type = "b")
 
-pdrat2 <- pdrat[18:38]
-save(pdrat2, file = "C:/Users/mandy.karnauskas/Desktop/Caribbean-ESR/indicator_data/PDRatioPR.RData")
+save(pdrat, file = "C:/Users/mandy.karnauskas/Desktop/Caribbean-ESR/indicator_data/PDRatioPR.RData")
 
 # make indicator object and plot P:D ratio ------------------
-datdata <- 2000:2020
-inddata <- data.frame(cbind(pd2[, 2]/10^6, pd2[, 1]/10^6, pdrat2))
+datdata <- styear:enyear
+inddata <- data.frame(cbind(pd[, 2]/10^6, pd[, 1]/10^6, pdrat))
 labs <- c("Total pelagic landings", "millions of pounds" , "Puerto Rico", 
           "Total demersal landings", "millions of pounds", "Puerto Rico", 
           "Pelagic:demersal ratio", "ratio of landings", "Puerto Rico")
@@ -174,13 +162,17 @@ 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, coltoplot = 3, plotrownum = 1, sublabel = T, widadj = 1.2, outtype = "png")
+plotIndicatorTimeSeries(s, coltoplot = 1:3,plotrownum = 3, sublabel = T, widadj = 1.2, trendAnalysis = F)
+
+dev.off()
+
+#######################  END P:D ratio  ################################
 
 # Lmax calculations -----------------------------------------
 
 tab <- tapply(dbf$ADJUSTED_POUNDS, list(dbf$ITIS_SCIENTIFIC_NAME, dbf$YEAR_LANDED), sum, na.rm = T)
 tab[is.na(tab)] <- 0
+tab
 
 splis <- data.frame(rownames(tab))
 names(splis) <- "SCIname"
@@ -197,7 +189,6 @@ Lmax5 <- colSums(tab[which(splisref$Lmax_cat == llis[5]), ])
 #Lmax6 <- colSums(tab[which(splisref$Lmax_cat == llis[6]), ])
 
 Lmaxcl <- cbind(Lmax1, Lmax2, Lmax3, Lmax4, Lmax5) #, Lmax6)
-Lmaxcl <- Lmaxcl[18:38, ]
 
 Lmaxcl2 <- Lmaxcl
 for (i in 1:nrow(Lmaxcl)) {  Lmaxcl2[i, ] <- Lmaxcl[i, ] / sum(Lmaxcl[i, ], na.rm = T)   }
@@ -208,7 +199,7 @@ barplot(t(Lmaxcl2), col = 1:6, las = 2, main = "Distribution of catch in Lmax si
 legend = c("<40 cm", "40-60 cm", "60-100 cm", "100-200 cm", ">200 cm"), args.legend = c(x = "right", bty = "n"))
 
 # format indicator objects and plot ---------------------------------------
-datdata <- 2000:2020
+datdata <- styear:enyear
 inddata <- data.frame(Lmaxcl/10^6)
 labs <- c(rep("Total landings in Lmax class", 5), 
           rep("millions of pounds", 5),
@@ -216,7 +207,7 @@ labs <- c(rep("Total landings in Lmax class", 5),
 indnames <- data.frame(matrix(labs, nrow = 3, byrow = T))
 s <- list(labels = indnames, indicators = inddata, datelist = datdata) #, ulim = ulidata, llim = llidata)
 class(s) <- "indicatordata"
-plotIndicatorTimeSeries(s, coltoplot = 1:5, plotrownum = 5, sublabel = T, widadj = 1.5) #, outtype = "png")
+plotIndicatorTimeSeries(s, coltoplot = 1:5, plotrownum = 5, sublabel = T, widadj = 1.5, trendAnalysis = F) #, outtype = "png")
 
 inddata <- data.frame(Lmaxcl2)
 labs <- c(rep("Proportion of landings in Lmax class",5), 
@@ -225,7 +216,7 @@ labs <- c(rep("Proportion of landings in Lmax class",5),
 indnames <- data.frame(matrix(labs, nrow = 3, byrow = T))
 s <- list(labels = indnames, indicators = inddata, datelist = datdata) #, ulim = ulidata, llim = llidata)
 class(s) <- "indicatordata"
-plotIndicatorTimeSeries(s, coltoplot = 1:5, plotrownum = 5, sublabel = T, widadj = 1.5, hgtadj = 0.6) # outtype = "png", sameYscale = F)
+plotIndicatorTimeSeries(s, coltoplot = 1:5, plotrownum = 5, sublabel = T, widadj = 1.5, hgtadj = 0.6, trendAnalysis = F) # outtype = "png", sameYscale = F)
 
 # understand what is driving the trends -------------------------------------------
 
@@ -235,7 +226,7 @@ splisref$COMname[which(splisref$Lmax_cat == "(100,200]")]
 sort(table(splisref$famcode[which(splisref$Lmax_cat == "(60,100]")]))
 sort(table(splisref$famcode[which(splisref$Lmax_cat == "(100,200]")]))
 
-splisref$recLand <- rowSums(tab[, 18:38])
+splisref$recLand <- rowSums(tab)
 
 plate <- splisref[which(splisref$Lmax_cat == "(60,100]"), ]
 head(plate[order(plate$recLand, decreasing = T), ], 15)
@@ -248,17 +239,17 @@ head(big[order(big$recLand, decreasing = T), ], 15)
 lmax <- rep(NA, ncol(tab))
 for (i in 1:ncol(tab))  {  lmax[i] <- weighted.mean(splisref$Lmax, tab[, i])  }
 
-plot(1983:2020, lmax)
-plot(2000:2020, lmax[18:38], type = "b", las = 2)
-axis(1, 2000:2020)
-out <- lm(lmax[18:38] ~ c(2000:2020))
+plot(styear:enyear, lmax)
+plot(styear:enyear, lmax, type = "b", las = 2)
+axis(1, styear:enyear)
+out <- lm(lmax ~ c(styear:enyear))
 summary(out)
 abline(out)
 
-plot(1983:2020, tapply(dbf$Lmax, dbf$YEAR_LANDED, mean, na.rm = T))
+plot(styear:enyear, tapply(dbf$Lmax, dbf$YEAR_LANDED, mean, na.rm = T))
 
-datdata <- 2000:2020
-inddata <- data.frame(lmax[18:38])
+datdata <- styear:enyear
+inddata <- data.frame(lmax)
 labs <- c("Average maximum length of catch", "length (cm)" , "Puerto Rico")
 indnames <- data.frame(matrix(labs, nrow = 3, byrow = F))
 s <- list(labels = indnames, indicators = inddata, datelist = datdata) #, ulim = ulidata, llim = llidata)
@@ -266,22 +257,25 @@ class(s) <- "indicatordata"
  
 plotIndicatorTimeSeries(s, coltoplot = 1, sublabel = T) #, outtype = "png")
 
+dev.off()
+
 # figure out what is going on in 2018 with spike in pelagics -------------------
 
-plot(1983:2020, pdrat, type = "b", las = 2)
-matplot(1983:2020, pd, axes = F, type = "b")
-axis(1, at = 1983:2020, las = 2)
+plot(styear:enyear, pdrat, type = "b", las = 2, pch = 19)
+matplot(styear:enyear, pd, axes = F, type = "b", pch = 19)
+axis(1, at = styear:enyear, las = 2)
 
 tabp <- tab[grep("pelagic", splisref$PD), ]
 tabp <- tabp[order(rowSums(tabp), decreasing = T), ]
 
-matplot(1983:2020, t(tabp[1:10, ]), type = "l", col = glasbey(10), lwd = 2, lty = 1)
+matplot(styear:enyear, t(tabp[1:10, ]), type = "l", col = glasbey(10), lwd = 2, lty = 1)
 legend("topright", rownames(tabp)[1:10], col = glasbey(10), lwd = 2, lty = 1)
+abline(v = c(2008, 2019))
 
 tabd <- tab[-grep("pelagic", splisref$PD), ]
 tabd <- tabd[order(rowSums(tabd), decreasing = T), ]
 
-matplot(1983:2020, t(tabd[1:10, ]), type = "b", col = glasbey(10), lwd = 2, lty = 1, pch = 19, las = 2)
+matplot(styear:enyear, t(tabd[1:10, ]), type = "b", col = glasbey(10), lwd = 2, lty = 1, pch = 19, las = 2)
 legend("topright", rownames(tabd)[1:10], col = glasbey(10), lwd = 2, lty = 1)
-
+abline(v = c(2008, 2019))
 
diff --git a/indicator_processing/fishery_dependent/INDICATOR_fishery_indicators_STT.R b/indicator_processing/fishery_dependent/INDICATOR_fishery_indicators_STT.R
index db08dec..c6870b9 100644
--- a/indicator_processing/fishery_dependent/INDICATOR_fishery_indicators_STT.R
+++ b/indicator_processing/fishery_dependent/INDICATOR_fishery_indicators_STT.R
@@ -1,10 +1,23 @@
+# M. Karnauskas 10/19/2023
+# code for calculating pelagic:demersal ratio and Lmax indicators
+# uses logbook data for PR and USVI 
 
 rm(list = ls())
-
 library(pals)
+
 setwd("C:/Users/mandy.karnauskas/Desktop/Caribbean-ESR/indicator_processing/fishery_dependent/")
 
-d <- read.csv("C:/Users/mandy.karnauskas/Desktop/CONFIDENTIAL/CaribbeanData/STT_landings.csv")
+dat <- read.csv("C:/Users/mandy.karnauskas/Desktop/CONFIDENTIAL/CaribbeanData/STT_landings.csv")
+
+# define start and end years ---------------------------
+styear <- 1990
+enyear <- 2020
+table(dat$TRIP_YEAR)
+
+# subset years------------------------------------------
+
+d <- dat[which(dat$TRIP_YEAR >= styear & dat$TRIP_YEAR <= enyear), ]
+table(d$TRIP_YEAR)
 
 # take a look at data fields ----------------------------
 
@@ -40,43 +53,26 @@ par(mar = c(4, 4, 1, 1))
 tab <- tapply(d$POUNDS_LANDED, list(d$gear, d$TRIP_YEAR), sum, na.rm = T)
 barplot(tab, col = c(2, gray(0.1), gray(0.5), gray(0.8)), args.legend = c(x = "topright"), legend.text = rownames(tab))
 
-tab2 <- tab
-for (i in 1:ncol(tab)) {   tab2[, i] <- tab[, i] / sum(tab[, i])   }
-
-barplot(tab2, col = c(2, gray(0.2), gray(0.4), gray(0.6)))
-
-numtrap <- tab[2,] / (10^6)
-pertrap <- 1 - tab2[1, ]
-
-# format into indicator object --------------------------------------
-#datdata <- 1974:2021
-#inddata <- data.frame(cbind(numtrap, pertrap))
-#labs <- c("Number of traps", "millions of traps" , "St. Thomas", "Proportion of effort from traps", "proportion", "St. Thomas")
-#indnames <- data.frame(matrix(labs, nrow = 3, byrow = F))
-
-#s <- list(labels = indnames, indicators = inddata, datelist = datdata) #, ulim = ulidata, llim = llidata)
-#class(s) <- "indicatordata"
-
-#plotIndicatorTimeSeries(s, coltoplot = 1:2, plotrownum = 2, sublabel = T) # outtype = "png")
-
 # look at main species landed --------------------------------
 tab <- sort(tapply(d$POUNDS_LANDED, d$SPECIES_NM, sum, na.rm = T), decreasing = T)
 par(mar = c(15, 5, 2, 2))
 barplot(tab, las = 2)
-barplot(tab[1:100], las = 2)
 barplot(tab[1:50], las = 2)
+barplot(tab[1:25], las = 2)
 
 tab <- tapply(d$POUNDS_LANDED, list(d$SPECIES_NM, d$TRIP_YEAR), sum, na.rm = T)
 tab <- tab[order(rowSums(tab, na.rm = T), decreasing = T), ]
 
 par(mar = c(4, 4, 1, 1)+1)
-matplot(1974:2021, t(tab[1:10, ]), type = "l", col = 1:10, lty = c(1, 1, 1, 1, 1, 2, 2, 2, 2, 2), lwd = 2)
+matplot(styear:enyear, t(tab[1:10, ]), type = "l", col = 1:10, lty = c(1, 1, 1, 1, 1, 2, 2, 2, 2, 2), lwd = 2)
 legend("topright", rownames(tab)[1:10], col = 1:10, lty = c(1, 1, 1, 1, 1, 2, 2, 2, 2, 2), lwd = 2)
 abline(v = 1999)   # 2000 is data-rich period
 
-d <- d[which(d$TRIP_YEAR >= 2010), ]  # cut out because very little data beforehand
+styear <- 2010
 
-yrs <- 2010:2021
+d <- d[which(d$TRIP_YEAR >= styear), ]  # cut out because very little data beforehand
+
+yrs <- styear:enyear
 
 tab <- tapply(d$POUNDS_LANDED, list(d$SPECIES_NM, d$TRIP_YEAR), sum, na.rm = T)
 tab <- tab[order(rowSums(tab, na.rm = T), decreasing = T), ]
@@ -103,6 +99,7 @@ ref$Lmax_cat <- cut(ref$Lmax, breaks = c(0, 40, 60, 100, 200, 2000))
 db <- merge(d, ref, by.x = "SPECIES_NM", by.y = "COMname", all.x = TRUE)
 dim(d)
 dim(db)
+which(is.na(db$SPECIES_NM))
 
 tab <- tapply(db$POUNDS_LANDED, list(db$famcode, db$TRIP_YEAR), sum, na.rm = T)
 tab <- tab[order(rowSums(tab, na.rm = T), decreasing = T), ]
@@ -172,30 +169,7 @@ class(s) <- "indicatordata"
 
 plotIndicatorTimeSeries(s, coltoplot = 1:3, plotrownum = 3, sublabel = T, trendAnalysis = F) #, outtype = "png")
 
-# combine with PR ------------------------
-
-load("C:/Users/mandy.karnauskas/Desktop/Caribbean-ESR/indicator_data/PDRatioPR.RData")
-pdrat
-length(pdrat)
-pdrat[1:11]
-pdrat2
-
-pdrat <- c(rep(NA, 10), pdrat[1:11])
-#pdrat2 <- c(pdrat2, NA)
-
-datdata <- 2000:2020
-inddata <- data.frame(cbind(pdrat2, pdrat))
-labs <- c("Pelagic to demersal ratio", "ratio of landings", "Puerto Rico", 
-          "Pelagic to demersal ratio", "ratio of landings", "St. Thomas and St. John")
-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, coltoplot = 1:2, plotrownum = 2, sublabel = T, sameYscale = T)
-plotIndicatorTimeSeries(s, coltoplot = 1:2, plotrownum = 2, sublabel = T, sameYscale = T, hgtadj = 0.8, outtype = "png")
-
+############################  END PD ratio  #######################################
 
 # Lmax calculations -----------------------------------------
 
@@ -237,7 +211,7 @@ labs <- c(rep("Total landings in Lmax class", 5),
 indnames <- data.frame(matrix(labs, nrow = 3, byrow = T))
 s <- list(labels = indnames, indicators = inddata, datelist = datdata) #, ulim = ulidata, llim = llidata)
 class(s) <- "indicatordata"
-plotIndicatorTimeSeries(s, coltoplot = 1:5, plotrownum = 5, sublabel = T, widadj = 1.5) #, outtype = "png")
+plotIndicatorTimeSeries(s, coltoplot = 1:5, plotrownum = 5, sublabel = T, widadj = 1.5, trendAnalysis = F) #, outtype = "png")
 
 inddata <- data.frame(Lmaxcl2)
 labs <- c(rep("Proportion of landings in Lmax class",5), 
@@ -246,7 +220,7 @@ labs <- c(rep("Proportion of landings in Lmax class",5),
 indnames <- data.frame(matrix(labs, nrow = 3, byrow = T))
 s <- list(labels = indnames, indicators = inddata, datelist = datdata) #, ulim = ulidata, llim = llidata)
 class(s) <- "indicatordata"
-plotIndicatorTimeSeries(s, coltoplot = 1:5, plotrownum = 5, sublabel = T, widadj = 1.5, hgtadj = 0.6) # outtype = "png", sameYscale = F)
+plotIndicatorTimeSeries(s, coltoplot = 1:5, plotrownum = 5, sublabel = T, widadj = 1.5, hgtadj = 0.6, trendAnalysis = F) # outtype = "png", sameYscale = F)
 
 # understand what is driving the trends -------------------------------------------
 
@@ -270,39 +244,30 @@ lmax <- rep(NA, ncol(tab))
 for (i in 1:ncol(tab))  {  lmax[i] <- weighted.mean(splisref$Lmax, tab[, i])  }
 
 plot(yrs, lmax)
-#plot(2000:2020, lmax[18:38], type = "b", las = 2)
-axis(1, 2000:2021)
+plot(yrs, lmax, type = "b", las = 2)
+axis(1, yrs)
 out <- lm(lmax ~ yrs)
 summary(out)
 abline(out)
 
-plot(yrs, tapply(dbf$Lmax, dbf$TRIP_YEAR, mean, na.rm = T))
-
-datdata <- yrs
-inddata <- data.frame(lmax)
-labs <- c("Average maximum length of catch", "length (cm)" , "Puerto Rico")
-indnames <- data.frame(matrix(labs, nrow = 3, byrow = F))
-s <- list(labels = indnames, indicators = inddata, datelist = datdata) #, ulim = ulidata, llim = llidata)
-class(s) <- "indicatordata"
-
-plotIndicatorTimeSeries(s, coltoplot = 1, sublabel = T) # outtype = "png")
+plot(yrs, tapply(dbf$Lmax, dbf$TRIP_YEAR, mean, na.rm = T), type = "l")
 
 # figure out what is going on in 2018 with spike in pelagics -------------------
 
-plot(1983:2020, pdrat, type = "b", las = 2)
-matplot(1983:2020, pd, axes = F, type = "b")
-axis(1, at = 1983:2020, las = 2)
+plot(yrs, pdrat, type = "b", las = 2)
+matplot(yrs, pd, axes = F, type = "b")
+axis(1, at = yrs, las = 2)
 
 tabp <- tab[grep("pelagic", splisref$PD), ]
 tabp <- tabp[order(rowSums(tabp), decreasing = T), ]
 
-matplot(1983:2020, t(tabp[1:20, ]), type = "l", col = glasbey(10), lwd = 2, lty = 1)
-legend("topright", rownames(tabp)[1:20], col = glasbey(10), lwd = 2, lty = 1)
+matplot(yrs, t(tabp[1:20, ]), type = "l", col = glasbey(10), lwd = 2, lty = 1)
+legend("topleft", rownames(tabp)[1:20], col = glasbey(10), lwd = 2, lty = 1, cex = 0.7)
 
 tabd <- tab[-grep("pelagic", splisref$PD), ]
 tabd <- tabd[order(rowSums(tabd), decreasing = T), ]
 
-matplot(1983:2020, t(tabd[1:20, ]), type = "b", col = glasbey(10), lwd = 2, lty = 1, pch = 19, las = 2)
+matplot(yrs, t(tabd[1:20, ]), type = "b", col = glasbey(10), lwd = 2, lty = 1, pch = 19, las = 2)
 legend("topright", rownames(tabd)[1:20], col = glasbey(10), lwd = 2, lty = 1)
 
 
diff --git a/indicator_processing/fishery_dependent/INDICATOR_fishery_indicators_STX.R b/indicator_processing/fishery_dependent/INDICATOR_fishery_indicators_STX.R
new file mode 100644
index 0000000..31fd06c
--- /dev/null
+++ b/indicator_processing/fishery_dependent/INDICATOR_fishery_indicators_STX.R
@@ -0,0 +1,276 @@
+# M. Karnauskas 10/19/2023
+# code for calculating pelagic:demersal ratio and Lmax indicators
+# uses logbook data for PR and USVI 
+
+rm(list = ls())
+library(pals)
+
+setwd("C:/Users/mandy.karnauskas/Desktop/Caribbean-ESR/indicator_processing/fishery_dependent/")
+
+dat <- read.csv("C:/Users/mandy.karnauskas/Desktop/CONFIDENTIAL/CaribbeanData/STX_072011_present_LANDINGS_trip_2021-03-11.csv")
+
+# define start and end years ---------------------------
+styear <- 2011
+enyear <- 2020
+table(dat$TRIP_YEAR)
+
+# subset years------------------------------------------
+
+d <- dat[which(dat$TRIP_YEAR >= styear & dat$TRIP_YEAR <= enyear), ]
+table(d$TRIP_YEAR)
+
+# take a look at data fields ----------------------------
+
+table(d$TRIP_ID, useNA = "always")
+table(d$TRIP_YEAR, useNA = "always")
+table(d$TRIP_MONTH, useNA = "always")
+table(d$TRIP_DAY, useNA = "always")
+table(d$VESSEL_CD, useNA = "always")
+table(d$FISHER_PERMIT, useNA = "always")
+table(d$FISHER_FIRST_NAME, useNA = "always")
+table(d$FISHER_LAST_NAME, useNA = "always")
+table(d$CHARTER_TRIP, useNA = "always")
+table(d$NUM_PARTNERS_OR_HELPERS, useNA = "always")
+table(d$IS_CATCH_SPLIT, useNA = "always")
+table(d$ISLAND, useNA = "always")
+table(d$TOTAL_TRAPS_IN_WATER, useNA = "always")
+table(d$FAD_CD, useNA = "always")
+table(d$GEAR_TYPE_NM, useNA = "always")
+table(d$GEAR1_NAME, useNA = "always")
+table(d$SPECIES_CD, useNA = "always")
+table(d$SPECIES_NM, useNA = "always")
+table(d$TRIP_MAY_BE_DUPLICATE, useNA = "always")
+
+# look at gear trends ---------------------------------
+tapply(d$POUNDS_LANDED, d$GEAR_TYPE_NM, sum, na.rm = T)
+par(mar = c(5, 15, 2, 2))
+barplot(tapply(d$POUNDS_LANDED, d$GEAR_TYPE_NM, sum, na.rm = T), horiz = T, las = 2)
+
+d$gear <- "other" 
+d$gear[which(d$GEAR_TYPE_NM == "TRAPS")] <- "traps"
+
+par(mar = c(4, 4, 1, 1))
+tab <- tapply(d$POUNDS_LANDED, list(d$gear, d$TRIP_YEAR), sum, na.rm = T)
+barplot(tab, col = c(2, gray(0.1), gray(0.5), gray(0.8)), args.legend = c(x = "topright"), legend.text = rownames(tab))
+
+# look at main species landed --------------------------------
+tab <- sort(tapply(d$POUNDS_LANDED, d$SPECIES_NM, sum, na.rm = T), decreasing = T)
+par(mar = c(15, 5, 2, 2))
+barplot(tab, las = 2)
+barplot(tab[1:50], las = 2)
+barplot(tab[1:25], las = 2)
+
+tab <- tapply(d$POUNDS_LANDED, list(d$SPECIES_NM, d$TRIP_YEAR), sum, na.rm = T)
+tab <- tab[order(rowSums(tab, na.rm = T), decreasing = T), ]
+
+par(mar = c(4, 4, 1, 1)+1)
+matplot(styear:enyear, t(tab[1:10, ]), type = "l", col = 1:10, lty = c(1, 1, 1, 1, 1, 2, 2, 2, 2, 2), lwd = 2)
+legend("topright", rownames(tab)[1:10], col = 1:10, lty = c(1, 1, 1, 1, 1, 2, 2, 2, 2, 2), lwd = 2)
+#abline(v = 1999)   
+
+#styear <- 2010
+#d <- d[which(d$TRIP_YEAR >= styear), ]  # cut out because very little data beforehand
+
+yrs <- styear:enyear
+
+tab <- tapply(d$POUNDS_LANDED, list(d$SPECIES_NM, d$TRIP_YEAR), sum, na.rm = T)
+tab <- tab[order(rowSums(tab, na.rm = T), decreasing = T), ]
+
+par(mar = c(4, 4, 1, 1)+1)
+matplot(yrs, t(tab[1:10, ]), type = "l", col = 1:10, lty = c(1, 1, 1, 1, 1, 2, 2, 2, 2, 2), lwd = 2)
+legend("topright", rownames(tab)[1:10], col = 1:10, lty = c(1, 1, 1, 1, 1, 2, 2, 2, 2, 2), lwd = 2)
+
+tab2 <- tab[1:10, ]
+for (i in 1:ncol(tab)) {  tab2[, i] <- tab[1:10, i] / sum(tab[1:10, i], na.rm = T)   }
+barplot(tab2, col = glasbey(10), xlim = c(0, 56), legend.text = rownames(tab2), args.legend = c(x = "right"))
+
+
+# pull in reference file and merge ------------------------------------
+
+ref <- read.csv("spp_ref_STT_manualEdit.csv")
+head(ref)
+head(d)
+
+hist(ref$Lmax)
+table(cut(ref$Lmax, breaks = c(0, 40, 60, 100, 200, 2000)))
+ref$Lmax_cat <- cut(ref$Lmax, breaks = c(0, 40, 60, 100, 200, 2000))
+
+db <- merge(d, ref, by.x = "SPECIES_NM", by.y = "COMname", all.x = TRUE)
+dim(d)
+dim(db)
+which(is.na(db$SPECIES_NM))
+
+tab <- tapply(db$POUNDS_LANDED, list(db$famcode, db$TRIP_YEAR), sum, na.rm = T)
+tab <- tab[order(rowSums(tab, na.rm = T), decreasing = T), ]
+head(tab, 15)
+
+nsp <- 15
+cols <- glasbey(nsp)
+
+par(mar = c(4, 4, 1, 1)+1)
+matplot(yrs, t(tab[1:nsp, ]), type = "l", col = cols, lty = rep(1:3, (nsp/3)), lwd = 2)
+legend("topright", rownames(tab)[1:nsp], col = cols, lty = rep(1:3, (nsp/3)), lwd = 2)
+
+tab2 <- tab[1:(nsp-1), ]
+tab2 <- rbind(tab2, colSums(tab[nsp:nrow(tab), ], na.rm = T))
+rownames(tab2)[nsp] <- "other"
+
+tab3 <- tab2
+for (i in 1:ncol(tab2)) {  tab3[, i] <- tab2[, i] / sum(tab2[, i], na.rm = T)   }
+
+cols[which(rownames(tab3) == "UNID")] <- "white"
+barplot(tab3, col = cols, xlim = c(0, 18), legend.text = rownames(tab2), args.legend = c(x = "right"), las = 2)
+
+#png(filename="PR_proportion_landings.png", units="in", width=8, height=5, pointsize=12, res=72*10)
+#barplot(tab3[, 23:38], col = cols, xlim = c(0, 23), legend.text = rownames(tab2), args.legend = c(x = "right", bty = "n"), las = 2, space = 0, border = NA, 
+#        xlab = "", ylab = "proportion of catch", main = "Puerto Rico landings composition by year                                ")
+#abline(h = 0)
+dev.off()
+
+#table(db$SPECIES_NM == db$COMname)
+table(db$SPECIES_NM[which(is.na(db$Lmax))])
+
+length(which(db$ITIS_COMMON_NAME == "FISHES,BONY,UNSPECIFIED"))
+table(db$PD, useNA = "always")
+table(db$Lmax, useNA = "always")
+
+# remove invertebrates and unidentified fish ----------------
+dbf <- db[which(db$PD != "invert"), ]
+dbf <- dbf[which(dbf$SPECIES_NM != "FISHES,BONY,UNSPECIFIED"), ]
+dbf <- dbf[which(dbf$SPECIES_NM != "OTHER SPECIES"), ]
+table(dbf$SPECIES_NM[which(is.na(dbf$Lmax))])
+
+head(dbf)
+dim(dbf)
+table(dbf$PD, useNA = "always")
+table(dbf$Lmax, useNA = "always")
+
+# calculate P:D ratio ---------------------------------------
+
+dbf$PD2 <- "demersal"
+dbf$PD2[grep("pelagic", dbf$PD)] <- "pelagic"
+table(dbf$PD, dbf$PD2)
+
+pd <- tapply(dbf$POUNDS_LANDED, list(dbf$TRIP_YEAR, dbf$PD2), sum, na.rm = T)
+pdrat <- pd[, 2] / pd[, 1]
+plot(yrs, pdrat, type = "b")
+
+# make indicator object and plot P:D ratio ------------------
+datdata <- yrs
+inddata <- data.frame(cbind(pd[, 2]/10^6, pd[, 1]/10^6, pdrat))
+labs <- c("Total pelagic landings", "millions of pounds" , "St. Croix", 
+          "Total demersal landings", "millions of pounds", "St. Croix", 
+          "Pelagic:demersal ratio", "ratio of landings", "St. Croix")
+indnames <- data.frame(matrix(labs, nrow = 3, byrow = F))
+
+s <- list(labels = indnames, indicators = inddata, datelist = datdata) #, ulim = ulidata, llim = llidata)
+class(s) <- "indicatordata"
+
+plotIndicatorTimeSeries(s, coltoplot = 1:3, plotrownum = 3, sublabel = T, trendAnalysis = F) #, outtype = "png")
+
+############################  END PD ratio  #######################################
+
+dev.off()
+
+# Lmax calculations -----------------------------------------
+
+tab <- tapply(dbf$POUNDS_LANDED, list(dbf$SPECIES_NM, dbf$TRIP_YEAR), sum, na.rm = T)
+tab[is.na(tab)] <- 0
+
+splis <- data.frame(rownames(tab))
+names(splis) <- "COMname"
+
+splisref <- merge(splis, ref, by = "COMname", all.x = TRUE)
+dim(splis)
+dim(splisref)
+table(rownames(tab) == splisref$COMname)
+
+llis <- levels(splisref$Lmax_cat)
+Lmax1 <- colSums(tab[which(splisref$Lmax_cat == llis[1]), ])
+Lmax2 <- colSums(tab[which(splisref$Lmax_cat == llis[2]), ])
+Lmax3 <- colSums(tab[which(splisref$Lmax_cat == llis[3]), ])
+Lmax4 <- colSums(tab[which(splisref$Lmax_cat == llis[4]), ])
+Lmax5 <- colSums(tab[which(splisref$Lmax_cat == llis[5]), ])
+#Lmax6 <- colSums(tab[which(splisref$Lmax_cat == llis[6]), ])
+
+Lmaxcl <- cbind(Lmax1, Lmax2, Lmax3, Lmax4, Lmax5) #, Lmax6)
+
+Lmaxcl2 <- Lmaxcl
+for (i in 1:nrow(Lmaxcl)) {  Lmaxcl2[i, ] <- Lmaxcl[i, ] / sum(Lmaxcl[i, ], na.rm = T)   }
+
+matplot(Lmaxcl)
+barplot(t(Lmaxcl2), col = 1:6, las = 2, main = "Distribution of catch in Lmax size classes                      ", 
+        ylab = "proportion", xlim = c(1, 30), 
+        legend = c("<40 cm", "40-60 cm", "60-100 cm", "100-200 cm", ">200 cm"), args.legend = c(x = "right", bty = "n"))
+
+# format indicator objects and plot ---------------------------------------
+datdata <- yrs
+inddata <- data.frame(Lmaxcl/10^6)
+labs <- c(rep("Total landings in Lmax class", 5), 
+          rep("millions of pounds", 5),
+          "<40 cm", "40-60 cm", "60-100 cm", "100-200 cm", ">200 cm")
+indnames <- data.frame(matrix(labs, nrow = 3, byrow = T))
+s <- list(labels = indnames, indicators = inddata, datelist = datdata) #, ulim = ulidata, llim = llidata)
+class(s) <- "indicatordata"
+plotIndicatorTimeSeries(s, coltoplot = 1:5, plotrownum = 5, sublabel = T, widadj = 1.5, trendAnalysis = F) #, outtype = "png")
+
+inddata <- data.frame(Lmaxcl2)
+labs <- c(rep("Proportion of landings in Lmax class",5), 
+          rep("proportion", 5),
+          "<40 cm", "40-60 cm", "60-100 cm", "100-200 cm", ">200 cm")
+indnames <- data.frame(matrix(labs, nrow = 3, byrow = T))
+s <- list(labels = indnames, indicators = inddata, datelist = datdata) #, ulim = ulidata, llim = llidata)
+class(s) <- "indicatordata"
+plotIndicatorTimeSeries(s, coltoplot = 1:5, plotrownum = 5, sublabel = T, widadj = 1.5, hgtadj = 0.6, trendAnalysis = F) # outtype = "png", sameYscale = F)
+
+# understand what is driving the trends -------------------------------------------
+
+splisref$COMname[which(splisref$Lmax_cat == "(0,40]")]
+splisref$COMname[which(splisref$Lmax_cat == "(100,200]")]
+
+sort(table(splisref$famcode[which(splisref$Lmax_cat == "(0,40]")]))
+sort(table(splisref$famcode[which(splisref$Lmax_cat == "(100,200]")]))
+
+splisref$recLand <- rowSums(tab)
+
+small <- splisref[which(splisref$Lmax_cat == "(0,40]"), ]
+head(small[order(small$recLand, decreasing = T), ], 15)
+
+big <- splisref[which(splisref$Lmax_cat == "(100,200]"), ]
+head(big[order(big$recLand, decreasing = T), ], 15)
+
+dev.off()
+
+# calculate average Lmax indicator -----------------------------------------------
+
+lmax <- rep(NA, ncol(tab))
+for (i in 1:ncol(tab))  {  lmax[i] <- weighted.mean(splisref$Lmax, tab[, i])  }
+
+plot(yrs, lmax)
+plot(yrs, lmax, type = "b", las = 2)
+axis(1, yrs)
+out <- lm(lmax ~ yrs)
+summary(out)
+abline(out)
+
+plot(yrs, tapply(dbf$Lmax, dbf$TRIP_YEAR, mean, na.rm = T), type = "l")
+
+# look at what is driving PD ratio  -------------------
+
+plot(yrs, pdrat, type = "b", las = 2)
+matplot(yrs, pd, axes = F, type = "b")
+axis(1, at = yrs, las = 2)
+
+tabp <- tab[grep("pelagic", splisref$PD), ]
+tabp <- tabp[order(rowSums(tabp), decreasing = T), ]
+
+matplot(yrs, t(tabp), type = "l", col = glasbey(10), lwd = 2, lty = 1)
+legend("topleft", rownames(tabp), col = glasbey(10), lwd = 2, lty = 1, cex = 0.7)
+
+tabd <- tab[-grep("pelagic", splisref$PD), ]
+tabd <- tabd[order(rowSums(tabd), decreasing = T), ]
+
+matplot(yrs, t(tabd[1:20, ]), type = "b", col = glasbey(10), lwd = 2, lty = 1, pch = 19, las = 2)
+legend("topright", rownames(tabd)[1:20], col = glasbey(10), lwd = 2, lty = 1)
+
+
diff --git a/indicator_processing/fishery_dependent/INDICATOR_total_landings.R b/indicator_processing/fishery_dependent/INDICATOR_total_landings.R
index 215f409..e17b429 100644
--- a/indicator_processing/fishery_dependent/INDICATOR_total_landings.R
+++ b/indicator_processing/fishery_dependent/INDICATOR_total_landings.R
@@ -1,8 +1,13 @@
+# M. Karnauskas 10/19/2023
 # code for calculating total landings or revenues
 # uses logbook data for PR and USVI 
 
 rm(list = ls())
 
+# define start and end years ---------------------------
+styear <- 1990
+enyear <- 2020
+
 # input data for Puerto Rico ---------------------------
 setwd("C:/Users/mandy.karnauskas/Desktop/Caribbean-ESR/indicator_processing/fishery_dependent/")
 dat <- read.csv("C:/Users/mandy.karnauskas/Desktop/CONFIDENTIAL/CaribbeanData/Jun2022/PR_landings_83_20_wSC_2005cor.csv")
@@ -13,18 +18,15 @@ table(dat$YEAR_LANDED)
 
 summary(dat$xADJ == dat$POUNDS_LANDED * 1/dat$CORRECTION_FACTOR)
 
-# define start and end years ---------------------------
-styear <- 1990
-enyear <- 2021
-
+# subset years------------------------------------------
 d <- dat[which(dat$YEAR_LANDED >= styear & dat$YEAR_LANDED <= enyear), ]
 
-# check field codes -----------------
+# check field codes ---------------------------------
 
 table(d$YEAR_LANDED, useNA = "always")
 table(d$CORRECTION_FACTOR, d$YEAR_LANDED)
 
-# look at top landings -----------------------
+# look at top landings -------------------------------
 
 sort(tapply(d$ADJUSTED_POUNDS, d$ITIS_COMMON_NAME, sum, na.rm = T))
 sort(table(d$ITIS_COMMON_NAME[grep("LOBSTER", d$ITIS_COMMON_NAME)]))
@@ -47,10 +49,12 @@ table(table(d$sppgrp, d$YEAR_LANDED) < 3)
 totland_pr <- tapply(d$ADJUSTED_POUNDS, list(d$YEAR_LANDED, d$sppgrp), sum, na.rm = T) / 10^3
 dim(totland_pr)
 totland_pr
-matplot(1990:2020, totland_pr, type = "l", lty = 1, lwd = 2)
+matplot(styear:enyear, totland_pr, type = "l", lty = 1, lwd = 2)
 
 rm(list = ls()[-match(c("totland_pr", "styear", "enyear"), ls())])
 
+#################     END PR    ########################
+
 # calculate for STT --------------------------------------
 
 dat <- read.csv("C:/Users/mandy.karnauskas/Desktop/CONFIDENTIAL/CaribbeanData/STT_landings.csv")
@@ -91,6 +95,8 @@ totland_st
 totland_st[is.na(totland_st)] <- 0
 matplot(totland_st, type = "l")
 
+#################     END STT    ########################
+
 # calculate for STX --------------------------------------
 
 rm(list = ls()[-match(c("totland_st", "totland_pr", "styear", "enyear"), ls())])
@@ -135,6 +141,9 @@ totland_sx
 #totland_sx[is.na(totland_sx)] <- 0
 matplot(totland_sx, type = "l")
 
+#################     END STX    ########################
+
+
 # summarize and plot ---------------------
 
 ls()[grep("totland", ls())]
@@ -162,7 +171,7 @@ setwd("C:/Users/mandy.karnauskas/Desktop/Caribbean-ESR/indicator_plots/")
 plotIndicatorTimeSeries(s, coltoplot = 1:9, plotrownum = 3, plotcolnum = 3, sublabel = T, sameYscale = F, 
                         widadj = 0.8, hgtadj = 0.7, trendAnalysis = F)   # outtype = "png")
 
-#inddata <- s
-#save(inddata, file = "C:/Users/mandy.karnauskas/Desktop/Caribbean-ESR/indicator_objects/landings.RData")
+inddata <- s
+save(inddata, file = "C:/Users/mandy.karnauskas/Desktop/Caribbean-ESR/indicator_objects/total_landings.RData")
 
 
diff --git a/indicator_processing/fishery_dependent/PREPROCESS_build_reference_file.R b/indicator_processing/fishery_dependent/PREPROCESS_build_reference_file.R
index 052df17..11f2f5c 100644
--- a/indicator_processing/fishery_dependent/PREPROCESS_build_reference_file.R
+++ b/indicator_processing/fishery_dependent/PREPROCESS_build_reference_file.R
@@ -261,14 +261,12 @@ table(d$SPECIES_ITIS, useNA = "always")
 table(d$ITIS_COMMON_NAME, useNA = "always")
 table(d$ITIS_SCIENTIFIC_NAME, useNA = "always")
 hist(d$POUNDS_LANDED)
+hist(d$ADJUSTED_POUNDS)
 hist(d$VALUE_IN_DOLLARS)
 hist(d$PRICE_PER_LB)
-table(d$CORRECTION_FACTOR, useNA = "always")
-table(d$CORRECTION_FACTOR, d$YEAR_LANDED)
-hist(d$ADJUSTED_POUNDS)
 table(d$DISTANCE, useNA = "always")
 table(d$DISTANCE_DESCRIBE, useNA = "always")
-table(d$TRIP_TICKET_NUMBER_ED, useNA = "always")
+#table(d$TRIP_TICKET_NUMBER_ED, useNA = "always")
 
 # compile list of species names into dataframe ----------------------------
 tab <- table(d$ITIS_COMMON_NAME, d$ITIS_SCIENTIFIC_NAME)