temperature_july.R

`%>%` <- magrittr::`%>%`

# shapefile setup
crs <- "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"
# equal area crs
new_crs <- "+proj=utm +zone=12 +datum=NAD83 +no_defs +ellps=GRS80"

# different geom fix
sf::sf_use_s2(FALSE)

mab <- NEesp::shape %>%
  dplyr::select(STRATA, geometry) %>%
  sf::st_transform(proj4string = new_crs) %>% 
  # try geom fix?
  # dplyr::mutate(geometry = geometry %>% 
  #                 s2::s2_rebuild() %>%
  #                 sf::st_as_sfc()) %>%
  dplyr::summarise(geometry = sf::st_union(geometry)) %>%  
  sf::st_crop(y = c(xmin = -80, xmax = -69, 
                    ymax = 41.5, ymin =  35.8327))

strata <- readxl::read_excel(here::here("data-raw/BLF_STRATA.xlsx"),
                             skip = 3)
colnames(strata) <- c("alb", "big")
strata <- c(strata$alb, strata$big) %>%
  unique() 
print(strata)

bf_strata <- NEesp::shape %>%
  dplyr::select(STRATA, geometry) %>%
  dplyr::filter(STRATA %in% strata) %>%
  sf::st_transform(proj4string = new_crs) %>% 
  dplyr::summarise(geometry = sf::st_union(geometry))

years <- 1982:2021
prop_july <- c()
prop_july_bf <- c()

prop_july_cold <- c()
prop_july_bf_cold <- c()

prop_july_warm <- c()
prop_july_bf_warm <- c()

for(j in years) {
  message(paste("starting", j))
  # download data ----
  dir.create(here::here("data-raw","gridded", "sst_data"), recursive = TRUE)
  url <- paste0("https://downloads.psl.noaa.gov/Datasets/noaa.oisst.v2.highres/sst.day.mean.", j, ".v2.nc")
  download.file(url, destfile = "test.nc")
  
  # R can't open the file (will have to do this in a gh action...)
  # download file manually for testing on desktop
  
   # name <- paste0(j, ".nc")
  name <- "test.nc"
  
  data <- ecopull::nc_to_raster(nc = name, varname = 'sst') # converts to NAD83
  data <- raster::rotate(data)
  message("converted to raster...")
  
  # make sure all days are there ----
  
  if(raster::nlayers(data) < 365) {
    message(j, " does not have a full year of data! skipping!")
  } else {
  
    # filter to just july
    july <- c(paste0("X", j, ".07.0", 1:9),
              paste0("X", j, ".07.", 10:31))
    
  # crop to MAB ----
  # ndays <- raster::nlayers(data) # account for leap years
  
  mab_temp <- raster::mask(x = data[[july]], 
                           mask = mab
                           )
  message("cropped to July MAB...")

  # calculate total area ----
  total_area <- raster::area(mab_temp, na.rm = TRUE) %>%
    raster::as.data.frame(xy = TRUE) %>%
    dplyr::select(-c(x, y)) %>%
    colSums(na.rm = TRUE)
  
  # calculate area 18-25.6C ----
  mab_temp@data@values[which(mab_temp@data@values < 18)] <- NA
  # also remove areas that are too warm
  mab_temp@data@values[which(mab_temp@data@values > 25.6)] <- NA
  
  july_area <- raster::area(mab_temp, na.rm = TRUE) %>%
    raster::as.data.frame(xy = TRUE) %>%
    dplyr::select(-c(x, y)) %>%
    colSums(na.rm = TRUE)
  
  # * calculate proportion ----
  
  props <- july_area/total_area
  this_prop <- mean(props)
  
  prop_july <- c(prop_july, this_prop)
  message("calculated MAB proportion...")
  
  # calculate area under 18C ----
  mab_temp <- raster::mask(x = data[[july]], 
                           mask = mab
  )
  mab_temp@data@values[which(mab_temp@data@values >= 18)] <- NA
  
  july_area <- raster::area(mab_temp, na.rm = TRUE) %>%
    raster::as.data.frame(xy = TRUE) %>%
    dplyr::select(-c(x, y)) %>%
    colSums(na.rm = TRUE)
  
  # * calculate proportion ----
  
  props <- july_area/total_area
  this_prop <- mean(props)
  
  prop_july_cold <- c(prop_july_cold, this_prop)
  message("calculated MAB proportion...")
  
  # calculate area over 25.6C ----
  mab_temp <- raster::mask(x = data[[july]], 
                           mask = mab
  )
  mab_temp@data@values[which(mab_temp@data@values <= 25.6)] <- NA
  
  july_area <- raster::area(mab_temp, na.rm = TRUE) %>%
    raster::as.data.frame(xy = TRUE) %>%
    dplyr::select(-c(x, y)) %>%
    colSums(na.rm = TRUE)
  
  # * calculate proportion ----
  
  props <- july_area/total_area
  this_prop <- mean(props)
  
  prop_july_warm <- c(prop_july_warm, this_prop)
  message("calculated MAB proportion...")
  
  # crop to bluefish strata ----
  
  mab_temp <- raster::mask(x = data[[july]], 
                           mask = bf_strata
  )
  message("cropped to July MAB...")
  
  # calculate total area ----
  total_area <- raster::area(mab_temp, na.rm = TRUE) %>%
    raster::as.data.frame(xy = TRUE) %>%
    dplyr::select(-c(x, y)) %>%
    colSums(na.rm = TRUE)
  
  # calculate area 18-25.6C ----
  mab_temp@data@values[which(mab_temp@data@values < 18)] <- NA
  # also remove areas that are too warm
  mab_temp@data@values[which(mab_temp@data@values > 25.6)] <- NA
  
  july_area <- raster::area(mab_temp, na.rm = TRUE) %>%
    raster::as.data.frame(xy = TRUE) %>%
    dplyr::select(-c(x, y)) %>%
    colSums(na.rm = TRUE)
  
  # * calculate proportion ----
  
  props <- july_area/total_area
  this_prop <- mean(props)
  
  prop_july_bf <- c(prop_july_bf, this_prop)
  message("calculated bluefish strata proportion...")
  
  # calculate area under 18C ----
  mab_temp <- raster::mask(x = data[[july]], 
                           mask = bf_strata
  )
  mab_temp@data@values[which(mab_temp@data@values >= 18)] <- NA
  
  july_area <- raster::area(mab_temp, na.rm = TRUE) %>%
    raster::as.data.frame(xy = TRUE) %>%
    dplyr::select(-c(x, y)) %>%
    colSums(na.rm = TRUE)
  
  # * calculate proportion ----
  
  props <- july_area/total_area
  this_prop <- mean(props)
  
  prop_july_bf_cold <- c(prop_july_bf_cold, this_prop)
  message("calculated bluefish strata proportion...")
  
  # calculate area over 25.6C ----
  mab_temp <- raster::mask(x = data[[july]], 
                           mask = bf_strata
  )
  mab_temp@data@values[which(mab_temp@data@values <= 25.6)] <- NA
  
  july_area <- raster::area(mab_temp, na.rm = TRUE) %>%
    raster::as.data.frame(xy = TRUE) %>%
    dplyr::select(-c(x, y)) %>%
    colSums(na.rm = TRUE)
  
  # * calculate proportion ----
  
  props <- july_area/total_area
  this_prop <- mean(props)
  
  prop_july_bf_warm <- c(prop_july_bf_warm, this_prop)
  message("calculated bluefish strata proportion...")
  
  }
  message(paste("done with", j))
}

# prop_july

out_data <- tibble::tibble(years, prop_july, prop_july_bf, 
                           prop_july_cold, prop_july_bf_cold, 
                           prop_july_warm, prop_july_bf_warm)
write.csv(out_data, here::here("data-raw/temperature_july_bfstrata.csv"))