diff --git a/R/dat_proc.R b/R/dat_proc.R index 204b2b2..eaa47e4 100644 --- a/R/dat_proc.R +++ b/R/dat_proc.R @@ -220,8 +220,74 @@ dcgdat <- dcgraw %>% save(dcgdat, file = here('data/dcgdat.RData')) -# k brevis data ------------------------------------------------------------------------------- +# pyro data ----------------------------------------------------------------------------------- + +# https://f50006a.eos-intl.net/F50006A/OPAC/Details/Record.aspx?BibCode=5635517 +datall <- read.csv('https://f50006a.eos-intl.net/ELIBSQL12_F50006A_Documents/OTBMP_Pyrodinium_Chl_2011-2020_v101922.csv') %>% + select( + yr = Year, + date = Sample_Date, + Latitude, + Longitude, + pyro = P..bahamense.Abundance..cells.L. + ) %>% + mutate(date = mdy(date)) + +# 2021 only +dat2021 <- read.csv(url('https://raw.githubusercontent.com/tbep-tech/tbep-os-presentations/master/data/Pyrodinium_Chl_2021_OTBMP_mbeck.csv')) %>% + select( + date = Sample_Date, + Latitude, + Longitude, + pyro = Pbahamense..cells.L. + ) %>% + mutate( + date = case_when( + grepl('[a-z]', date) ~ dmy(date), + T ~ mdy(date) + ) + ) + +# 2022 only +dat2022 <- read.csv(url('https://raw.githubusercontent.com/tbep-tech/tbep-os-presentations/master/data/Pyrodinium_Chla_OTBMP_2022.csv')) %>% + select( + date = Date, + Latitude, + Longitude, + pyro = Pyrodinium..Cells.L. + ) %>% + mutate(date = mdy(date)) + +# 2023 only +tmpfile <- tempfile(fileext = '.xlsx') +download.file('https://github.com/tbep-tech/tbep-os-presentations/raw/master/data/2023%20OTB%20Pyrodinium%20bahamense%20abundance%20data.xlsx', + tmpfile, + mode = 'wb') +dat2023raw <- read_excel(tmpfile) +unlink(tmpfile) +dat2023 <- dat2023raw %>% + select( + date = `Sample Date`, + Latitude, + Longitude, + pyro = `Pyrodinium bahamense abundance (cells/L)` + ) %>% + mutate(date = ymd(date)) + +brks <- c(-Inf, 1e4, 1e5, 1e6, Inf) +labs <- c('No bloom', 'Low', 'Medium', 'High') + +pyrdat <- bind_rows(datall, dat2021, dat2022, dat2023) %>% + mutate( + yr = year(date), + doy = yday(date), + pyro = ifelse(pyro == 0, NA, pyro), + pyrocat = cut(pyro, breaks = brks, labels = labs) + ) +save(pyrdat, file = here::here('data/pyrdat.Rdata')) + +# k brevis data ------------------------------------------------------------------------------- # query api path <- 'https://gis.ncdc.noaa.gov/arcgis/rest/services/ms/HABSOS_CellCounts/MapServer/0/query?' diff --git a/R/hab_fib_scratch.R b/R/hab_fib_scratch.R deleted file mode 100644 index 3a8f7ee..0000000 --- a/R/hab_fib_scratch.R +++ /dev/null @@ -1,138 +0,0 @@ -library(tbeptools) -library(tidyverse) -library(readxl) -library(httr) -library(jsonlite) -library(sf) - -# pyro ---------------------------------------------------------------------------------------- - -# https://f50006a.eos-intl.net/F50006A/OPAC/Details/Record.aspx?BibCode=5635517 -datall <- read.csv('https://f50006a.eos-intl.net/ELIBSQL12_F50006A_Documents/OTBMP_Pyrodinium_Chl_2011-2020_v101922.csv') %>% - select( - yr = Year, - date = Sample_Date, - Latitude, - Longitude, - pyro = P..bahamense.Abundance..cells.L. - ) %>% - mutate(date = mdy(date)) - -# 2021 only -dat2021 <- read.csv(url('https://raw.githubusercontent.com/tbep-tech/tbep-os-presentations/master/data/Pyrodinium_Chl_2021_OTBMP_mbeck.csv')) %>% - select( - date = Sample_Date, - Latitude, - Longitude, - pyro = Pbahamense..cells.L. - ) %>% - mutate( - date = case_when( - grepl('[a-z]', date) ~ dmy(date), - T ~ mdy(date) - ) - ) - -# 2022 only -dat2022 <- read.csv(url('https://raw.githubusercontent.com/tbep-tech/tbep-os-presentations/master/data/Pyrodinium_Chla_OTBMP_2022.csv')) %>% - select( - date = Date, - Latitude, - Longitude, - pyro = Pyrodinium..Cells.L. - ) %>% - mutate(date = mdy(date)) - -# 2023 only -tmpfile <- tempfile(fileext = '.xlsx') -download.file('https://github.com/tbep-tech/tbep-os-presentations/raw/master/data/2023%20OTB%20Pyrodinium%20bahamense%20abundance%20data.xlsx', - tmpfile, - mode = 'wb') -dat2023raw <- read_excel(tmpfile) -unlink(tmpfile) -dat2023 <- dat2023raw %>% - select( - date = `Sample Date`, - Latitude, - Longitude, - pyro = `Pyrodinium bahamense abundance (cells/L)` - ) %>% - mutate(date = ymd(date)) - -brks <- c(-Inf, 1e4, 1e5, 1e6, Inf) -labs <- c('No bloom', 'Low', 'Medium', 'High') - -dat <- bind_rows(datall, dat2021, dat2022, dat2023) %>% - mutate( - yr = year(date), - doy = yday(date), - pyro = ifelse(pyro == 0, NA, pyro), - pyrocat = cut(pyro, breaks = brks, labels = labs), - pyro = pmin(3e6, pyro) - ) - -toplo <- dat %>% - filter(yr > 2011) - -toplomed <- toplo %>% - reframe( - medv = median(pyro, na.rm = T), - .by = yr - ) %>% - mutate( - pyrocat = cut(medv, breaks = brks, labels = labs), - ) - -p1 <- ggplot(toplo, aes(x = yr, y = pyro, group = yr)) + - geom_point(position = position_jitter(width = 0.2), alpha = 0.75, color = 'grey') + - scale_y_log10(breaks = c(1e3, 1e4, 1e5, 1e6), labels = parse(text = c('10^3', 'Low~(10^4)', 'Medium~(10^5)', 'High~(10^6)'))) + - scale_x_continuous(breaks = seq(min(toplo$yr), max(toplo$yr), 1)) + - geom_segment(data = toplomed, - aes(x = yr - 0.25, xend = yr + 0.25, - y = medv, yend = medv), - linewidth = 1 - ) + - theme_minimal() + - theme( - panel.grid.major.x = element_blank(), - panel.grid.minor = element_blank() - ) + - labs( - y = 'Bloom intensity (cells / L)', - x = NULL, - title = expression(paste(italic('Pyrodinium bahamense'), ' bloom intensity in Old Tampa Bay')), - subtitle = 'Observed cell counts and annual medians', - caption = 'Data source: Florida Fish and Wildlife Conservation Commission' - ) - - -# karenia ------------------------------------------------------------------------------------- - -load(file = here::here('data/kbrdat.RData')) - -# habdat p1 -toplo <- kbrdat %>% - st_set_geometry(NULL) %>% - filter(var == 'kb') %>% - filter(date < as.Date('2024-01-01') & date > as.Date('1960-01-01')) %>% - # filter(month(date) > 3 & month(date) < 10) %>% - mutate( - yr = year(date) - ) - -# plot -p2 <- ggplot(toplo, aes(x = yr, y = 1 + val)) + - geom_point(position = position_jitter(width = 0.1), alpha = 0.9) + - scale_y_log10(breaks = c(1e3, 1e4, 1e5, 1e6), labels = parse(text = c('10^3', 'Low~(10^4)', 'Medium~(10^5)', 'High~(10^6)')), limits= c(1000, NA)) + - # scale_x_discrete(breaks = seq(1950, 2025, by = 5)) + - theme_minimal() + - theme( - axis.text.x = element_text(angle = 45, size = 8, hjust = 1), - panel.grid.major.x = element_blank(), - panel.grid.minor = element_blank() - ) + - labs( - x = NULL - y = 'Cells / L (log-scale)' - ) - diff --git a/createfigs.R b/createfigs.R index c03164c..0ef5170 100644 --- a/createfigs.R +++ b/createfigs.R @@ -178,3 +178,92 @@ dev.off() png('figures/seagrasscov.png', height = 3.25, width = 6, res = 300, unit = 'in') sgcov_plo(seagrass, family = fml) dev.off() + +# habs ---------------------------------------------------------------------------------------- + +load(file = here::here('data/pyrdat.RData')) +load(file = here::here('data/kbrdat.RData')) + +# pyro +toplo1 <- pyrdat %>% + filter(yr > 2011) + +toplomed1 <- toplo1 %>% + reframe( + medv = median(pyro, na.rm = T), + .by = yr + ) + +p1 <- ggplot(toplo1, aes(x = yr, y = pyro, group = yr)) + + geom_point(position = position_jitter(width = 0.2), alpha = 0.75, color = 'grey') + + scale_y_log10(breaks = c(1e3, 1e4, 1e5, 1e6), labels = parse(text = c('10^3', 'Low~(10^4)', 'Medium~(10^5)', 'High~(10^6)')), + limits = c(1e3, NA)) + + scale_x_continuous(breaks = seq(min(toplo1$yr), max(toplo1$yr), 1)) + + geom_segment(data = toplomed1, + aes(x = yr - 0.25, xend = yr + 0.25, + y = medv, yend = medv), + linewidth = 1 + ) + + theme_minimal() + + theme( + panel.grid.major.x = element_blank(), + panel.grid.minor = element_blank(), + axis.text.x = element_text(angle = 45, size = 8, hjust = 1) + ) + + labs( + y = 'Bloom intensity (cells / L)', + x = NULL, + title = expression(paste(italic('Pyrodinium bahamense'), ' bloom intensity in Old Tampa Bay')), + subtitle = 'Observed cell counts and annual medians', + caption = 'Source: Florida Fish and Wildlife Conservation Commission' + ) + +# karenia +toplo2 <- kbrdat %>% + st_set_geometry(NULL) %>% + filter(var == 'kb') %>% + filter(date < as.Date('2024-01-01') & date > as.Date('1960-01-01')) %>% + # filter(month(date) > 3 & month(date) < 10) %>% + mutate( + yr = year(date), + val = ifelse(val == 0, NA, val) + ) %>% + filter(yr >= 1990) + +toplomed2 <- toplo2 %>% + reframe( + medv = median(val, na.rm = T), + .by = yr + ) + +# plot +p2 <- ggplot(toplo2, aes(x = yr, y = val)) + + geom_point(position = position_jitter(width = 0.2), alpha = 0.75, color = 'grey') + + scale_y_log10(breaks = c(1e3, 1e4, 1e5, 1e6), labels = parse(text = c('10^3', 'Low~(10^4)', 'Medium~(10^5)', 'High~(10^6)')), + limits= c(1e3, NA)) + + scale_x_continuous(breaks = seq(min(toplo2$yr), max(toplo2$yr), 1)) + + geom_segment(data = toplomed2, + aes(x = yr - 0.25, xend = yr + 0.25, + y = medv, yend = medv), + linewidth = 1 + ) + + theme_minimal() + + theme( + axis.text.x = element_text(angle = 45, size = 8, hjust = 1), + panel.grid.major.x = element_blank(), + panel.grid.minor = element_blank() + ) + + labs( + x = NULL, + y = 'Bloom intensity (cells / L)', + title = expression(paste(italic('Karenia brevis'), ' bloom intensity in Tampa Bay')), + subtitle = 'Observed cell counts and annual medians', + caption = 'Source: NOAA NCEI Harmful Algal BloomS Observing System (HABSOS)' + ) + +p <- p1 + p2 + plot_layout(ncol = 1) + +jpeg('figures/habs.jpg', family = fml, height = 5, width = 9, units = 'in', res = 300) +print(p) +dev.off() + diff --git a/data/pyrdat.Rdata b/data/pyrdat.Rdata new file mode 100644 index 0000000..1cae73d Binary files /dev/null and b/data/pyrdat.Rdata differ diff --git a/figures/habs.jpg b/figures/habs.jpg new file mode 100644 index 0000000..70645c3 Binary files /dev/null and b/figures/habs.jpg differ