Merge branch 'main' of https://github.com/Gulf-IEA/Caribbean-ESR

indicator_processing/automated_download/DHW.R

@@ -11,7 +11,7 @@ load("indicator_processing/spec_file.RData")
 
 options(download.file.method="libcurl")
 
-styear <- 2012
+styear <- 1980
 enyear <- terminal_year
 
 # download data and calculate mean -------------------

indicator_processing/automated_download/hotel_occupancy.R

@@ -7,8 +7,20 @@ rm(list = ls())
 library(pdftools)
 library(xml2)
 library(rvest)
+library(tidyverse)
+library(stringr)
+library(purrr)
 
 URL <- "https://usviber.org/archived-data/"
+
+page <- read_html(URL)
+
+raw_list <- page %>% # takes the page above for which we've read the html
+html_nodes("a") %>%  # find all links in the page
+html_attr("href") %>% # get the url for these links
+str_subset("\\.pdf") %>% # find those that end in pdf only
+
+
 pg <- read_html(URL)
 lis <- html_attr(html_nodes(pg, "a"), "href")
 lablis <- lis[9:32]

indicator_processing/fishery_dependent/INDICATOR_disturbance.R

@@ -19,7 +19,7 @@ enyear <- 2022
 
 # input data for Puerto Rico ---------------------------
 
-dat <- read.csv(paste0(confpath, "wrkeithly_pr_com_data_2000_2022_20240501_C.csv"))
+dat <- read.csv(paste0(confpath, "wrkeithly_pr_com_data_2000_2022_20240625_C.csv"))
 
 d <- dat[which(dat$YEAR_LANDED >= styear & dat$YEAR_LANDED <= enyear), ]
 
@@ -319,19 +319,19 @@ av <- apply(an2, 1, mean, na.rm = T)
 er <- apply(an2, 1, sd, na.rm = T)*1.96
 
 # identify anomalous years and highlight -----------------------
-anomlis <- c(2017, 2020)
+anomlis <- c(2017, 2019)
 cols <- rep("#00000060", dim(an)[2])
 cols[which(colnames(an2) == 2017)] <- 3
-cols[which(colnames(an2) == 2020)] <- 4
+cols[which(colnames(an2) == 2019)] <- 4
 lwds <- rep(1, dim(an)[2])
 lwds[which(colnames(an2) == 2017)] <- 3
-lwds[which(colnames(an2) == 2020)] <- 3
+lwds[which(colnames(an2) == 2019)] <- 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[c(7:12, 1:6)])
 axis(2, las = 2); box()
-text(2, 0.11, col = 4, "2020", cex = 1.2, font = 2)
+text(2, 0.11, col = 4, "2019", cex = 1.2, font = 2)
 text(8, 0.025, col = 3, "2017", cex = 1.2, font = 2)
 
 # calculate disturbance indicator ---------------------------

indicator_processing/fishery_dependent/INDICATOR_fishery_indicators_PR.R

@@ -20,7 +20,7 @@ enyear <- 2022
 
 # input data for Puerto Rico ---------------------------
 
-dat <- read.csv(paste0(confpath, "wrkeithly_pr_com_data_2000_2022_20240501_C.csv"))
+dat <- read.csv(paste0(confpath, "wrkeithly_pr_com_data_2000_2022_20240625_C.csv"))
 
 table(dat$YEAR_LANDED)
 

indicator_processing/fishery_dependent/INDICATOR_gearchanges.R

@@ -547,7 +547,7 @@ dev.off()
 ###########################  START PR  #################################
 rm(list = ls()[-match(c("styear", "enyear", "confpath"), ls())])
 
-dat <- read.csv(paste0(confpath, "wrkeithly_pr_com_data_2000_2022_20240501_C.csv"))
+dat <- read.csv(paste0(confpath, "wrkeithly_pr_com_data_2000_2022_20240625_C.csv"))
 
 # subset years------------------------------------------
 
@@ -742,6 +742,7 @@ save(byc, file ="indicator_data/prop_trips_nonselective_PR.RData")
 
 # plot trips and percent bycatch gears -------------------
 
+mat
 mat <- mat[, c(1, 3, 4, 2, 5)]
 cols <- cols25(7)
 cols <- cols[c(7, 2:5)]

indicator_processing/fishery_dependent/INDICATOR_gini.R

@@ -27,7 +27,7 @@ enyear <- 2022
 
 # input data for Puerto Rico ---------------------------
 
-dat <- read.csv(paste0(confpath, "wrkeithly_pr_com_data_2000_2022_20240501_C.csv"))
+dat <- read.csv(paste0(confpath, "wrkeithly_pr_com_data_2000_2022_20240625_C.csv"))
 
 d <- dat[which(dat$YEAR_LANDED >= styear & dat$YEAR_LANDED <= enyear), ]
 
@@ -37,6 +37,22 @@ table(d$YEAR_LANDED, useNA = "always")
 table(d$PR_ID_CODE_ED, useNA = "always")
 table(d$PR_ID_CODE_ED, d$YEAR_LANDED)
 
+# remove land crab trips -------------------
+
+sort(table(d$ITIS_COMMON_NAME[grep("CRAB", d$ITIS_COMMON_NAME)]))
+sort(table(d$ERDMAN_GEAR_NAME[grep("CRAB,BLUE", d$ITIS_COMMON_NAME)]))
+sort(table(d$FAMILY[grep("CRAB,BLUE", d$ITIS_COMMON_NAME)]))
+sort(table(d$ERDMAN_GEAR_NAME[which(d$FAMILY == "LAND CRABS")]))
+par(mar = c(12, 4, 1, 1))
+barplot(sort(table(d$ERDMAN_GEAR_NAME[which(d$FAMILY == "LAND CRABS")])), las = 2)
+sort(table(d$ITIS_COMMON_NAME[which(d$ERDMAN_GEAR_NAME == "BY HAND")]))  # don't remove, contains conch as well
+# filtering by LAND CRAB TRAP gear takes out vast majority (96%) of land crab trips
+
+d[which(d$ERDMAN_GEAR_NAME == "LAND CRAB TRAP"), ]
+dim(d)
+d <- d[which(d$ERDMAN_GEAR_NAME != "LAND CRAB TRAP"), ]
+dim(d)
+
 # remove bad price values ------------------------------
 #Notes from Juan: 
 # individual catch/revenues assignments to a single vessel are likely to be better after 2012 when reporting improved 
@@ -54,9 +70,21 @@ d$PRICE_PER_LB[(which(d$ITIS_COMMON_NAME != "CRAB,BLUE LAND" & d$PRICE_PER_LB >
 hist(d$PRICE_PER_LB[(which(d$ITIS_COMMON_NAME == "CRAB,BLUE LAND"))])
 max(d$PRICE_PER_LB[(which(d$ITIS_COMMON_NAME == "CRAB,BLUE LAND"))])
 
+# insert missing prices -------------------------------
+
+table(is.na(d$PRICE_PER_LB), d$ITIS_COMMON_NAME)
+table(is.na(d$PRICE_PER_LB))
+
+for (i in unique(d$ITIS_COMMON_NAME))  {
+  m <- mean(d$PRICE_PER_LB[which(d$ITIS_COMMON_NAME == i)], na.rm = T)
+  d$PRICE_PER_LB[which(d$ITIS_COMMON_NAME == i & is.na(d$PRICE_PER_LB))] <- m
+}
+
+table(is.na(d$PRICE_PER_LB), d$ITIS_COMMON_NAME)
+table(is.na(d$PRICE_PER_LB))
+
 # calculate revenue, sum by permit and year ------------
 
-#d$REV <- d$POUNDS_LANDED * d$PRICE_PER_LB
 d$REV <- d$ADJUSTED_POUNDS * d$PRICE_PER_LB
 totrev <- tapply(d$REV, list(d$PR_ID_CODE_ED, d$YEAR_LANDED), sum, na.rm = T)
 dim(totrev)
@@ -118,6 +146,13 @@ 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$TRIP_MAY_BE_DUPLICATE, useNA = "always")
+
+dim(d)
+d <- d[-which(d$TRIP_MAY_BE_DUPLICATE == "Y"), ]
+dim(d)
+d <- d[-which(d$CHARTER_TRIP == "Y"), ]
+dim(d)
 
 # remove bad price values ------------------------------
 
@@ -127,6 +162,28 @@ d$SPECIES_NM[which(d$PRICE > 15)]
 hist(d$PRICE[which(d$PRICE > 15)])
 d$PRICE[which(d$PRICE > 15)] <- NA
 
+# insert missing prices -------------------------------
+
+table(is.na(d$PRICE), d$SPECIES_NM)
+table(is.na(d$PRICE))
+
+for (i in unique(d$SPECIES_NM))  {
+  m <- mean(d$PRICE[which(d$SPECIES_NM == i)], na.rm = T)
+  d$PRICE[which(d$SPECIES_NM == i & is.na(d$PRICE))] <- m
+}
+
+table(is.na(d$PRICE))
+table(d$SPECIES_NM, is.na(d$PRICE))
+table(d$FAMILY, is.na(d$PRICE))
+
+for (i in unique(d$FAMILY))  {
+  m <- mean(d$PRICE[which(d$FAMILY == i)], na.rm = T)
+  d$PRICE[which(d$FAMILY == i & is.na(d$PRICE))] <- m
+}
+
+table(is.na(d$PRICE))
+table(d$FAMILY, is.na(d$PRICE))
+
 # calculate revenue, sum by permit and year ------------
 
 d$REV <- d$POUNDS_LANDED * d$PRICE
@@ -152,11 +209,12 @@ dim(totland)
 
 # calculate gini index --------------------------------
 
+par(mfrow = c(2, 1))
 gini_rev_stt <- apply(totrev, 2, calcGini)
 plot(names(gini_rev_stt), gini_rev_stt, type = "b")
 
 gini_land_stt <- apply(totland, 2, calcGini)
-plot(names(gini_land_stt), gini_land_stt, col = 2)
+plot(names(gini_land_stt), gini_land_stt, col = 2, type = "b")
 
 plot(gini_land_stt, gini_rev_stt)
 cor(gini_land_stt, gini_rev_stt)
@@ -192,8 +250,35 @@ 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$TRIP_MAY_BE_DUPLICATE, useNA = "always")
+
+dim(d)
+#d <- d[-which(d$TRIP_MAY_BE_DUPLICATE == "Y"), ]
+#dim(d)
+d <- d[-which(d$CHARTER_TRIP == "Y"), ]
+dim(d)
+
+# insert missing prices -------------------------------
 
-#d <- d[which(d$CHARTER_TRIP == "Y"), ]
+table(is.na(d$PRICE), d$SPECIES_NM)
+table(is.na(d$PRICE))
+
+for (i in unique(d$SPECIES_NM))  {
+  m <- mean(d$PRICE[which(d$SPECIES_NM == i)], na.rm = T)
+  d$PRICE[which(d$SPECIES_NM == i & is.na(d$PRICE))] <- m
+}
+
+table(is.na(d$PRICE))
+table(d$SPECIES_NM, is.na(d$PRICE))
+table(d$FAMILY, is.na(d$PRICE))
+
+for (i in unique(d$FAMILY))  {
+  m <- mean(d$PRICE[which(d$FAMILY == i)], na.rm = T)
+  d$PRICE[which(d$FAMILY == i & is.na(d$PRICE))] <- m
+}
+
+table(is.na(d$PRICE))
+table(d$FAMILY, is.na(d$PRICE))
 
 # remove bad price values ------------------------------
 
@@ -244,7 +329,7 @@ cor(gini_land_stx, gini_rev_stx)
 ls()[grep("gini", ls())]
 
 datdata <- styear:enyear
-inddata <- data.frame(gini_rev_pr, gini_rev_stt, gini_land_stx)
+inddata <- data.frame(gini_rev_pr, gini_rev_stt, gini_rev_stx)
 labs <- c("Inequality in revenues" , "Gini index", "Puerto Rico", 
           "Inequality in revenues" , "Gini index", "St. Thomas and St. John",
           "Inequality in revenues" , "Gini index", "St. Croix")
@@ -255,9 +340,28 @@ class(s) <- "indicatordata"
 ind <- s 
 
 plotIndicatorTimeSeries(ind, coltoplot = 1:3, plotrownum = 3, sublabel = T, sameYscale = T, 
-                        widadj = 1.3, hgtadj = 1, trendAnalysis = T)
+                        widadj = 1.5, hgtadj = 0.8, trendAnalysis = T, outtype = "png")
+
+save(ind, file = "indicator_objects/gini_revenue.RData")
+
+
+ls()[grep("gini", ls())]
+
+datdata <- styear:enyear
+inddata <- data.frame(gini_land_pr, gini_land_stt, gini_land_stx)
+labs <- c("Inequality in landings" , "Gini index", "Puerto Rico", 
+          "Inequality in landings" , "Gini index", "St. Thomas and St. John",
+          "Inequality in landings" , "Gini index", "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"
+
+ind <- s 
+
+plotIndicatorTimeSeries(ind, coltoplot = 1:3, plotrownum = 3, sublabel = T, sameYscale = T, 
+                        widadj = 1.3, hgtadj = 0.8, trendAnalysis = T)
 
-save(ind, file = "indicator_objects/gini.RData")
+save(ind, file = "indicator_objects/gini_landings.RData")
 
 
 # hurricane Maria - 2017 - STX and PR

indicator_processing/fishery_dependent/INDICATOR_total_landings.R

@@ -18,7 +18,7 @@ enyear <- 2022
 
 # input data for Puerto Rico ---------------------------
 
-dat <- read.csv(paste0(confpath, "wrkeithly_pr_com_data_2000_2022_20240501_C.csv"))
+dat <- read.csv(paste0(confpath, "wrkeithly_pr_com_data_2000_2022_20240625_C.csv"))
 
 table(dat$YEAR_LANDED)
 
@@ -37,6 +37,22 @@ d <- dat[which(dat$YEAR_LANDED >= styear & dat$YEAR_LANDED <= enyear), ]
 table(d$YEAR_LANDED, useNA = "always")
 table(d$CORRECTION_FACTOR, d$YEAR_LANDED)
 
+# remove land crab trips -------------------
+
+sort(table(d$ITIS_COMMON_NAME[grep("CRAB", d$ITIS_COMMON_NAME)]))
+sort(table(d$ERDMAN_GEAR_NAME[grep("CRAB,BLUE", d$ITIS_COMMON_NAME)]))
+sort(table(d$FAMILY[grep("CRAB,BLUE", d$ITIS_COMMON_NAME)]))
+sort(table(d$ERDMAN_GEAR_NAME[which(d$FAMILY == "LAND CRABS")]))
+par(mar = c(12, 4, 1, 1))
+barplot(sort(table(d$ERDMAN_GEAR_NAME[which(d$FAMILY == "LAND CRABS")])), las = 2)
+sort(table(d$ITIS_COMMON_NAME[which(d$ERDMAN_GEAR_NAME == "BY HAND")]))  # don't remove, contains conch as well
+# filtering by LAND CRAB TRAP gear takes out vast majority (96%) of land crab trips
+
+d[which(d$ERDMAN_GEAR_NAME == "LAND CRAB TRAP"), ]
+dim(d)
+d <- d[which(d$ERDMAN_GEAR_NAME != "LAND CRAB TRAP"), ]
+dim(d)
+
 # look at top landings -------------------------------
 
 sort(tapply(d$ADJUSTED_POUNDS, d$ITIS_COMMON_NAME, sum, na.rm = T))
@@ -222,7 +238,7 @@ class(s) <- "indicatordata"
 ind <- s
 
 plotIndicatorTimeSeries(ind, coltoplot = 1:9, plotrownum = 3, plotcolnum = 3, sublabel = T, sameYscale = F, 
-                        widadj = 0.8, hgtadj = 0.7, trendAnalysis = F)   # outtype = "png")
+                        widadj = 0.8, hgtadj = 0.7, trendAnalysis = T) #, outtype = "png")
 
 save(ind, file = "indicator_objects/total_landings.RData")
 

indicator_processing/non_automated/Sargassum_inundation.R

@@ -37,8 +37,11 @@ res <- griddap(afai, fields = 'AFAI',
                     longitude = c(min_lon, max_lon), 
                     latitude = c(min_lat, max_lat))
 
+hist(res$data$AFAI)
+quantile(res$data$AFAI, na.rm = T, probs = 0.95)
+
 res$data$lim <- 0
-res$data$lim [which(res$data$AFAI > 0.002)] <- 1
+res$data$lim[which(res$data$AFAI > 0.00063)] <- 1
 
 res$data$year <- substr(res$data$time, 1, 4)
 res$data$tim <- substr(res$data$time, 1, 7)
@@ -49,8 +52,7 @@ mon_afai <- tapply(res$data$lim, res$data$tim, sum, na.rm = T)
 yr_afai <- tapply(res$data$lim, res$data$year, sum, na.rm = T)
 
 plot(mon_afai, type = "l")
-plot(2017:2023, 
-     yr_afai, type = "l")
+plot(2017:2023, yr_afai, type = "l")
 
 # create indicator object --------------