add STX data to disturbance indicator

indicator_plots/.Rhistory

@@ -0,0 +1,177 @@
+rm(list=ls())
+if (!"maps" %in% installed.packages()) install.packages("maps", repos='http://cran.us.r-project.org')
+library(maps)
+source("C:/Users/mandy.karnauskas/Desktop/RS_FATEproject/MASTER_codes/findMinDepth.R")
+load("C:/Users/mandy.karnauskas/Desktop/RS_FATEproject/MASTER_codes/GOMreleaseForScaling.RData")
+rm(list=ls())
+if (!"maps" %in% installed.packages()) install.packages("maps", repos='http://cran.us.r-project.org')
+library(maps)
+source("C:/Users/mandy.karnauskas/Desktop/completed_manuscripts/RS_FATEproject/MASTER_codes/findMinDepth.R")
+load("C:/Users/mandy.karnauskas/Desktop/completed_manuscripts/RS_FATEproject/MASTER_codes/GOMreleaseForScaling.RData")
+GOM <- matfinGOM
+load("C:/Users/mandy.karnauskas/Desktop/completed_manuscripts/RS_FATEproject/MASTER_codes/ATLreleaseForScaling.RData")
+ATL <- matS
+sum(GOM$V5)/10^5
+sum(ATL$V5)/10^5
+areaGOM <- length(unique(paste(GOM$V2, GOM$V3)))
+areaATL <- length(unique(paste(ATL$V2, ATL$V3)))
+areaGOM
+areaATL
+arearat <- areaGOM/areaATL
+arearat
+# calculate constants for scaling -------------------------------------
+const_scaler <- sum(GOM$V5)/sum(ATL$V5)
+const_scaler
+unscaled <- rbind(GOM, ATL)
+dim(unscaled)
+GOMsc <- GOM
+ATLsc <- ATL
+GOMsc$V5 <- round(GOMsc$V5 * 10)
+ATLsc$V5 <- round(ATLsc$V5 * 10 * const_scaler / arearat)
+scaled <- rbind(GOMsc, ATLsc)
+dim(scaled)
+which(scaled$V5==0)
+par(mfrow=c(1,2), mex=0.8)
+barplot(c(sum(GOM$V5),   sum(ATL$V5)),   names.arg=c("GoM", "Atl"), main="unscaled", ylab=""); abline(0,0)
+barplot(c(sum(GOMsc$V5), sum(ATLsc$V5)), names.arg=c("GoM", "Atl"), main="scaled", ylab=""); abline(0,0)
+# check ratios -----------------------------------------
+(sum(GOMsc$V5)) / (sum(ATLsc$V5))
+arearat
+(sum(GOM$V5)) / (sum(ATL$V5))
+dev.off()
+par(mar = c(5, 7, 1, 1))
+sc2 <- scaled[which(scaled$V6==2013 & scaled$V7==6 & scaled$V8==27),]
+x <- sc2$V5 / 130
+pos <- c(0.005,  0.05, 0.1, 0.2, 0.5, 1, 2, 4, 5, 10, 20, 50, 100, 1000)
+a <- floor(min(x))
+b <- max((x)-a)*1.03
+pind <- round((x-a)/b*100+1); print(min(pind)); print(max(pind))
+cols <- c(rainbow(30, start=0.82, end=0.99), rainbow(70, start=0.01, end=0.17))[100:1]
+map('state', fill = 1, interior=F, col = gray(0.85), ylim=c(22.5, 35), xlim=c(-88,-76))
+#mtext(side = 1, line = 2.5, "longitude")
+#mtext(side = 2, line = 2.5, "latitude")
+points(sc2$V2, sc2$V3, col=cols[pind], pch=15, cex=0.5)
+#box(); axis(1); axis(2, las = 2)
+xloc <- seq(-86, -78, length.out=100)
+for (j in 1:100) {   polygon(c(xloc[j], xloc[j+1],xloc[j+1], xloc[j]), c(23.0,23.0,23.4,23.4), col=cols[j], border=NA) }
+w <- which.min(abs(((max(x)-min(x))/6) - pos))
+if(-pos[w]<min(x)) { xx <- seq(0, max(x), pos[w]); xx <- xx[xx>min(x)] } else {  xx <- c(seq(-pos[w], min(x), -pos[w]), seq(0, max(x), pos[w])) }
+text(xloc[round((xx-a)/b*100+1)], y=22.75, xx, pos=2)
+text(-82, 24.25, "relative index of spawning output", pos = 1)
+degs = seq(88, 76, -2)
+a = sapply(degs, function(x) bquote(.(x)*degree ~ W))
+rm(list = ls())
+library(rfishbase)
+setwd("C:/Users/mandy.karnauskas/Desktop/CONFIDENTIAL/CaribbeanData/Jun2022/")
+rm(list = ls())
+setwd("C:/Users/mandy.karnauskas/Desktop/Caribbean-ESR/indicator_processing/fishery_dependent")
+sc <- read.csv("C:/Users/mandy.karnauskas/Desktop/CONFIDENTIAL/CaribbeanData/Jun2022/shellcatch_pr_data_req_02152023_C.csv")   # original shellcatch data
+sc$POUNDS_LANDED2 <- as.numeric(as.vector(sc$POUNDS_LANDED))
+sc$ADJUSTED_POUNDS2 <- as.numeric(as.vector(sc$ADJUSTED_POUNDS))
+sc$CORRECTION_FACTOR <- as.numeric(as.vector(sc$CORRECTION_FACTOR))
+sc$PRICE <- as.numeric(as.vector(sc$PRICE))
+summary(sc$POUNDS_LANDED == sc$POUNDS_LANDED2)
+summary(sc$ADJUSTED_POUNDS == sc$ADJUSTED_POUNDS2)
+apply(sc[1:4], 2, table, useNA = "always")
+table(sc$LANDING_AREA_COUNTY_OR_MUNICIPALITY)
+table(sc$COAST)
+table(sc$GEAR_NAME)
+hist(sc$POUNDS_LANDED2)
+hist(sc$ADJUSTED_POUNDS2)
+table(sc$CORRECTION_FACTOR)
+co <- sc$POUNDS_LANDED2 / sc$CORRECTION_FACTOR
+table(round(co- sc$ADJUSTED_POUNDS2))
+hist(sc$PRICE)
+table(sc$AREA_CD1)
+sc$SPECIES_NM <- as.character(sc$SPECIES_NM)
+# replace duplicate names -----------------------
+sc$SPECIES_NM[which(sc$SPECIES_NM == "CROAKERS")] <- "CROAKER"
+sc$SPECIES_NM[which(sc$SPECIES_NM == "HERRING")] <- "THREAD HERRING"
+sc$SPECIES_NM[which(sc$SPECIES_NM == "SMOOTHTAIL SPINY LOBSTER")] <- "SPINY LOBSTER"
+sc$SPECIES_NM[which(sc$SPECIES_NM == "KING MACKAREL, KINGFISH")] <- "KINGFISH MACKEREL"
+sc$SPECIES_NM[which(sc$SPECIES_NM == "WENCHMAN")] <- "CARDINAL"
+sc$SPECIES_NM[which(sc$SPECIES_NM == "SPOTTED TRUNKFISH")] <- "TRUNKFISH"
+sc$SPECIES_NM[which(sc$SPECIES_NM == "SMOOTH TRUNKFISH")] <- "TRUNKFISH"
+sc$SPECIES_NM[which(sc$SPECIES_NM == "SEA BASSES")] <- "GROUPERS"
+nam <- read.csv("name_matches_edited.csv")
+nam
+head(nam)
+dim(nam)
+unique(nam$shellcatch)
+length(unique(nam$shellcatch))
+length(unique(nam$logbook))
+length(unique(nam$logbook))
+unique(nam$second_name)
+tail(nam)
+unique(nam$third_name)
+length(unique(nam$shellcatch))
+unique(nam$shellcatch)
+unique(nam$logbook)
+135/300
+dat
+nam
+rm(list = ls())
+# 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.csv")
+table(dat$YEAR_LANDED)
+dat$xADJ <- dat$POUNDS_LANDED * 1/dat$CORRECTION_FACTOR
+hist(dat$ADJUSTED_POUNDS - dat$xADJ)
+max(abs(dat$ADJUSTED_POUNDS - dat$xADJ), na.rm = T)
+table(round(dat$ADJUSTED_POUNDS) - round(dat$xADJ))
+# define start and end years ---------------------------
+styear <- 1990
+enyear <- 2020
+d <- dat[which(dat$YEAR_LANDED >= styear & dat$YEAR_LANDED <= enyear), ]
+table(d$YEAR_LANDED, useNA = "always")
+table(d$CORRECTION_FACTOR, d$YEAR_LANDED)
+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)]))
+sort(table(d$ITIS_COMMON_NAME[grep("CONCH", d$ITIS_COMMON_NAME)]))
+d$sppgrp <- "other"
+#d$sppgrp[grep("LOBSTER", d$ITIS_COMMON_NAME)] <- "lobster"
+d$sppgrp[grep("LOBSTERS,SPINY", d$ITIS_COMMON_NAME)] <- "lobster"
+d$sppgrp[grep("CONCH", d$ITIS_COMMON_NAME)] <- "conch"
+table(d$ITIS_COMMON_NAME, d$sppgrp)
+totland_pr <- tapply(d$ADJUSTED_POUNDS, list(d$YEAR_LANDED, d$sppgrp), sum, na.rm = T) / 10^3
+dim(totland_pr)
+totland_pr
+matplot(totland_pr, type = "l", lty = 1, lwd = 2)
+ls()
+dat <- read.csv("C:/Users/mandy.karnauskas/Desktop/CONFIDENTIAL/CaribbeanData/STT_landings.csv")
+head(dat)
+table(dat$TRIP_YEAR)
+d <- dat[which(dat$TRIP_YEAR >= styear & dat$TRIP_YEAR <= enyear), ]
+table(d$TRIP_YEAR, useNA = "always")
+sort(tapply(d$POUNDS_LANDED, d$SPECIES_NM, sum, na.rm = T))
+sort(table(d$SPECIES_NM[grep("LOBSTER", d$SPECIES_NM)]))
+sort(table(d$SPECIES_NM[grep("CONCH", d$SPECIES_NM)]))
+d$sppgrp <- "other"
+d$sppgrp[grep("LOBSTER", d$SPECIES_NM)] <- "lobster"
+d$sppgrp[grep("CONCH", d$SPECIES_NM)] <- "conch"
+table(d$SPECIES_NM, d$sppgrp)
+totland_st <- tapply(d$POUNDS_LANDED, list(d$TRIP_YEAR, d$sppgrp), sum, na.rm = T) / 10^3
+dim(totland_st)
+totland_st
+totland_st[is.na(totland_st)] <- 0
+matplot(totland_st, type = "l")
+ls()[grep("totland", ls())]
+datdata <- styear:enyear
+inddata <- data.frame(cbind(totland_pr[, 2], totland_st[, 2],
+totland_pr[, 1], totland_st[, 1],
+totland_pr[, 3], totland_st[, 3]))
+labs <- c("Lobster landings", "thousands of pounds", "Puerto Rico",
+"Lobster landings", "thousands of pounds", "St. Thomas and St. John",
+"Conch landings", "thousands of pounds", "Puerto Rico",
+"Conch landings", "thousands of pounds", "St. Thomas and St. John",
+"Landings of all other species", "thousands of pounds", "Puerto Rico",
+"Landings of all other species",  "thousands of pounds", "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:6, plotrownum = 2, sublabel = T, sameYscale = F,
+widadj = 0.8, hgtadj = 0.7, trendAnalysis = F)   # outtype = "png")
+library(plotTimeSeries)
+plotIndicatorTimeSeries(s, coltoplot = 1:6, plotrownum = 2, sublabel = T, sameYscale = F,
+widadj = 0.8, hgtadj = 0.7, trendAnalysis = F)   # outtype = "png")