Skip to content

Commit

Permalink
update disturbance indicator
Browse files Browse the repository at this point in the history
  • Loading branch information
@MandyKarnauskas-NOAA
MandyKarnauskas-NOAA committed Jun 5, 2024
1 parent d0a926f commit adb9d1e
Show file tree
Hide file tree
Showing 5 changed files with 74 additions and 60 deletions.
4 changes: 2 additions & 2 deletions indicator_processing/CalculateAllIndicators.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -51,10 +51,10 @@ dir("indicator_processing/fishery_dependent/")
source("indicator_processing/fishery_dependent/INDICATOR_total_landings.R") # calc total landings by group
source("indicator_processing/fishery_dependent/INDICATOR_gearchanges.R") # produces plots but no indicator objects
source("indicator_processing/fishery_dependent/INDICATOR_gini.R") # calculate gini index based on revenues
source("indicator_processing/fishery_dependent/INDICATOR_disturbance.R") # calculate fishery disturbances based on distribution of landings over year
# not done below

#
source("indicator_processing/fishery_dependent/INDICATOR_disturbance.R") #

source("indicator_processing/fishery_dependent/INDICATOR_fishery_indicators_PR.R") #
source("indicator_processing/fishery_dependent/INDICATOR_fishery_indicators_STT.R") #
source("indicator_processing/fishery_dependent/INDICATOR_fishery_indicators_STX.R") #
Expand Down
116 changes: 64 additions & 52 deletions indicator_processing/fishery_dependent/INDICATOR_disturbance.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3,29 +3,30 @@
# uses logbook data for PR and USVI

# specification file and libraries -----------------------------

rm(list = ls())

library(maps)
library(plotTimeSeries)

load("indicator_processing/spec_file.RData")

confpath <- "C:/Users/mandy.karnauskas/Desktop/CONFIDENTIAL/CaribbeanData/"

# define years --------------------------------
confpath <- "C:/Users/mandy.karnauskas/Desktop/CONFIDENTIAL/CaribbeanData/MOST_RECENT/"

styear <- 1985
enyear <- 2020
# define start and end years ---------------------------
styear <- 2000
enyear <- 2022

# for PUERTO RICO -----------------------------------------------
# input data for Puerto Rico ---------------------------

dat <- read.csv(paste0(confpath, "/Jun2022/PR_landings_83_20_wSC_2005cor.csv"))
dat <- read.csv(paste0(confpath, "wrkeithly_pr_com_data_2000_2022_20240501_C.csv"))

d <- dat[which(dat$YEAR_LANDED >= styear & dat$YEAR_LANDED <= enyear), ]

# look at seasonal cycles ---------------------------------------

an <- tapply(d$ADJUSTED_POUNDS, list(d$MONTH_LANDED, d$YEAR_LANDED), sum, na.rm = T)
an
dim(an)
matplot(an, type = "l", col = rainbow(dim(an)[2]))

Expand Down Expand Up @@ -69,14 +70,14 @@ plot(styear:enyear, dst, type = "b")
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[1:20], las = 2)
lis <- names(tab)[1:9]
lis <- lis[-which(lis == "FISHES,BONY,UNSPECIFIED")]
lis <- names(tab)[1:10]
#lis <- lis[-which(lis == "FISHES,BONY,UNSPECIFIED")]

# recalculate index for individual spp --------------
dstfin <- c()
par(mfrow = c(2, 4), mex = 0.5)
par(mfrow = c(2, 5), mex = 0.5)

for (j in 1:8) {
for (j in 1:10) {
d1 <- d[which(d$ITIS_COMMON_NAME == lis[j]), ]
d1$YEAR_LANDED <- factor(d1$YEAR_LANDED, levels = styear:enyear)
an <- tapply(d1$ADJUSTED_POUNDS, list(d1$MONTH_LANDED, d1$YEAR_LANDED), sum, na.rm = T)
Expand All @@ -89,7 +90,7 @@ for (j in 1:8) {
dstfin <- cbind(dstfin, dst)
}

selec <- c(1, 4, 5, 6)
selec <- c(1, 3, 5, 7, 8, 9) # snappers, groupers subject to seasonal closure
lis[selec]

# look at correlations across spp --------------------
Expand All @@ -104,7 +105,7 @@ matplot(styear: enyear, dstfin, type = "b")
rownames(dstfin) <- styear: enyear
dstfin

avdst <- rowMeans(dstfin[, selec]) # snappers, groupers subject to seasonal closure
avdst <- rowMeans(dstfin[, selec])
sddst <- apply(dstfin[, selec], 1, sd)
plot(styear: enyear, avdst, type = "b")

Expand All @@ -122,30 +123,34 @@ names(final_PR) <- c("year", "ind", "lli", "uli")

rm(list = ls()[-match(c("final_PR", "styear", "enyear", "confpath"), ls())])

# load data --------------------------------
# calculate for STT --------------------------------------

dat <- read.csv(paste0(confpath, "STT_landings.csv"))
dat <- read.csv(paste0(confpath, "STT_2024.csv"))

table(dat$TRIP_YEAR, dat$TRIP_MONTH)

# adjust to fishing year ---------------------------------------
# adjust year to fishing year (Jul 1 - Jun 30) -------------

dat$TRIP_YEARold <- dat$TRIP_YEAR
dat$TRIP_MONTHold <- dat$TRIP_MONTH
aa <- which(dat$TRIP_MONTH < 7)
table(dat$TRIP_MONTH[aa])
table(dat$TRIP_YEAR[aa])
dat$TRIP_YEAR[aa] <- dat$TRIP_YEAR[aa] - 1
dat$TRIP_MONTH[aa] <- dat$TRIP_MONTH[aa] + 12
table(dat$TRIP_YEAR, dat$TRIP_MONTH)
tab <- table(dat$TRIP_YEAR, dat$TRIP_MONTH)

head(table(dat$TRIP_YEARold, dat$TRIP_MONTHold))
head(table(dat$TRIP_YEAR, dat$TRIP_MONTH))
tail(table(dat$TRIP_YEARold, dat$TRIP_MONTHold))
tail(table(dat$TRIP_YEAR, dat$TRIP_MONTH))
# take out incomplete years -----------------------------

d <- dat[which(dat$TRIP_YEAR >= styear & dat$TRIP_YEAR <= enyear), ]
table(dat$TRIP_YEAR)
lis <- as.numeric(names(which(apply(tab, 1, min) == 0)))
lis
dat <- dat[!(dat$TRIP_YEAR %in% lis), ]
table(dat$TRIP_YEAR)
table(dat$TRIP_YEAR, dat$TRIP_MONTH)

d <- dat[which(dat$TRIP_YEAR >= styear & dat$TRIP_YEAR <= enyear), ]
d$TRIP_YEAR <- factor(d$TRIP_YEAR, levels = c(styear: enyear))
table(d$TRIP_YEAR)


names(d)[which(names(d) == "TRIP_YEAR")] <- "YEAR_LANDED"
names(d)[which(names(d) == "TRIP_MONTH")] <- "MONTH_LANDED"
names(d)[which(names(d) == "SPECIES_NM")] <- "ITIS_COMMON_NAME"
Expand Down Expand Up @@ -176,10 +181,10 @@ lwds[which(colnames(an2) == 2020)] <- 3

matplot(an2, col = cols, type = "l", lty = 1, lwd = lwds, axes = F, ylab = "proportion of landings",
main = "Distribution of landings throughout the year")
axis(1, at = 1:12, lab = month.abb)
axis(1, at = 1:12, lab = month.abb[c(7:12, 1:6)])
axis(2, las = 2); box()
text(2, 0.11, col = 4, "2020", cex = 1.2, font = 2)
text(8, 0.025, col = 3, "2017", cex = 1.2, font = 2)
text(4, 0.025, col = 3, "2017", cex = 1.2, font = 2)
#lines(1:12, av, lwd = 3)
#polygon(c(1:12, 12:1), c((av + er), (av - er)[12:1]), col = "#00000020", border = NA)
#legend(10, 0.17, colnames(an)[which(pc %in% anomlis)], lty = 1, lwd = 2, col = 3:4, bty = "n")
Expand Down Expand Up @@ -209,14 +214,14 @@ table(db$famcode, useNA = "always")
tab <- sort(tapply(db$POUNDS_LANDED, db$famcode, sum, na.rm = T), decreasing = T)
par(mar = c(15, 5, 2, 2))
barplot(tab[1:20], las = 2)
lis <- names(tab)[1:9]
lis <- lis[-which(lis == "UNID")]
lis <- names(tab)[1:10]
#lis <- lis[-which(lis == "UNID")]

# recalculate index for individual spp --------------
dstfin <- c()
par(mfrow = c(2, 4), mex = 0.5)
par(mfrow = c(2, 5), mex = 0.5)

for (j in 1:8) {
for (j in 1:10) {
d1 <- db[which(db$famcode == lis[j]), ]
d1$YEAR_LANDED <- factor(d1$YEAR_LANDED, levels = styear:enyear)
an <- tapply(d1$POUNDS_LANDED, list(d1$MONTH_LANDED, d1$YEAR_LANDED), sum, na.rm = T)
Expand Down Expand Up @@ -245,6 +250,7 @@ rownames(dstfin) <- styear:enyear
dstfin # looks only reasonable starting in 2000

dstfin[which(styear: enyear < 2000), ] <- NA
dstfin[which(styear:enyear == 2022), ] <- NA
dstfin

avdst <- rowMeans(dstfin[, selec]) # snappers, groupers subject to seasonal closure
Expand All @@ -259,33 +265,38 @@ ulidata <- data.frame(avdst + sddst)
final_ST <- data.frame(inddata, llidata, ulidata)
names(final_ST) <- c("ind", "lli", "uli")

##################### END STT ######################

# for STX --------------------------------------------------------

rm(list = ls()[-match(c("final_PR", "final_ST", "styear", "enyear", "confpath"), ls())])

# load data --------------------------------

dat <- read.csv(paste0(confpath, "STX_072011_present_LANDINGS_trip_2021-03-11.csv"))
dat <- read.csv(paste0(confpath, "STX_2024.csv"))
head(dat)
table(dat$TRIP_YEAR, dat$TRIP_MONTH)

# adjust to fishing year ---------------------------------------
# adjust year to fishing year (Jul 1 - Jun 30) -------------

dat$TRIP_YEARold <- dat$TRIP_YEAR
dat$TRIP_MONTHold <- dat$TRIP_MONTH
aa <- which(dat$TRIP_MONTH < 7)
table(dat$TRIP_MONTH[aa])
table(dat$TRIP_YEAR[aa])
dat$TRIP_YEAR[aa] <- dat$TRIP_YEAR[aa] - 1
dat$TRIP_MONTH[aa] <- dat$TRIP_MONTH[aa] + 12
table(dat$TRIP_YEAR, dat$TRIP_MONTH)
tab <- table(dat$TRIP_YEAR, dat$TRIP_MONTH)

head(table(dat$TRIP_YEARold, dat$TRIP_MONTHold))
head(table(dat$TRIP_YEAR, dat$TRIP_MONTH))
tail(table(dat$TRIP_YEARold, dat$TRIP_MONTHold))
tail(table(dat$TRIP_YEAR, dat$TRIP_MONTH))
# take out incomplete years -----------------------------

d <- dat[which(dat$TRIP_YEAR >= styear & dat$TRIP_YEAR <= enyear), ]
table(dat$TRIP_YEAR)
lis <- as.numeric(names(which(apply(tab, 1, min) == 0)))
lis
dat <- dat[!(dat$TRIP_YEAR %in% lis), ]
table(dat$TRIP_YEAR)

# take a look at data fields ----------------------------

table(dat$LANDING_AREA_NAME)
table(dat$TRIP_YEAR)
d <- dat[which(dat$TRIP_YEAR >= styear & dat$TRIP_YEAR <= enyear), ]
d$TRIP_YEAR <- factor(d$TRIP_YEAR, levels = c(styear: enyear))
table(d$TRIP_YEAR)

names(d)[which(names(d) == "TRIP_YEAR")] <- "YEAR_LANDED"
Expand Down Expand Up @@ -318,7 +329,7 @@ lwds[which(colnames(an2) == 2020)] <- 3

matplot(an2, col = cols, type = "l", lty = 1, lwd = lwds, axes = F, ylab = "proportion of landings",
main = "Distribution of landings throughout the year")
axis(1, at = 1:12, lab = month.abb)
axis(1, at = 1:12, lab = month.abb[c(7:12, 1:6)])
axis(2, las = 2); box()
text(2, 0.11, col = 4, "2020", cex = 1.2, font = 2)
text(8, 0.025, col = 3, "2017", cex = 1.2, font = 2)
Expand All @@ -340,20 +351,21 @@ dim(d)
dim(db)
head(db)
which(is.na(db$famcode))
db$ITIS_COMMON_NAME[which(is.na(db$famcode))]
table(db$famcode, useNA = "always")

#tab <- sort(tapply(d$POUNDS_LANDED, d$ITIS_COMMON_NAME, sum, na.rm = T), decreasing = T)
tab <- sort(tapply(db$POUNDS_LANDED, db$famcode, sum, na.rm = T), decreasing = T)
par(mar = c(15, 5, 2, 2))
barplot(tab[1:20], las = 2)
lis <- names(tab)[1:8]
lis <- names(tab)[1:10]
#lis <- lis[-which(lis == "UNID")]

# recalculate index for individual spp --------------
dstfin <- c()
par(mfrow = c(2, 4), mex = 0.5)
par(mfrow = c(2, 5), mex = 0.5)

for (j in 1:8) {
for (j in 1:10) {
d1 <- db[which(db$famcode == lis[j]), ]
d1$YEAR_LANDED <- factor(d1$YEAR_LANDED, levels = styear:enyear)
an <- tapply(d1$POUNDS_LANDED, list(d1$MONTH_LANDED, d1$YEAR_LANDED), sum, na.rm = T)
Expand All @@ -366,7 +378,7 @@ for (j in 1:8) {
dstfin <- cbind(dstfin, dst)
}

selec <- c(1, 2, 4, 5)
selec <- c(2, 4, 5, 8)
lis[selec]

# look at correlations across spp --------------------
Expand All @@ -379,9 +391,10 @@ dev.off()

matplot(styear:enyear, dstfin, type = "b")
rownames(dstfin) <- styear:enyear
dstfin # only have data starting in 2012
dstfin # only have data starting in 2011

dstfin[which(styear:enyear < 2011), ] <- NA
dstfin[which(styear:enyear == 2022), ] <- NA

avdst <- rowMeans(dstfin[, selec]) #
sddst <- apply(dstfin[, selec], 1, sd)
Expand Down Expand Up @@ -418,7 +431,6 @@ plotIndicatorTimeSeries(s, sublabel = T, coltoplot = 1:3, plotrownum = 3, yposad
ind <- s
save(ind, file = "indicator_objects/disturbance.RData")



########################## END #############################


9 changes: 4 additions & 5 deletions indicator_processing/fishery_dependent/INDICATOR_fishery_indicators_PR.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -12,15 +12,14 @@ load("indicator_processing/spec_file.RData")

confpath <- "C:/Users/mandy.karnauskas/Desktop/CONFIDENTIAL/CaribbeanData/MOST_RECENT/"

# define start and end years ---------------------------
styear <- 2000
enyear <- 2022

# input data for Puerto Rico ---------------------------
setwd("C:/Users/mandy.karnauskas/Desktop/Caribbean-ESR/indicator_processing/fishery_dependent/")

dat <- read.csv(paste0(confpath, "Jun2022/PR_landings_83_20_wSC_2005cor.csv"))
dat <- read.csv(paste0(confpath, "wrkeithly_pr_com_data_2000_2022_20240501_C.csv"))

# define start and end years ---------------------------
styear <- 1990
enyear <- 2020
table(dat$YEAR_LANDED)

# subset years------------------------------------------
Expand Down
3 changes: 2 additions & 1 deletion indicator_processing/fishery_dependent/INDICATOR_gini.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -260,4 +260,5 @@ plotIndicatorTimeSeries(ind, coltoplot = 1:3, plotrownum = 3, sublabel = T, same
save(ind, file = "indicator_objects/gini.RData")



# hurricane Maria - 2017 - STX and PR
# hurricane Irma - 2017 - STT
2 changes: 2 additions & 0 deletions indicator_processing/fishery_dependent/PREPROCESS_build_reference_file.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -24,6 +24,8 @@ table(olddat$YEAR_LANDED)

# use only most recent pull b/c most indicators only are good back to 2020
# updated pull has higher sample numbers in recent years.
# it appears that issue with 2005 expansion factor has been resolved? - check
# does this new data set include shell catch?


# STT ---------------------------------
Expand Down

0 comments on commit adb9d1e

Please sign in to comment.