From c0c0d965314b63914eb79fee3c7caeb4404d1536 Mon Sep 17 00:00:00 2001 From: SarahJWeisberg Date: Wed, 8 Jan 2025 11:41:45 -0500 Subject: [PATCH] EPU-wide avg of daily SST - test --- R/avg_sst_EPU.R | 36 ++++++++++++++++++++++++++++++++++++ 1 file changed, 36 insertions(+) create mode 100644 R/avg_sst_EPU.R diff --git a/R/avg_sst_EPU.R b/R/avg_sst_EPU.R new file mode 100644 index 0000000..b118a05 --- /dev/null +++ b/R/avg_sst_EPU.R @@ -0,0 +1,36 @@ +##SW version +#goal is to get EPU-wide estimates of average daily sst (from OISST) +#used kevins methods to go from global netcdf to global raster then local (NE shelf) raster + +#load packages +library(terra) +library(sf) +#load raster file +load("~/Desktop/CurrentConditions/oisst/RAST_NESREG_2024.12.01.AV.TEMP.1DAY.000082150.RData") +data = terra::rast(masked.raster) +#not sure why here didn't work but ok + +#load shapefiles +#EPUs <- read_sf("data-raw/shp_files/adv rep MAB GOM GBK NES SCSPoly.shp") +#index.raster <- masked.raster +#index.raster[] <- 1:ncell(index.raster) #not sure what this is for, but there are 5120 cells apparently + +shp.vect = terra::vect("data-raw/shp_files/adv rep MAB GOM GBK NES SCSPoly.shp") +shp.str = as.data.frame(shp.vect) + +#pull out EPU names +area.names <- as.vector(shp.str$Primary) +which.att = which(apply(shp.str,2,function(x) all(area.names %in% x))) +which.area = match(area.names,shp.str[,which.att]) +#testing using terra instead +avg_sst <- c() +for(j in 1:length(area.names)){ + area.data = terra::mask(data,shp.vect[which.area[j],], touches = T) + area.agg = terra::global(area.data,fun = mean,na.rm = T) + area.agg <- cbind(area.agg$mean,area.names[j]) + avg_sst <- rbind(avg_sst,area.agg) +} + +#make into tidy data frame +avg_sst <- as.data.frame(avg_sst) %>% + dplyr::rename(SST = V1, EPU = V2)