0))
+ ACLIM_weekly_hind[rr,]$warming_level<- GWLS[i,]$warming_level
+ rm(rr)
+}
+
+ggplot(ACLIM_weekly_hind)+
+ geom_point(aes(x=))
+
+
+
+
+
+
+hind <- ACLIM_seasonal_hind%>%filter(var%in%plot_vars)%>%
+ select(all_of(c(grplist,"year","season","mn_val","mnVal_hind","sdVal_hind","mnDate","jday")))
+tt <- hind%>%ungroup()%>%filter(year<2019)%>%
+ group_by(!!!syms(c(grplist,"season")))%>%
+ summarise(mnVal_hind = mean(mn_val,na.rm=T),
+ sdVal_hind = sd(mn_val,na.rm=T))
+hind <- hind%>%rename(mnVal_hind_old = mnVal_hind,sdVal_hind_old =sdVal_hind)%>%left_join(tt)
+rm(tt)
+
+hind_mo <- ACLIM_monthly_hind%>%filter(var%in%plot_vars)%>%rename(units = units.x)%>%
+ select(all_of(c(grplist,"year","season","mo","sdVal_hind","mn_val","mnVal_hind","mnDate","jday")))
+tt <- hind_mo%>%ungroup()%>%filter(year<2019)%>%
+ group_by(!!!syms(c(grplist,"season","mo")))%>%
+ summarise(mnVal_hind = mean(mn_val,na.rm=T),
+ sdVal_hind = sd(mn_val,na.rm=T))
+hind_mo <- hind_mo%>%rename(mnVal_hind_old = mnVal_hind)%>%left_join(tt)
+
+ggplot(hind%>%filter(var=="aice",basin=="SEBS",season=="Spring"))+
+ geom_line(aes(x=year,y=mnVal_hind_old,color="old"))+
+ geom_line(aes(x=year,y=mnVal_hind,color="new"))
+
+ggplot(hind_mo%>%filter(var=="aice",basin=="SEBS",mo==12))+
+ geom_line(aes(x=year,y=mnVal_hind_old,color="old"))+
+ geom_line(aes(x=year,y=mnVal_hind,color="new"))+
+ geom_line(aes(x=year,y=mnVal_hind+sdVal_hind))+
+ geom_line(aes(x=year,y=mnVal_hind-sdVal_hind))
+
+hind <- ACLIM_seasonal_hind%>%filter(var%in%plot_vars)%>%
+ select(all_of(c(grplist,"year","season","sd_val","mn_val","val_raw","mnVal_hind","mnDate","jday")))
+tt <- hind%>%ungroup()%>%filter(year<2019)%>%
+ group_by(!!!syms(c(grplist,"season")))%>%
+ summarise(mnVal_hind = mean(mn_val,na.rm=T),
+ sdVal_hind = sd(mn_val,na.rm=T))
+hind <- hind%>%rename(mnVal_hind_old = mnVal_hind)%>%left_join(tt)
+rm(tt)
+
+hind_mo <- ACLIM_monthly_hind%>%filter(var%in%plot_vars)%>%
+ select(all_of(c(grplist,"year","season","mo","sd_val","mn_val","val_raw","mnVal_hind","mnDate","jday")))
+ tt <- hind_mo%>%ungroup()%>%filter(year<2019)%>%
+ group_by(!!!syms(c(grplist,"season","mo")))%>%
+ summarise(mnVal_hind = mean(mn_val,na.rm=T),
+ sdVal_hind = sd(mn_val,na.rm=T))
+hind_mo <- hind_mo%>%rename(mnVal_hind_old = mnVal_hind)%>%left_join(tt)
+
+
+
+
diff --git a/R/sub_scripts/load_maps.R b/R/sub_scripts/load_maps.R
old mode 100644
new mode 100755
diff --git a/R/sub_scripts/make_CEATTLE_dat_files.R b/R/sub_scripts/make_CEATTLE_dat_files.R
new file mode 100755
index 0000000..1db0ff5
--- /dev/null
+++ b/R/sub_scripts/make_CEATTLE_dat_files.R
@@ -0,0 +1,297 @@
+#'
+#'
+#'make_CEATTLE_dat_files.R
+#'
+#'first run "R/sub_scripts/make_ceattle_indices.R"
+#'
+
+# "Data/out/CEATTLE_indices/ceattle_vars_wide_op.Rdata"
+
+ # Load setup
+ #--------------------------------------
+ suppressMessages(source("R/make.R"))
+
+ # load the datafile:
+ load(file="Data/out/CEATTLE_indices/ceattle_vars_op.Rdata")
+ load(file="Data/out/CEATTLE_indices/ceattle_vars_wide_op.Rdata")
+
+ # switches
+ varall <- varlist <- unique(ceattle_vars_op$var)
+ shortlist <- c("Summer_temp_bottom5m","Winter_temp_bottom5m",
+ "Winter_pH_depthavg","Summer_oxygen_bottom5m",
+ "Summer_temp_surface5m",
+ "Spring_Cop_integrated",
+ "Fall_largeZoop_integrated")
+ thisYr <- as.numeric(format(Sys.time(), "%Y"))
+ today <- format(Sys.time(), "%b %d, %Y")
+ lastyr_hind <- thisYr #as.numeric(thisYr) #2021
+ hind_yrs <- 1979:(lastyr_hind) # define the years of your estimation model (noting that recruitment covars are for cohort year)
+ fut_yrs <- (lastyr_hind+1):2099 # define the years of your projections
+ stitchDate <- "2019-12-30" # last date of the ACLIM hindcast
+ stitchDate_op <- "2022-05-16" #last operational hindcast date
+ log_adj <- 1e-4
+ zscore_years <- 1980:2010 # years to recenter z score on
+ plotbasin <- "SEBS"
+
+ # Define the name for the .dat file
+ fn <- paste0("ACLIM2_CMIP6_short_delta_bc",thisYr,".dat")
+ fn_long <- paste0("ACLIM2_CMIP6_delta_bc",thisYr,".dat")
+
+ archive_old <- T # Archive the older version of the .dat file?
+ normlist <- read.csv(file=file.path(Rdata_path,"../normlist.csv"))
+
+ CMIPS <- c("K20P19_CMIP6","K20P19_CMIP5")
+
+ # create outpath folder
+ outpath <- "Data/out/ADMB_datfiles"
+ if(!dir.exists(outpath)) dir.create(outpath)
+
+ # define hind and fut data files
+ fndat_hind <- file.path(outpath,paste("KKHhind_",fn,sep=""))
+ fndat_fut <- file.path(outpath,paste("KKHfut_",fn,sep=""))
+
+ fndat_hind_long <- file.path(outpath,paste("KKHhind_",fn_long,sep=""))
+ fndat_fut_long <- file.path(outpath,paste("KKHfut_",fn_long,sep=""))
+
+ fndat_hind2 <- file.path(outpath,paste("hind_",fn,sep=""))
+ fndat_fut2 <- file.path(outpath,paste("fut_",fn,sep=""))
+
+ # create and archive .dat files
+ outfile <- fndat_fut
+ if(file.exists(outfile)&archive_old){
+ # archive older version
+ archivefl <- paste0(substr(outfile,start=1,stop=nchar(outfile)-4),
+ format(Sys.time(), "%Y%m%d_%H%M%S"),".dat")
+ file.rename(outfile, archivefl)
+ #file.remove(outfile)
+ }
+
+ file.create(outfile)
+ outfile <- fndat_hind
+ if(file.exists(outfile)&archive_old){
+ # archive older version
+ archivefl <- paste0(substr(outfile,start=1,stop=nchar(outfile)-4),
+ format(Sys.time(), "%Y%m%d_%H%M%S"),".dat")
+ file.rename(outfile, archivefl)
+ #file.remove(outfile)
+ }
+ file.create(outfile)
+
+ # ----------------------------------------------
+ # CREATE INDICES
+ # ----------------------------------------------
+ # 2 -- rescale (Z-score) data and get variables
+
+ # preview possible variables
+ varall
+
+ grpby2 <- c("type","var","basin",
+ "year","sim","sim_type",
+ "bc")
+
+ #define hind and fut:
+ hind <- ceattle_vars_op%>%filter(year%in%hind_yrs,sim_type=="hind")%>%
+ ungroup()
+
+ mnhind <- hind%>%
+ group_by(season,var,basin)%>%
+ summarize(mnhind = mean(val_use,na.rm=T),
+ sdhind = sd(val_use, na.rm=T))%>%ungroup()
+ hind <- hind%>%left_join(mnhind)%>%
+ mutate(val_use_scaled = (val_use-mnhind)/sdhind,
+ zscoreyrs = paste0( zscore_years[1],":", rev(zscore_years)[1]))%>%
+ select(all_of(c(grpby2,"val_use","val_use_scaled")))%>%distinct()%>%
+ ungroup()
+
+ hind <- hind%>%full_join(expand.grid(year= hind_yrs, var = unique(hind$var))) %>%ungroup()
+
+ fut <- ceattle_vars_op%>%
+ filter(year%in%fut_yrs)%>%ungroup()
+
+ #temporary fix: Add Large Zoop for CESM RCP85
+ tmp_fut <- fut%>%filter(GCM_scen=="cesm_rcp85")
+ grplist <- c("NCaS_integrated","EupS_integrated")
+
+ if(length(grep("largeZoop_integrated",unique(tmp_fut$var_raw)) )==0){
+ tmp_futA <- tmp_fut[grep(grplist[1],tmp_fut$var),]
+ tmp_futB <- tmp_fut[grep(grplist[2],tmp_fut$var),]
+ tmp_futA$var2 <- grplist[1]
+ tmp_futB$var2 <- grplist[2]
+ sumat<-c("val_use","mnVal_hind","val_delta","val_biascorrected",
+ "val_raw","mn_val")
+ tmp_var_zoop <- rbind(tmp_futA,tmp_futB)%>%
+ dplyr::filter(var2%in%c("NCaS_integrated","EupS_integrated"))%>%
+ dplyr::group_by(
+ var,
+ season,
+ type,
+ basin,
+ year,
+ sim,
+ gcmcmip,
+ GCM,
+ scen,
+ sim_type,
+ bc,
+ GCM_scen,
+ GCM_scen_sim,
+ CMIP,
+ GCM2,
+ GCM2_scen_sim,
+ jday,mnDate,var_raw,lognorm,var2)%>%
+ summarise_at(all_of(sumat),mean,na.rm=T)%>%
+ mutate(var_raw="largeZoop_integrated",
+ var = paste0(season,"_largeZoop_integrated"))%>%select(-var2)%>%
+ summarise_at(all_of(sumat),sum,na.rm=T)%>%relocate(names(fut))%>%
+ distinct()
+
+ fut <- rbind(fut%>%ungroup(), tmp_var_zoop%>%ungroup())
+
+ }
+ fut <- fut%>%left_join(mnhind)%>%
+ mutate(val_use_scaled = (val_use-mnhind)/sdhind,
+ zscoreyrs = paste0( zscore_years[1],":", rev(zscore_years)[1]))%>%ungroup()
+
+ data.frame(fut%>%filter(var=="Winter_largeZoop_integrated",GCM_scen=="cesm_rcp85"))
+ data.frame(fut%>%ungroup()%>%group_by(GCM_scen)%>%summarise(count = length(var)))
+
+ fut <- fut%>%full_join(expand.grid(year= fut_yrs,
+ var = unique(fut$var),
+ GCM_scen = unique(fut$GCM_scen)))
+ data.frame(fut%>%ungroup()%>%group_by(GCM_scen)%>%summarise(count = length(var)))
+ data.frame(fut%>%filter(var=="Winter_largeZoop_integrated",GCM_scen=="cesm_rcp85"))
+
+ # now identify which covars are highly correlated
+
+ d_wide <- ceattle_vars_op%>%
+ filter(year<=lastyr_hind)%>%
+ left_join(mnhind)%>%
+ mutate(val_use_scaled = (val_use-mnhind)/sdhind)%>%ungroup()%>%
+ group_by(across(all_of(grpby2)))%>%
+ summarize_at(all_of(c("val_use_scaled")), mean, na.rm=T)%>%
+ tidyr::pivot_wider(names_from = "var", values_from = "val_use_scaled")
+
+ # calculate correlations and display in column format
+
+ # define columns with meta data:
+ col_meta <- which(colnames(d_wide)%in%grpby2)
+ d_wide_dat <-d_wide[,-col_meta]
+ num_col <- ncol(d_wide[,-col_meta])
+ out_indx <- which(upper.tri(diag(num_col)))
+ cor_cols <- d_wide_dat%>%
+ do(melt(cor(.,
+ method="spearman",
+ use="pairwise.complete.obs"),
+ value.name="cor")[out_indx,])
+
+ corr <- cor(na.omit(d_wide_dat))
+
+ long_dat <- reshape2::melt(corr,variable.name = "variable") %>%
+ as.data.frame()
+
+ # plot co variation between variables
+ corplot <- long_dat %>%arrange(value)%>%
+ ggplot(aes(x=Var1, y=Var2, fill=value)) +
+ geom_raster() +
+ scale_fill_viridis_c()+
+ theme_minimal()+
+ theme(axis.text.x = element_text(angle = 90))
+
+ jpeg(filename = file.path("Data/out/CEATTLE_indices/CEATTLE_indicescorplot.jpg"),
+ width=8,height=7,units="in",res=350)
+ print(corplot)
+ dev.off()
+
+ # # remove those where cov is high (temp by season and cold pool by season)
+ # subset <- long_dat$>$filter(abs(value)<0.6)
+
+ # 3 -- write data to hind .dat file
+ # ------------------------------------
+
+
+ # CEATTLE uses a spp overlap index - you can skip this
+
+ overlapdat <- data.frame(
+ atf_OL=c(0.9391937,0.8167094,0.808367,0.5926875,0.7804481,0.5559549,
+ 0.4006931,0.5881404,0.7856776,0.511565,0.6352048,0.5583476,
+ 0.5792738,0.5417657,0.8212887,0.6287613,0.4536608,0.6587292,
+ 0.4884194,0.8289379,0.4399257,0.5950167,0.6388434,0.6111834,
+ 0.8742649,0.7868746,0.8024257,0.6227457,0.4956742,0.4347917,
+ 0.4791108,0.4369006,0.5613625,0.4353015),
+ south_OL=c(0.9980249,0.9390368,0.9959974,0.6130846,0.951234,0.5851891,
+ 0.4934879,0.641471,0.9809618,0.5596813,0.7196964,0.6754698,
+ 0.5774808,0.6041351,0.9406521,0.7949525,0.5306435,0.7977694,
+ 0.5345031,0.9879945,0.5079171,0.7148121,0.8997132,0.7340859,
+ 0.9962068,0.9627235,0.998043,0.8111,0.6087638,0.513057,0.5492621,
+ 0.4971361,0.665453,0.5969653)
+ )
+
+
+ includeOverlap <- F
+ overlap <- matrix(1,3,length(sort(unique(hind$year))))
+ overlap_fut <- array(1,c(3,length(unique(fut$GCM_scen))+1,length(sort(unique(fut$year)))))
+ if(includeOverlap){
+ overlap[3,] <- overlapIN
+ overlap[3,][overlap[3,]>1]<-1 #covs$BT2to6/covs$BT0to6
+ }
+
+ # replace NA values with the mean
+ hind_short <- hind%>%mutate(long_name=var)%>%
+ filter(var%in%shortlist)
+
+ cat("making short hind and future dat files\n")
+
+
+ # Kir's .dat file
+ makeDat_hind(datIN = hind_short,
+ outfile = fndat_hind,
+ value2use = "val_use",
+ value2use_scaled = "val_use_scaled",
+ NAVal = "mean",
+ nsppIN = 3,
+ overlapIN = overlap,
+ nonScaled_covlist = c("Summer_temp_bottom5m","Summer_temp_surface5m"),
+ Scaled_covlist = unique(hind_short$var))
+
+ makeDat_fut( datIN = fut%>%mutate(long_name=var)%>%
+ filter(var%in%shortlist),
+ hinddatIN = hind%>%mutate(long_name=var)%>%
+ filter(var%in%shortlist),
+ outfile = fndat_fut,
+ value2use = "val_use",
+ value2use_scaled = "val_use_scaled",
+ NAVal = "mean",
+ nsppIN = 3,
+ last_nyrs_avg = 10,
+ overlapIN = overlap_fut, #(nspp,nsim+1,nyrs_fut)
+ overlap_hind = overlap,
+ nonScaled_covlist = c("Summer_temp_bottom5m","Summer_temp_surface5m"),
+ Scaled_covlist = unique(hind_short$var))
+
+ # Kir's .dat file
+ cat("making long hind and future dat files\n")
+
+ makeDat_hind(datIN = hind%>%mutate(long_name=var),
+ outfile = fndat_hind_long,
+ value2use = "val_use",
+ value2use_scaled = "val_use_scaled",
+ NAVal = "mean",
+ nsppIN = 3,
+ overlapIN = overlap,
+ nonScaled_covlist = c("Summer_temp_bottom5m","Summer_temp_surface5m"),
+ Scaled_covlist = unique(hind$var))
+
+ makeDat_fut( datIN = fut%>%mutate(long_name=var),
+ hinddatIN = hind%>%mutate(long_name=var),
+ outfile = fndat_fut_long,
+ value2use = "val_use",
+ value2use_scaled = "val_use_scaled",
+ NAVal = "mean",
+ nsppIN = 3,
+ last_nyrs_avg = 10,
+ overlapIN = overlap_fut, #(nspp,nsim+1,nyrs_fut)
+ overlap_hind = overlap,
+ nonScaled_covlist = c("Summer_temp_bottom5m","Summer_temp_surface5m"),
+ Scaled_covlist = unique(hind$var))
+
+ message(paste0("dat files successfully created (e.g., ",fndat_fut_long,")"))
diff --git a/R/sub_scripts/make_NRS_indices.R b/R/sub_scripts/make_NRS_indices.R
old mode 100644
new mode 100755
diff --git a/R/sub_scripts/make_NRS_indices_raw.R b/R/sub_scripts/make_NRS_indices_raw.R
old mode 100644
new mode 100755
diff --git a/R/sub_scripts/make_RKC_Indices_Andre.R b/R/sub_scripts/make_RKC_Indices_Andre.R
new file mode 100755
index 0000000..7b3337f
--- /dev/null
+++ b/R/sub_scripts/make_RKC_Indices_Andre.R
@@ -0,0 +1,450 @@
+#'
+#'
+#'
+#'make_RKC_Indices_Andre.R
+#'
+#'This script generates the indices for the NRS and pcod papers
+#'for Punt et al. 2024
+#'
+#' staion info for
+#' I-11
+#' K-14
+#' Z-05
+#' D-10
+#'
+#'
+#'
+
+# rm(list=ls())
+
+# load ACLIM packages and functions
+ suppressMessages(source("R/make.R"))
+
+# Set up stations
+ RKC_immature_males_stations <- c(
+ "D-10 ",
+ "E-12 ",
+ "F-11 ",
+ "G-10 ",
+ "H-10 ",
+ "H-11 ",
+ "H-13 ",
+ "I-10 ",
+ "I-11 ",
+ "I-12 ",
+ "J-11 ",
+ "J-13 "
+ )
+
+
+ TC_small_males_stations <- c(
+ "A-02 ",
+ "C-18 ",
+ "D-18 ",
+ "F-19 ",
+ "F-24 ",
+ "G-25 ",
+ "I-21 ",
+ "I-26 ",
+ "I-20 ",
+ "H-19 ",
+ "J-21 ",
+ "I-20 ",
+ "L-27 ",
+ "M-26 ")
+
+
+# set up paths
+ outfldr <- "Data/out/RKC_indices"
+ if(!dir.exists(outfldr)) dir.create(outfldr)
+
+# specific the stations, variables, and projection sets
+ station_set <- c("D-10 ","Z-05 ", "K-14 ","I-11 ")
+ varlist <- c("pH_bottom5m","pH_depthavg","pH_integrated","pH_surface5m","temp_surface5m", "temp_bottom5m")
+ CMIPset <- c("K20P19_CMIP6","K20P19_CMIP5")
+ stitchDate_op <- "2021-12-30" #last operational hindcast date
+
+# preview possible variables
+ load(file = "Data/out/weekly_vars.Rdata")
+ load(file = "Data/out/srvy_vars.Rdata")
+
+# save setup for plotting
+ save(list=ls(),file=file.path(outfldr,"RKC_indices_setup.Rdata"))
+
+
+# ------------------------------------
+# Hindcast First
+# ------------------------------------
+
+ # Load and select stations (survey replicated hindcast):
+
+ load("Data/out/K20P19_CMIP6/BC_ACLIMsurveyrep/ACLIMsurveyrep_B10K-K20P19_CORECFS_BC_hind.Rdata")
+
+ # first for RCK Immature Males
+ # preview the data
+ station_set <- RKC_immature_males_stations
+ head(hind%>%filter(station_id%in%station_set,var%in%varlist)%>%data.frame())
+
+ # select the subset for RKC
+ hind_RKC_immature_males <- hind%>%filter(station_id%in%station_set,var%in%varlist)%>%
+ mutate(val_use = val_raw)%>%ungroup()
+
+ # output the data as Rdata and CSV
+ write_csv(hind_RKC_immature_males, file.path(outfldr,"hind_RKC_immature_males.csv"))
+ save(hind_RKC_immature_males, file=file.path(outfldr,"hind_RKC_immature_males.Rdata"))
+
+ # Then for TC_small_males
+ # preview the data
+ station_set <- TC_small_males_stations
+ head(hind%>%filter(station_id%in%station_set,var%in%varlist)%>%data.frame())
+
+ # select the subset for RKC
+ hind_TC_small_males <- hind%>%filter(station_id%in%station_set,var%in%varlist)%>%
+ mutate(val_use = val_raw)%>%ungroup()
+
+ # output the data as Rdata and CSV
+ write_csv(hind_TC_small_males, file.path(outfldr,"hind_TC_small_males.csv"))
+ save(hind_TC_small_males, file=file.path(outfldr,"hind_TC_small_males.Rdata"))
+
+ # ------------------------------------
+ # now get monthly values:
+ # ------------------------------------
+ # SST, Ph and SBT values. predictions for 1/15, 2/15 .. 12/15..
+ # Load bias corrected weekly strata data:
+
+ rm(hind)
+ load("Data/out/K20P19_CMIP6/BC_ACLIMregion/ACLIMregion_B10K-K20P19_CORECFS_BC_hind.Rdata")
+
+ for(outType in c("RKC_immature_males","TC_small_males")){
+ if(outType == "TC_small_males")
+ hindDat <- hind_TC_small_males
+
+ if(outType == "RKC_immature_males")
+ hindDat <- hind_RKC_immature_males
+
+ unique(hindDat$strata)
+ unique(hind$mo)
+ unique(hind$year)
+ tmp_mat <- expand_grid(year =unique(hind$year), mo = unique(hind$mo))%>%
+ rowwise()%>%mutate(targetdate = as.Date(paste0(year,"-",mo,"-",15)))
+ tt <- strptime(tmp_mat$targetdate,format = "%Y-%m-%d")
+ tmp_mat$target_jday <- tt$yday+1
+ tmp_mat$target_mday <- tt$mday
+
+ hindDat_monthly <- hind%>%
+ filter(strata%in%unique(hindDat$strata),
+ var%in%varlist)%>%
+ mutate(val_use = val_raw)%>%left_join(tmp_mat)%>%
+ rowwise()%>%mutate(mday=strptime(mnDate,format = "%Y-%m-%d")$mday)%>%ungroup()%>%
+ mutate(keepA = mday>=target_mday-3,keepB = mday<=target_mday+3)%>%
+ filter(keepA,keepB)%>%
+ ungroup()
+
+ grp_names <- names(hindDat_monthly)[!names(hindDat_monthly)%in%c("strata","strata_area_km2", "keepA","keepB")]
+ nms <- c("val_use","val_raw","mn_val","sd_val",
+ "n_val","jday","mnDate",
+ "mnVal_hind","sdVal_hind","nVal_hind",
+ "seVal_hind","sdVal_hind_mo", "sdVal_hind_yr")
+
+ nms2 <- grp_names[!grp_names%in%nms]
+
+ hindDat_monthly_avg <- hindDat_monthly%>%
+ select(!!!syms(grp_names))%>%group_by(!!!syms(nms2))%>%
+ summarize_at(nms,mean,na.rm=T)%>%ungroup()
+
+ if(outType == "TC_small_males"){
+ hind_TC_small_males_monthly <- hindDat_monthly
+ hind_TC_small_males_monthly_avg <- hindDat_monthly_avg
+
+ # output the data as Rdata and CSV
+ write_csv(hind_TC_small_males_monthly, file.path(outfldr,"hind_TC_small_males_monthly.csv"))
+ save(hind_TC_small_males_monthly, file=file.path(outfldr,"hind_TC_small_males_monthly.Rdata"))
+
+ write_csv(hind_TC_small_males_monthly_avg, file.path(outfldr,"hind_TC_small_males_monthly_avg.csv"))
+ save(hind_TC_small_males_monthly_avg, file=file.path(outfldr,"hind_TC_small_males_monthly_avg.Rdata"))
+ }
+
+ if(outType == "RKC_immature_males"){
+ hind_RKC_immature_males_monthly <- hindDat_monthly
+ hind_RKC_immature_males_monthly_avg <- hindDat_monthly_avg
+ # output the data as Rdata and CSV
+ write_csv(hind_RKC_immature_males_monthly, file.path(outfldr,"hind_RKC_immature_males_monthly.csv"))
+ save(hind_RKC_immature_males_monthly, file=file.path(outfldr,"hind_RKC_immature_males_monthly.Rdata"))
+
+ write_csv(hind_RKC_immature_males_monthly_avg, file.path(outfldr,"hind_RKC_immature_males_monthly_avg.csv"))
+ save(hind_RKC_immature_males_monthly_avg, file=file.path(outfldr,"hind_RKC_immature_males_monthly_avg.Rdata"))
+ }
+ rm(tmp_mat)
+ rm(hindDat_monthly)
+ rm(hindDat_monthly_avg)
+
+ }
+
+# ------------------------------------
+# Projections Next
+# ------------------------------------
+ select_list <- unique(c(names(hind_RKC_immature_males),"long_name","qry_date","station_id",
+ "sim_type","GCM","RCP","mod", "CMIP", "val_raw","mnVal_hind","sdVal_hind","val_delta"))
+ jj <- 0
+
+for (CMIP in c("CMIP6","CMIP5")){
+
+ if(CMIP == "CMIP6"){
+
+ fl_base_srv <- "ACLIMsurveyrep_B10K-K20P19_CMIP6_"
+ fl_path_srv <- "Data/out/K20P19_CMIP6/BC_ACLIMsurveyrep"
+
+ fl_base_reg <- "ACLIMregion_B10K-K20P19_CMIP6_"
+ fl_path_reg <- "Data/out/K20P19_CMIP6/BC_ACLIMregion"
+
+
+ # CMIP6
+ ssps <- c("ssp126","ssp585")
+ gcms <- c("cesm","miroc","gfdl")
+ set <- expand.grid(gcms,ssps,stringsAsFactors = F)
+
+ }
+
+ if(CMIP == "CMIP5"){
+ #now get monthly
+ fl_base_srv <- "ACLIMsurveyrep_B10K-K20P19_CMIP5_"
+ fl_path_srv <- "Data/out/K20P19_CMIP5/BC_ACLIMsurveyrep"
+
+ fl_base_reg <- "ACLIMregion_B10K-K20P19_CMIP5_"
+ fl_path_reg <- "Data/out/K20P19_CMIP5/BC_ACLIMregion"
+
+ # now for CMIP5
+ gcms <- c("CESM","MIROC","GFDL")
+ ssps <- c("rcp45","rcp85")
+ set <- expand.grid(gcms,ssps,stringsAsFactors = F)
+
+ }
+
+ for(outType in c("RKC_immature_males","TC_small_males")){
+
+ eval(parse(text = paste0("station_set <- ",outType,"_stations") ))
+
+ # for each scenario
+ for(i in 1:dim(set)[1]){
+
+ cat("set = ",set[i,1],set[i,2],"\n")
+ fl <- paste0(fl_base_srv,set[i,1],"_",set[i,2],"_BC_fut.Rdata")
+ load(file.path(fl_path_srv,fl))
+
+ tmp <- fut%>%filter(station_id%in%station_set,var%in%varlist)%>%data.frame()
+
+ my_vars <- function() select_list[select_list%in%names(tmp)]
+
+ tmp <- tmp%>%ungroup()%>%
+ dplyr::select(my_vars())%>%
+ mutate(val_use = val_delta)%>%data.frame()
+
+ if(i==1){
+ fut_RKC_all <- tmp
+ }else{
+ fut_RKC_all <- rbind(fut_RKC_all,tmp)
+ }
+ rm(fut)
+
+ #now get monthly
+ fl <- paste0(fl_base_reg,set[i,1],"_",set[i,2],"_BC_fut.Rdata")
+ load(file.path(fl_path_reg,fl))
+
+ tmp_mat <- expand_grid(year =unique(tmp$year), mo = unique(substr(fut$mnDate,6,7)))%>%
+ rowwise()%>%mutate(targetdate = as.Date(paste0(year,"-",mo,"-",15)))
+ tt <- strptime(tmp_mat$targetdate,format = "%Y-%m-%d")
+ tmp_mat$target_jday <- tt$yday+1
+ tmp_mat$target_mday <- tt$mday
+
+ futDat_monthly <- fut%>%
+ filter(strata%in%unique(tmp$strata),
+ var%in%varlist)%>%mutate(mo = substr(mnDate,6,7))%>%
+ mutate(val_use = val_delta)%>%left_join(tmp_mat)%>%rowwise()%>%
+ rowwise()%>%mutate(mday=strptime(mnDate,format = "%Y-%m-%d")$mday)%>%ungroup()%>%
+ mutate(keepA = mday>=target_mday-3,keepB = mday<=target_mday+3)%>%
+ filter(keepA,keepB)%>%
+ ungroup()
+
+ if(i==1){
+ fut_RKC_all_monthly <- futDat_monthly
+ }else{
+ fut_RKC_all_monthly <- rbind(fut_RKC_all_monthly,futDat_monthly)
+ }
+ rm(futDat_monthly)
+ rm(fut)
+
+ rm(tmp)
+
+ }
+
+ jj <- jj + 1
+ if(jj==1){
+ fut_RKC_all_tmp <- fut_RKC_all
+ fut_RKC_all_tmp_monthly <- fut_RKC_all_monthly
+ }else{
+ fut_RKC_all <- rbind(fut_RKC_all_tmp,fut_RKC_all)
+ fut_RKC_all_monthly <- rbind(fut_RKC_all_tmp_monthly,fut_RKC_all_monthly)
+
+ }
+
+ # fix the caps issue:
+ fut_RKC <- fut_RKC_all%>%
+ mutate(GCM = gsub("MIROC","miroc",GCM))%>%
+ mutate(GCM = gsub("GFDL","gfdl",GCM))%>%
+ mutate(GCM = gsub("CESM","cesm",GCM))%>%
+ mutate(sim = gsub("MIROC","miroc",sim))%>%
+ mutate(sim = gsub("GFDL","gfdl",sim))%>%
+ mutate(sim = gsub("CESM","cesm",sim))%>%
+ mutate(GCM = factor(GCM, levels =c("hind","gfdl","cesm","miroc")))%>%
+ mutate(GCM_scen = paste0(GCM,"_",RCP))%>%data.frame()
+
+ # fix the caps issue:
+ fut_RKC_monthly <- fut_RKC_all_monthly%>%
+ mutate(GCM = gsub("MIROC","miroc",GCM))%>%
+ mutate(GCM = gsub("GFDL","gfdl",GCM))%>%
+ mutate(GCM = gsub("CESM","cesm",GCM))%>%
+ mutate(sim = gsub("MIROC","miroc",sim))%>%
+ mutate(sim = gsub("GFDL","gfdl",sim))%>%
+ mutate(sim = gsub("CESM","cesm",sim))%>%
+ mutate(GCM = factor(GCM, levels =c("hind","gfdl","cesm","miroc")))%>%
+ mutate(GCM_scen = paste0(GCM,"_",RCP))%>%data.frame()
+
+
+ if(outType =="RKC_immature_males"){
+ fut_RKC_immature_males <- fut_RKC
+ # output the data as Rdata and CSV
+ write_csv(fut_RKC_immature_males, file.path(outfldr,"fut_RKC_immature_males.csv"))
+ save(fut_RKC_immature_males, file=file.path(outfldr,"fut_RKC_immature_males.Rdata"))
+
+ fut_RKC_immature_males_monthly <- fut_RKC_monthly
+ # output the data as Rdata and CSV
+ write_csv(fut_RKC_immature_males_monthly, file.path(outfldr,"fut_RKC_immature_males_monthly.csv"))
+ save(fut_RKC_immature_males_monthly, file=file.path(outfldr,"fut_RKC_immature_males_monthly.Rdata"))
+
+
+ }else{
+ if(outType =="TC_small_males"){
+
+ fut_TC_small_males <- fut_RKC
+ # output the data as Rdata and CSV
+ write_csv(fut_TC_small_males, file.path(outfldr,"fut_TC_small_males.csv"))
+ save(fut_TC_small_males, file=file.path(outfldr,"fut_TC_small_males.Rdata"))
+
+ fut_TC_small_males_monthly <- fut_RKC_monthly
+ # output the data as Rdata and CSV
+ write_csv(fut_TC_small_males_monthly, file.path(outfldr,"fut_TC_small_males_monthly.csv"))
+ save(fut_TC_small_males_monthly, file=file.path(outfldr,"fut_TC_small_males_monthly.Rdata"))
+
+
+ }else{
+ stop("outType doesn't match options")
+ }
+ }
+
+
+ rm(fut_RKC)
+ rm(fut_RKC_monthly)
+ }
+
+}
+
+ if(1 == 10){
+ load(file=file.path(outfldr,"fut_RKC_immature_males.Rdata")) #fut_RKC_immature_males
+ load(file.path(outfldr,"fut_TC_small_males.Rdata")) #fut_TC_small_males
+
+ load(file=file.path(outfldr,"hind_RKC_immature_males.Rdata")) #hind_RKC_immature_males
+ load(file.path(outfldr,"hind_TC_small_males.Rdata")) #hind_TC_small_males
+ }
+
+ # now summarize across stations
+ AVG_hind_RKC_immature_males <- hind_RKC_immature_males%>%ungroup()%>%
+ group_by(var,units,year,basin,sim,type,sim_type)%>%
+ summarize(mn_val_use = mean(val_use,na.rm=T),
+ sd_val_use = sd(val_use,na.rm=T),
+ n_val_use = length.na(val_use),
+ mnDate = mean(mnDate, na.rm =T))%>%data.frame()
+ AVG_hind_TC_small_males <- hind_TC_small_males%>%ungroup()%>%
+ group_by(var,units,year,basin,sim,type,sim_type)%>%
+ summarize(mn_val_use = mean(val_use,na.rm=T),
+ sd_val_use = sd(val_use,na.rm=T),
+ n_val_use = length.na(val_use),
+ mnDate = mean(mnDate, na.rm =T))%>%data.frame()
+
+ AVG_fut_RKC_immature_males <- fut_RKC_immature_males%>%ungroup()%>%
+ group_by(var,units,year,basin,sim,type,sim_type,GCM,RCP,mod,CMIP,GCM_scen)%>%
+ summarize(mn_val_use = mean(val_use,na.rm=T),
+ sd_val_use = sd(val_use,na.rm=T),
+ n_val_use = length.na(val_use),
+ mnDate = mean(mnDate, na.rm =T))%>%data.frame()
+ AVG_fut_TC_small_males <- fut_TC_small_males%>%ungroup()%>%
+ group_by(var,units,year,basin,sim,type,sim_type,GCM,RCP,mod,CMIP,GCM_scen)%>%
+ summarize(mn_val_use = mean(val_use,na.rm=T),
+ sd_val_use = sd(val_use,na.rm=T),
+ n_val_use = length.na(val_use),
+ mnDate = mean(mnDate, na.rm =T))%>%data.frame()
+
+
+ # average across strata monhtly
+
+ futDat_monthly <- fut_TC_small_males_monthly
+
+ grp_names <- names(futDat_monthly)[!names(futDat_monthly)%in%c("strata","strata_area_km2", "keepA","keepB")]
+ nms <- c("val_use","val_raw","mnVal_hind","sdVal_hind",
+ "sdVal_hind_mo","sdVal_hind_yr","mnVal_hist",
+ "sdVal_hist","sdVal_hist_mo","sdVal_hist_yr",
+ "sf_wk","val_biascorrectedyr" ,"val_biascorrectedmo",
+ "sf_mo","val_delta_adj","val_biascorrected",
+ "val_biascorrectedwk", "val_delta",
+ "sf_yr","val_use" )
+
+ nms2 <- grp_names[!grp_names%in%nms]
+
+ futDat_monthly_avg <- futDat_monthly%>%
+ select(!!!syms(grp_names))%>%group_by(!!!syms(nms2))%>%
+ summarize_at(nms,mean,na.rm=T)%>%ungroup()
+
+ fut_TC_small_males_monthly_avg <- futDat_monthly_avg
+
+ futDat_monthly <- fut_RKC_immature_males_monthly
+
+ rm(futDat_monthly_avg)
+ grp_names <- names(futDat_monthly)[!names(futDat_monthly)%in%c("strata","strata_area_km2", "keepA","keepB")]
+ nms <- c("val_use","val_raw","mnVal_hind","sdVal_hind",
+ "sdVal_hind_mo","sdVal_hind_yr","mnVal_hist",
+ "sdVal_hist","sdVal_hist_mo","sdVal_hist_yr",
+ "sf_wk","val_biascorrectedyr" ,"val_biascorrectedmo",
+ "sf_mo","val_delta_adj","val_biascorrected",
+ "val_biascorrectedwk", "val_delta",
+ "sf_yr","val_use" )
+
+ nms2 <- grp_names[!grp_names%in%nms]
+
+ futDat_monthly_avg <- futDat_monthly%>%
+ select(!!!syms(grp_names))%>%group_by(!!!syms(nms2))%>%
+ summarize_at(nms,mean,na.rm=T)%>%ungroup()
+
+ fut_RKC_immature_males_monthly_avg <- futDat_monthly_avg
+
+
+ write_csv(AVG_hind_TC_small_males, file.path(outfldr,"AVG_hind_TC_small_males.csv"))
+ save(AVG_hind_TC_small_males, file=file.path(outfldr,"AVG_hind_TC_small_males.Rdata"))
+
+ write_csv(AVG_hind_RKC_immature_males, file.path(outfldr,"AVG_hind_RKC_immature_males.csv"))
+ save(AVG_hind_RKC_immature_males, file=file.path(outfldr,"AVG_hind_RKC_immature_males.Rdata"))
+
+ write_csv(AVG_fut_TC_small_males, file.path(outfldr,"AVG_fut_TC_small_males.csv"))
+ save(AVG_fut_TC_small_males, file=file.path(outfldr,"AVG_fut_TC_small_males.Rdata"))
+
+ write_csv(AVG_fut_RKC_immature_males, file.path(outfldr,"AVG_fut_RKC_immature_males.csv"))
+ save(AVG_fut_RKC_immature_males, file=file.path(outfldr,"AVG_fut_RKC_immature_males.Rdata"))
+
+ write_csv(fut_TC_small_males_monthly_avg, file.path(outfldr,"fut_TC_small_males_monthly_avg.csv"))
+ save(fut_TC_small_males_monthly_avg, file=file.path(outfldr,"fut_TC_small_males_monthly_avg.Rdata"))
+
+ write_csv(fut_TC_small_males_monthly_avg, file.path(outfldr,"fut_TC_small_males_monthly_avg.csv"))
+ save(fut_TC_small_males_monthly_avg, file=file.path(outfldr,"fut_TC_small_males_monthly_avg.Rdata"))
+
+ write_csv(fut_RKC_immature_males_monthly_avg, file.path(outfldr,"fut_RKC_immature_males_monthly_avg.csv"))
+ save(fut_RKC_immature_males_monthly_avg, file=file.path(outfldr,"fut_RKC_immature_males_monthly_avg.Rdata"))
+
+
diff --git a/R/sub_scripts/make_ceattle_indices.R b/R/sub_scripts/make_ceattle_indices.R
old mode 100644
new mode 100755
index 1ee0d1b..462c881
--- a/R/sub_scripts/make_ceattle_indices.R
+++ b/R/sub_scripts/make_ceattle_indices.R
@@ -59,7 +59,7 @@
adjIN = "val_delta",
ifmissingyrs = 5,
weekIN = NULL, #"Week"
- monthIN = NULL,
+ monthIN = NULL,
GCMIN = NULL,
scenIN = NULL,
facet_rowIN = "bc", # choices=c("bc","basin","scen")
@@ -392,6 +392,8 @@
save(ceattle_vars_op, file="Data/out/CEATTLE_indices/ceattle_vars_op.Rdata")
write.csv(ceattle_vars_op, file="Data/out/CEATTLE_indices/ceattle_vars_op.csv")
+ cat("Indices saved in : ACLIM2/Data/out/CEATTLE_indices/ceattle_vars_op.Rdata")
+
# recast with vars for each column:
ceattle_vars_wide_op<- ceattle_vars_op%>%
@@ -399,17 +401,10 @@
c("units","long_name","lognorm")])))%>%
summarize_at(all_of(c("val_use")), mean, na.rm=T)%>%
tidyr::pivot_wider(names_from = "var", values_from = "val_use")
-
- # ceattle_vars_wide_op$NE_winds <- getNE_winds(vNorth=ceattle_vars_wide_op$vNorth_surface5m,
- # uEast=ceattle_vars_wide_op$uEast_surface5m)
- #
- # ceattle_vars2_op<- ceattle_vars_wide_op%>%
- # tidyr::pivot_longer(cols = c(unique(NRS_vars$var),"NE_winds"),
- # names_to = "var",
- # values_to = "val_use")%>%mutate(mn_val=val_use)
-
+
save(ceattle_vars_wide_op, file="Data/out/CEATTLE_indices/ceattle_vars_wide_op.Rdata")
write.csv(ceattle_vars_wide_op, file="Data/out/CEATTLE_indices/ceattle_vars_wide_op.csv")
+ cat("Indices saved in : ACLIM2/Data/out/CEATTLE_indices/ceattle_vars_wide_op.Rdata")
pp<- ggplot(ceattle_vars_wide_op)+
geom_line(aes(x=year,y=annual_NE_winds,color= GCM_scen,linetype = basin),
@@ -427,7 +422,7 @@
pp
- pp<- ggplot(ceattle_vars_wide_op)+
+ pp_winds<- ggplot(ceattle_vars_wide_op)+
geom_line(aes(x=year,y=Spring_NE_winds,color= GCM_scen,linetype = basin),
alpha = 0.2,show.legend = FALSE)+
geom_smooth(aes(x=year,y=Spring_NE_winds,color= GCM_scen,
@@ -441,8 +436,8 @@
legend = "")+
scale_color_discrete()+ facet_grid(scen~bc)
- pp
+ pp_winds
jpeg(filename = file.path("Data/out/CEATTLE_indices/CEATTLE_indices2_spring_op.jpg"),
width=8,height=7,units="in",res=350)
- print(pp)
+ print(pp_winds)
dev.off()
\ No newline at end of file
diff --git a/R/sub_scripts/make_hind.R b/R/sub_scripts/make_hind.R
old mode 100644
new mode 100755
diff --git a/R/sub_scripts/make_hind_C.R b/R/sub_scripts/make_hind_C.R
old mode 100644
new mode 100755
diff --git a/R/sub_scripts/make_hist.R b/R/sub_scripts/make_hist.R
old mode 100644
new mode 100755
diff --git a/R/sub_scripts/make_hist_new.R b/R/sub_scripts/make_hist_new.R
old mode 100644
new mode 100755
diff --git a/R/sub_scripts/make_hist_old.R b/R/sub_scripts/make_hist_old.R
old mode 100644
new mode 100755
diff --git a/R/sub_scripts/make_pcod_indices_Andre.R b/R/sub_scripts/make_pcod_indices_Andre.R
new file mode 100755
index 0000000..3d0c296
--- /dev/null
+++ b/R/sub_scripts/make_pcod_indices_Andre.R
@@ -0,0 +1,534 @@
+#'
+#'
+#'
+#'make_pcod_indices_Andre.R
+#'
+#'This script generates the indices for the NRS and pcod papers
+#'for Punt et al. 20220
+#'
+#' SST SBT and pH in mar-April
+#' SST and SBT for June-Aug
+#'
+#'
+
+# rm(list=ls())
+suppressMessages(source("R/make.R"))
+
+CMIPset <- c("K20P19_CMIP6","K20P19_CMIP5")
+
+# preview possible variables
+load(file = "Data/out/weekly_vars.Rdata")
+load(file = "Data/out/srvy_vars.Rdata")
+load(paste0("Data/out/K20P19_CMIP6/allEBS_means/ACLIM_annual_hind_mn.Rdata"))
+
+varall <- unique(ACLIM_annual_hind$var)
+
+# varall
+load(paste0("Data/out/K20P19_CMIP6/allEBS_means/ACLIM_annual_fut_mn.Rdata"))
+load(paste0("Data/out/K20P19_CMIP6/allEBS_means/ACLIM_surveyrep_fut_mn.Rdata"))
+
+stitchDate <- "2019-12-30" # last date of the ACLIM hindcast
+#stitchDate_op <- "2022-05-16" #last operational hindcast date # can't be mid year for these
+stitchDate_op <- "2021-12-30" #last operational hindcast date
+# scens <- c("ssp126", "ssp585")
+# GCMs <- c("miroc", "gfdl", "cesm" )
+
+# Now compile the Pcod indices:
+#--------------------------------------
+grpby <- c("type","var","basin",
+ "year","sim","gcmcmip","GCM","scen","sim_type",
+ "bc","GCM_scen","GCM_scen_sim", "CMIP" )
+
+sumat <- c("jday","mnDate","val_use","mnVal_hind",
+ "val_delta","val_biascorrected","val_raw")
+
+# make PCOD_indices.csv using the ACLIM hindcast only as well as
+# PCOD_indices_op.csv, the operational hindcast filled in for 2019-2022
+varlist <- c("pH_bottom5m","pH_depthavg","pH_integrated","pH_surface5m")
+varlist <- c("temp_surface5m", "temp_bottom5m")
+i <- 0
+v <- varlist[1]
+mIN <- 3:4
+for(v in varlist){
+
+ i <- i + 1
+ cat("compiling indices : ",v,"\n")
+
+ # get the variable you want:
+ df <- get_var( typeIN = "monthly",
+ monthIN = mIN,
+ plotbasin = c("SEBS"),
+ plotvar = v,
+ bcIN = "bias corrected",
+ CMIPIN = CMIPset,
+ plothist = T, # ignore the hist runs
+ removeyr1 = T) # "Remove first year of proj ( burn in)")
+
+ tmpda <- df$dat%>%
+ group_by(across(all_of(grpby)))%>%
+ summarize_at(all_of(sumat), mean, na.rm=T)%>%
+ mutate(mn_val = val_use, var = paste0("Mar2Apr_",var))%>%ungroup()
+
+ tmpd <- stitchTS(dat = tmpda, maxD = stitchDate)
+ tmpd <- tmpd%>%mutate(type = "Pcod indices")
+ ggplot(tmpd)+
+ geom_line(aes(x=year, y=mn_val,color=GCM_scen_sim))+
+ facet_grid(.~scen)
+
+ # now for operational hindcasts:
+ dfop <- get_var_ophind( typeIN = "monthly",
+ monthIN = mIN,
+ stitchDateIN = stitchDate,
+ plotbasin = c("SEBS"),
+ plotvar = v,
+ bcIN = "bias corrected",
+ CMIPIN = CMIPset,
+ jday_rangeIN = c(0,365),
+ plothist = T, # ignore the hist runs
+ removeyr1 = T,
+ adjIN = "val_delta",
+ ifmissingyrs = 5,
+ weekIN = NULL, #"Week"
+ SeasonIN = NULL,
+ GCMIN = NULL,
+ scenIN = NULL,
+ facet_rowIN = "bc", # choices=c("bc","basin","scen")
+ facet_colIN = "scen")
+# dfop$dat%>%filter(var==v,year%in%(2022))%>%select(year,mo, basin, var,val_use,val_delta,val_raw)
+ tmpdop <- dfop$dat%>%
+ group_by(across(all_of(c(grpby,"GCM2","GCM2_scen_sim"))))%>%
+ summarize_at(all_of(sumat), mean, na.rm=T)%>%
+ mutate(mn_val = val_use, var = paste0("Mar2Apr_",var))%>%ungroup()
+
+ # tmpdop%>%filter(var==v,year%in%(2022))%>%select(year, basin, var,val_use,val_delta,val_raw)
+ tmpdop <- stitchTS(dat = tmpdop, maxD = stitchDate_op)
+ tmpdop <- tmpdop%>%mutate(type = "Pcod indices")
+ tmpdop%>%filter(var==v,year%in%(2022))%>%select(year, basin, var,val_use,val_delta,val_raw,GCM_scen,sim)
+
+ ggplot(tmpdop)+
+ geom_line(aes(x=year, y=mn_val,color=GCM_scen_sim))+
+ facet_grid(.~scen)
+
+ if(i==1){
+ PCOD_vars <- tmpd
+ PCOD_vars_op <- tmpdop
+ }else{
+ PCOD_vars <- rbind(PCOD_vars,tmpd)
+ PCOD_vars_op <- rbind(PCOD_vars_op,tmpdop)
+ }
+
+ rm(df)
+ rm(dfop)
+ rm(tmpd)
+ rm(tmpdop)
+
+}
+
+varlist <- c("pH_bottom5m","pH_depthavg","pH_integrated","pH_surface5m")
+i <- 0
+v <- varlist[1]
+mIN <- 3:4
+
+for(v in varlist){
+
+ cat("compiling indices : ",v,"\n")
+
+ # get the variable you want:
+ df <- get_var( typeIN = "monthly",
+ monthIN = mIN,
+ plotbasin = c("SEBS"),
+ plotvar = v,
+ bcIN = "bias corrected",
+ CMIPIN = CMIPset,
+ plothist = T, # ignore the hist runs
+ removeyr1 = T) # "Remove first year of proj ( burn in)")
+
+ tmpda <- df$dat%>%
+ group_by(across(all_of(grpby)))%>%
+ summarize_at(all_of(sumat), mean, na.rm=T)%>%
+ mutate(mn_val = val_use, var = paste0("Mar2Apr_",var))%>%ungroup()
+
+ tmpd <- stitchTS(dat = tmpda, maxD = stitchDate)
+ tmpd <- tmpd%>%mutate(type = "Pcod indices")
+ ggplot(tmpd)+
+ geom_line(aes(x=year, y=mn_val,color=GCM_scen_sim))+
+ facet_grid(.~scen)
+
+
+ PCOD_vars <- rbind(PCOD_vars,tmpd)
+
+
+ rm(df)
+ rm(tmpd)
+
+}
+
+
+
+varlist <- c("temp_surface5m", "temp_bottom5m")
+i <- 0
+v <- varlist[1]
+mIN <- 6:9
+for(v in varlist){
+ cat("compiling indices : ",v,"\n")
+
+ # get the variable you want:
+ df <- get_var( typeIN = "monthly",
+ monthIN = mIN,
+ plotbasin = c("SEBS"),
+ plotvar = v,
+ bcIN = "bias corrected",
+ CMIPIN = CMIPset,
+ plothist = T, # ignore the hist runs
+ removeyr1 = T) # "Remove first year of proj ( burn in)")
+
+ tmpda <- df$dat%>%
+ group_by(across(all_of(grpby)))%>%
+ summarize_at(all_of(sumat), mean, na.rm=T)%>%
+ mutate(mn_val=val_use, var = paste0("Jun2Aug_",var))%>%ungroup()
+
+ tmpd <- stitchTS(dat = tmpda, maxD = stitchDate)
+ tmpd <- tmpd%>%mutate(type = "Pcod indices")
+ ggplot(tmpd)+
+ geom_line(aes(x=year, y=mn_val,color=GCM_scen_sim))+
+ facet_grid(.~scen)
+
+ # now for operational hindcasts:
+ dfop <- get_var_ophind( typeIN = "monthly",
+ monthIN = mIN,
+ stitchDateIN = stitchDate,
+ plotbasin = c("SEBS"),
+ plotvar = v,
+ bcIN = "bias corrected",
+ CMIPIN = CMIPset,
+ jday_rangeIN = c(0,365),
+ plothist = T, # ignore the hist runs
+ removeyr1 = T,
+ adjIN = "val_delta",
+ ifmissingyrs = 5,
+ weekIN = NULL, #"Week"
+ SeasonIN = NULL,
+ GCMIN = NULL,
+ scenIN = NULL,
+ facet_rowIN = "bc", # choices=c("bc","basin","scen")
+ facet_colIN = "scen")
+ #dfop$dat%>%filter(var==v,year%in%(2022))%>%select(year,mo, basin, var,val_use,val_delta,val_raw)
+ tmpdop <- dfop$dat%>%
+ group_by(across(all_of(c(grpby,"GCM2","GCM2_scen_sim"))))%>%
+ summarize_at(all_of(sumat), mean, na.rm=T)%>%
+ mutate(mn_val=val_use, var = paste0("Jun2Aug_",var))%>%ungroup()
+# tmpdop%>%filter(var==v,year%in%(2022))%>%select(year, basin, var,val_use,val_delta,val_raw)
+ tmpdop <- stitchTS(dat = tmpdop, maxD = stitchDate_op)
+ tmpdop <- tmpdop%>%mutate(type = "Pcod indices")
+ tmpdop%>%filter(var==v,year%in%(2022))%>%select(year, basin, var,val_use,val_delta,val_raw,GCM_scen,sim)
+
+ ggplot(tmpdop)+
+ geom_line(aes(x=year, y=mn_val,color=GCM_scen_sim))+
+ facet_grid(.~scen)
+
+ PCOD_vars <- rbind(PCOD_vars,tmpd)
+ PCOD_vars_op <- rbind(PCOD_vars_op,tmpdop)
+
+ rm(df)
+ rm(dfop)
+ rm(tmpd)
+ rm(tmpdop)
+
+}
+
+if(any(unique(PCOD_vars_op$sim)=="ACLIMregion_B10K-K20P19_CORECFS"))
+ stop(paste("ACLIMregion_B10K-K20P19_CORECFS in operational hindcast!!!\n ",paste(varlist,"\n")))
+
+
+PCOD_vars <- PCOD_vars%>%ungroup()
+PCOD_vars_op <- PCOD_vars_op%>%ungroup()
+
+
+#Annual Indices
+# annual values
+varlist <- c("temp_surface5m", "temp_bottom5m")
+
+#
+# for(v in varlist){
+# cat("compiling indices : ",v,"\n")
+# # get the variable you want:
+# df <- get_var( typeIN = "annual",
+# plotbasin = c("SEBS"),
+# plotvar = v,
+# bcIN = "bias corrected",
+# CMIPIN = CMIPset,
+# plothist = T, # ignore the hist runs
+# removeyr1 = T) # "Remove first year of projection ( burn in)")
+# tmpd <- df$dat%>%mutate(season="annual")%>%
+# group_by(across(all_of(grpby)))%>%
+# summarize_at(all_of(sumat), mean, na.rm=T)%>%
+# mutate(mn_val=val_use, var = paste0("annual_",var))
+#
+# tmpd <- stitchTS(dat = tmpd, stitchDate)
+# tmpd <- tmpd%>%mutate(type = "Pcod indices")
+# PCOD_vars <- rbind(PCOD_vars,tmpd)
+# rm(df)
+# rm(tmpd)
+#
+# # now for operational hindcasts:
+# dfop <- get_var_ophind(stitchDateIN = stitchDate,
+# typeIN = "annual",
+# plotbasin = c("SEBS"),
+# plotvar = v,
+# bcIN = "bias corrected",
+# CMIPIN = CMIPset,
+# plothist = T, # ignore the hist runs
+# removeyr1 = T) # "Remove first year of proj ( burn in)")
+#
+# tmpdop <- dfop$dat%>%mutate(season="annual")%>%
+# group_by(across(all_of(c(grpby,"GCM2","GCM2_scen_sim"))))%>%
+# summarize_at(all_of(sumat), mean, na.rm=T)%>%
+# mutate(mn_val=val_use, var = paste0("annual_",var))
+#
+# tmpdop <- stitchTS(dat = tmpdop, stitchDate_op)
+# tmpdop <- tmpdop%>%mutate(type = "Pcod indices")
+# PCOD_vars_op <- rbind(PCOD_vars_op,tmpdop)
+# rm(dfop)
+# rm(tmpdop)
+# }
+#
+# if(any(unique(PCOD_vars_op$sim)=="ACLIMregion_B10K-K20P19_CORECFS"))
+# stop(paste("ACLIMregion_B10K-K20P19_CORECFS in operational hindcast!!!\n ",paste(varlist,"\n")))
+#
+# varlist <- c("pH_bottom5m","pH_depthavg","pH_integrated","pH_surface5m")
+#
+# for(v in varlist){
+# cat("compiling indices : ",v,"\n")
+# # get the variable you want:
+# df <- get_var( typeIN = "annual",
+# plotbasin = c("SEBS"),
+# plotvar = v,
+# bcIN = "bias corrected",
+# CMIPIN = CMIPset,
+# plothist = T, # ignore the hist runs
+# removeyr1 = T) # "Remove first year of projection ( burn in)")
+# tmpd <- df$dat%>%mutate(season="annual")%>%
+# group_by(across(all_of(grpby)))%>%
+# summarize_at(all_of(sumat), mean, na.rm=T)%>%
+# mutate(mn_val=val_use, var = paste0("annual_",var))
+#
+# tmpd <- stitchTS(dat = tmpd, stitchDate)
+# tmpd <- tmpd%>%mutate(type = "Pcod indices")
+# PCOD_vars <- rbind(PCOD_vars,tmpd)
+# rm(df)
+# rm(tmpd)
+#
+# }
+#
+# if(any(unique(PCOD_vars_op$sim)=="ACLIMregion_B10K-K20P19_CORECFS"))
+# stop(paste("ACLIMregion_B10K-K20P19_CORECFS in operational hindcast!!!\n ",paste(varlist,"\n")))
+#
+
+# fix the caps issue:
+PCOD_vars<-PCOD_vars%>%
+ mutate(GCM = gsub("MIROC","miroc",GCM))%>%
+ mutate(GCM = gsub("GFDL","gfdl",GCM))%>%
+ mutate(GCM = gsub("CESM","cesm",GCM))%>%
+ mutate(GCM_scen = gsub("MIROC","miroc",GCM_scen))%>%
+ mutate(GCM_scen = gsub("GFDL","gfdl",GCM_scen))%>%
+ mutate(GCM_scen = gsub("CESM","cesm",GCM_scen))%>%
+ mutate(GCM_scen_sim = gsub("MIROC","miroc",GCM_scen_sim))%>%
+ mutate(GCM_scen_sim = gsub("GFDL","gfdl",GCM_scen_sim))%>%
+ mutate(GCM_scen_sim = gsub("CESM","cesm",GCM_scen_sim))%>%
+ mutate(GCM = factor(GCM, levels =c("hind","gfdl","cesm","miroc")))
+
+# fix the caps issue:
+PCOD_vars_op<-PCOD_vars_op%>%
+ mutate(GCM = gsub("MIROC","miroc",GCM))%>%
+ mutate(GCM = gsub("GFDL","gfdl",GCM))%>%
+ mutate(GCM = gsub("CESM","cesm",GCM))%>%
+ mutate(GCM_scen = gsub("MIROC","miroc",GCM_scen))%>%
+ mutate(GCM_scen = gsub("GFDL","gfdl",GCM_scen))%>%
+ mutate(GCM_scen = gsub("CESM","cesm",GCM_scen))%>%
+ mutate(GCM_scen_sim = gsub("MIROC","miroc",GCM_scen_sim))%>%
+ mutate(GCM_scen_sim = gsub("GFDL","gfdl",GCM_scen_sim))%>%
+ mutate(GCM_scen_sim = gsub("CESM","cesm",GCM_scen_sim))%>%
+ mutate(GCM = factor(GCM, levels =c("hind","gfdl","cesm","miroc")))
+
+# save indices
+if(!dir.exists("Data/out/PCOD_indices"))
+ dir.create("Data/out/PCOD_indices")
+
+save(PCOD_vars, file="Data/out/PCOD_indices/PCOD_vars.Rdata")
+write.csv(PCOD_vars, file="Data/out/PCOD_indices/PCOD_vars.csv")
+save(PCOD_vars_op, file="Data/out/PCOD_indices/PCOD_vars_op.Rdata")
+write.csv(PCOD_vars_op, file="Data/out/PCOD_indices/PCOD_vars_op.csv")
+
+# recast with vars for each column:
+PCOD_vars_wide<- PCOD_vars%>%
+ group_by(across(all_of(grpby[!grpby%in%
+ c("units","long_name","lognorm")])))%>%
+ summarize_at(all_of(c("val_use")), mean, na.rm=T)%>%
+ tidyr::pivot_wider(names_from = "var", values_from = "val_use")
+
+# reorder to put annual indices at the end
+cc <- colnames(PCOD_vars_wide)
+cc <- cc[-grep("annual",cc)]
+PCOD_vars_wide<- PCOD_vars_wide%>%relocate(all_of(cc))
+
+PCOD_vars2<- PCOD_vars_wide%>%
+ tidyr::pivot_longer(cols = c(unique(PCOD_vars$var)),
+ names_to = "var",
+ values_to = "val_use")%>%mutate(mn_val=val_use)%>%ungroup()
+
+# recast with vars for each column:
+PCOD_vars_wide_op<- PCOD_vars_op%>%
+ group_by(across(all_of(grpby[!grpby%in%
+ c("units","long_name","lognorm")])))%>%
+ summarize_at(all_of(c("val_use")), mean, na.rm=T)%>%
+ tidyr::pivot_wider(names_from = "var", values_from = "val_use")
+
+# PCOD_vars_wide_op$NE_winds <-
+# getNE_winds(vNorth=PCOD_vars_wide_op$vNorth_surface5m,
+# uEast=PCOD_vars_wide_op$uEast_surface5m)
+
+# reorder to put annual indices at the end
+cc <- colnames(PCOD_vars_wide_op)
+cc <- cc[-grep("annual",cc)]
+PCOD_vars_wide_op<- PCOD_vars_wide_op%>%
+ relocate(all_of(cc))
+
+
+PCOD_vars2_op<- PCOD_vars_wide_op%>%
+ tidyr::pivot_longer(cols = c(unique(PCOD_vars_op$var)),
+ names_to = "var",
+ values_to = "val_use")%>%
+ mutate(mn_val=val_use)%>%ungroup()
+
+save(PCOD_vars_wide, file="Data/out/PCOD_indices/PCOD_vars_wide.Rdata")
+write.csv(PCOD_vars_wide, file="Data/out/PCOD_indices/PCOD_vars_wide.csv")
+save(PCOD_vars_wide_op, file="Data/out/PCOD_indices/PCOD_vars_wide_op.Rdata")
+write.csv(PCOD_vars_wide_op, file="Data/out/PCOD_indices/PCOD_vars_wide_op.csv")
+
+
+tt <- PCOD_vars_wide_op%>%filter(scen=="rcp45",year ==1976)
+unique(tt$sim)
+
+tt2 <- PCOD_vars_op%>%filter(scen=="rcp45",year ==1976,var=="fracbelow1")
+unique(tt2$sim)
+
+tt1 <- PCOD_vars_wide%>%filter(scen=="rcp45",year ==1976)
+unique(tt1$sim)
+
+tt12 <- PCOD_vars%>%filter(scen=="rcp45",year ==1976,var=="fracbelow1")
+unique(tt12$sim)
+
+
+sclr <- 2
+jpeg(filename = file.path("Data/out/PCOD_indices/PCOD_vars.jpg"),
+ width=6*sclr, height=9*sclr, units="in", res = 350)
+print(plotTS(PCOD_vars )+labs(title="Pcod indices, with ACLIM hind"))
+dev.off()
+
+
+
+jpeg(filename = file.path("Data/out/PCOD_indices/PCOD_vars_byGCM.jpg"),
+ width=6*sclr, height=9*sclr, units="in", res = 350)
+print(plotTS(PCOD_vars )+
+ labs(title="Pcod indices")+facet_grid(var~GCM,scales='free_y'))
+dev.off()
+
+
+jpeg(filename = file.path("Data/out/PCOD_indices/PCOD_vars_op_byGCM.jpg"),
+ width=6*sclr, height=9*sclr, units="in", res = 350)
+print(plotTS(PCOD_vars_op )+
+ labs(title="Pcod indices, with operational hind")+facet_grid(var~GCM,scales='free_y'))
+dev.off()
+
+jpeg(filename = file.path("Data/out/PCOD_indices/PCOD_vars_op.jpg"),
+ width=6*sclr, height=9*sclr, units="in", res = 350)
+print(plotTS(PCOD_vars_op )+labs(title="Pcod indices, with operational hind"))
+dev.off()
+
+
+plotTS(PCOD_vars_op%>%filter(var=="Mar2Apr_temp_bottom5m"), plotvalIN = "val_raw" )+
+ coord_cartesian(ylim=c(7.5,8.1))
+
+
+pp2 <- plotTS(PCOD_vars2%>%mutate(units="",mnDate=year) )+
+ facet_wrap(.~var,scales="free")+ylab("val")+
+ labs(title = "ACLIM hindcast +projections")
+jpeg(filename = file.path("Data/out/PCOD_indices/PCOD_indices.jpg"),
+ width=10,height=8,units="in",res=350)
+print(pp2)
+dev.off()
+
+
+pp2 <- plotTS(PCOD_vars2_op%>%mutate(units="",mnDate=year) )+
+ facet_wrap(.~var,scales="free")+ylab("val")+
+ labs(title = "operational hindcast + projections")
+jpeg(filename = file.path("Data/out/PCOD_indices/PCOD_indices_op.jpg"),
+ width=10,height=8,units="in",res=350)
+print(pp2)
+dev.off()
+
+#
+# pp<- ggplot(PCOD_vars_wide)+
+# geom_line(aes(x=year,y=NE_winds,color= GCM_scen,linetype = basin),
+# alpha = 0.2,show.legend = FALSE)+
+# geom_smooth(aes(x=year,y=NE_winds,color= GCM_scen,
+# fill=GCM_scen,linetype = basin),
+# alpha=0.1,method="loess",formula='y ~ x',span = .5)+
+# theme_minimal() +
+# labs(x="Year",
+# y="NE_winds (m s^-1)",
+# subtitle = "",
+# title = "NE_winds, ACLIM hindcast",
+# legend = "")+
+# scale_color_discrete()+ facet_grid(scen~bc)
+#
+# pp
+#
+# jpeg(filename = file.path("Data/out/PCOD_indices/PCOD_indices2.jpg"),
+# width=8,height=7,units="in",res=350)
+# print(pp)
+# dev.off()
+
+#
+# pp<- ggplot(PCOD_vars_wide_op)+
+# geom_line(aes(x=year,y=NE_winds,color= GCM_scen,linetype = basin),
+# alpha = 0.2,show.legend = FALSE)+
+# geom_smooth(aes(x=year,y=NE_winds,color= GCM_scen,
+# fill=GCM_scen,linetype = basin),
+# alpha=0.1,method="loess",formula='y ~ x',span = .5)+
+# theme_minimal() +
+# labs(x="Year",
+# y="NE_winds (m s^-1)",
+# subtitle = "",
+# title = "NE_winds, operational hindcast",
+# legend = "")+
+# scale_color_discrete()+ facet_grid(scen~bc)
+#
+# pp
+#
+# jpeg(filename = file.path("Data/out/PCOD_indices/PCOD_indices2_op.jpg"),
+# width=8,height=7,units="in",res=350)
+# print(pp)
+# dev.off()
+
+
+pp<- ggplot(PCOD_vars_wide_op)+
+ geom_hline(yintercept = 0,color="gray")+
+ geom_line(aes(x=year,y=Mar2Apr_temp_bottom5m,color= GCM_scen,linetype = basin),
+ alpha = 0.2,show.legend = FALSE)+
+ geom_smooth(aes(x=year,y=Mar2Apr_temp_bottom5m,color= GCM_scen,
+ fill=GCM_scen,linetype = basin),
+ alpha=0.1,method="loess",formula='y ~ x',span = .5)+
+ theme_minimal() +
+ labs(x="Year",
+ subtitle = "",
+ title = "Mar2Apr_temp_bottom5m, operational hindcast",
+ legend = "")+
+ scale_color_discrete()+ facet_grid(scen~bc)
+
+pp
+jpeg(filename = file.path("Data/out/PCOD_indices/PCOD_indices2_op_Mar2Apr_temp_bottom5m.jpg"),
+ width=8,height=7,units="in",res=350)
+print(pp)
+dev.off()
+
diff --git a/R/sub_scripts/make_proj.R b/R/sub_scripts/make_proj.R
old mode 100644
new mode 100755
diff --git a/R/sub_scripts/make_proj_C.R b/R/sub_scripts/make_proj_C.R
old mode 100644
new mode 100755
diff --git a/R/sub_scripts/make_proj_new.R b/R/sub_scripts/make_proj_new.R
old mode 100644
new mode 100755
diff --git a/R/sub_scripts/make_salmon_indices.R b/R/sub_scripts/make_salmon_indices.R
new file mode 100644
index 0000000..0682680
--- /dev/null
+++ b/R/sub_scripts/make_salmon_indices.R
@@ -0,0 +1,417 @@
+#'
+#'
+#' make_salmon_indices.R
+#' K. Holsman
+#' 2024
+#' generate indices for E. Yasumiishi
+#'
+#'
+
+#
+# 1. temperature_surface5m for months 9:10 and strata 70 & 71
+# 2. NW ditrrection for months 5:6 and strata 71
+# 3. temperature_surface5m for months 7:10 and strata 90,61,62
+
+# load ACLIM packages and functions
+rm(list=ls())
+suppressMessages(source("R/make.R"))
+
+cat("\n\nsetting things up\n")
+
+# set up paths
+outfldr <- "Data/out/salmon_indices"
+if(!dir.exists(outfldr)) dir.create(outfldr)
+
+strata_set <- list(c(70,71),c(71),c(90,61,62))
+mo_set <- list(c(9,10),c(5,6),7:10)
+var_set <- list(c("temp_surface5m"),c("NEwinds"),"temp_surface5m")
+names(var_set) <- c("fall_SST_70_71","Winds_71","7through10_SST_90_61_62")
+
+NEWinds <- c(0,1,0)
+
+# PCOD_vars_wide_op$NE_winds <-
+# getNE_winds(vNorth=PCOD_vars_wide_op$vNorth_surface5m,
+# uEast=PCOD_vars_wide_op$uEast_surface5m)
+
+#varlist <- c("pH_bottom5m","pH_depthavg","pH_integrated","pH_surface5m","temp_surface5m", "temp_bottom5m")
+CMIPset <- c("K20P19_CMIP6","K20P19_CMIP5")
+stitchDate_op <- "2021-12-30" #last operational hindcast date
+
+# preview possible variables
+load(file = "Data/out/weekly_vars.Rdata")
+load(file = "Data/out/srvy_vars.Rdata")
+# save setup for plotting
+save(list=ls(),file=file.path(outfldr,"salmon_indices_setup.Rdata"))
+
+
+# ------------------------------------
+# Hindcast First
+# ------------------------------------
+cat("running hindcast summary\n")
+
+load("Data/out/K20P19_CMIP6/BC_ACLIMregion/ACLIMregion_B10K-K20P19_CORECFS_BC_hind.Rdata")
+rm_cols <- c( "lognorm","type", "val_raw","mn_val" ,"sd_val", "n_val",
+ "sdVal_hind_mo", "sdVal_hind_yr")
+
+select_list <- c("sim","year","var_salmon","basin","strata","strata_area_km2","units",
+ "season","mo","wk","jday","mnDate","qry_date","sim_type","mday",
+ "val_use","mnVal_hind","sdVal_hind","nVal_hind","seVal_hind","val_raw","var")
+
+for(i in 1:length(strata_set)){
+
+ if(var_set[[i]]=="NEwinds"){
+ my_vars2 <- function(IN) select_list[select_list%in%names(IN)]
+ valset <- c("val_use","mnVal_hind","sdVal_hind","nVal_hind", "seVal_hind","val_raw")
+
+ hindDat_weekly_strata_tmp <- hind%>%
+ filter(strata%in%strata_set[[i]],
+ mo%in%mo_set[[i]],
+ var%in%c("vNorth_surface5m","uEast_surface5m"))%>%
+ mutate(val_use = val_raw, units = "meter second-1")%>%
+ rowwise()%>%mutate(mday=strptime(mnDate,format = "%Y-%m-%d")$mday)%>%ungroup()%>%
+ mutate(var_salmon = names(var_set)[i])
+
+ # remove all the extra stuff:
+ hindDat_weekly_strata_tmp <- hindDat_weekly_strata_tmp%>%select(my_vars2(hindDat_weekly_strata_tmp))
+
+ hindDat_weekly_strata_tmp <- hindDat_weekly_strata_tmp%>%
+ pivot_wider(names_from = var,
+ values_from = all_of(valset))%>%
+ ungroup()%>%data.frame()
+
+ for(k in 1:length(valset)){
+
+ eval(parse(text = paste0(" hindDat_weekly_strata_tmp <- hindDat_weekly_strata_tmp%>%
+ rowwise()%>%mutate(",valset[k], "= getNE_winds
+ (vNorth = ",valset[k],"_vNorth_surface5m,
+ uEast = ",valset[k],"_uEast_surface5m))%>%
+ select(-",valset[k],"_vNorth_surface5m,-",
+ valset[k],"_uEast_surface5m)%>%
+ ungroup()%>%data.frame()") ))
+ }
+ hindDat_weekly_strata_tmp <- hindDat_weekly_strata_tmp%>%mutate(var = "NEwinds")
+
+
+ }else{
+ my_vars <- function(IN) select_list[select_list%in%names(IN)]
+
+ hindDat_weekly_strata_tmp <- hind%>%
+ filter(strata%in%strata_set[[i]],
+ mo%in%mo_set[[i]],
+ var%in%var_set[[i]])%>%
+ mutate(val_use = val_raw)%>%
+ rowwise()%>%mutate(mday=strptime(mnDate,format = "%Y-%m-%d")$mday)%>%ungroup()%>%
+ mutate(var_salmon = names(var_set)[i])
+
+ # remove all the extra stuff:
+ hindDat_weekly_strata_tmp <- hindDat_weekly_strata_tmp%>%select(my_vars(hindDat_weekly_strata_tmp))
+
+
+ }
+
+ grp_names <- names(hindDat_weekly_strata_tmp)
+ nms3 <- grp_names[!grp_names%in%rm_cols]
+
+ hindDat_weekly_strata_tmp <-hindDat_weekly_strata_tmp%>%
+ select(!!!syms(nms3))%>%
+ relocate(val_use, .before = mnVal_hind)%>%
+ relocate(var_salmon, .before = var)%>%
+ ungroup()%>%data.frame()
+
+ grp_names <- names(hindDat_weekly_strata_tmp)[!names(hindDat_weekly_strata_tmp)%in%c("strata","strata_area_km2","wk","season","basin","mday")]
+ nms <- c("val_use","jday","mnDate",
+ "mnVal_hind","sdVal_hind","nVal_hind",
+ "seVal_hind")
+
+ nms2 <- grp_names[!grp_names%in%nms]
+ nms3 <- grp_names[!grp_names%in%rm_cols]
+
+ hindDat_monthly_tmp <- hindDat_weekly_strata_tmp%>%
+ select(!!!syms(grp_names))%>%group_by(!!!syms(nms2))%>%
+ summarize_at(nms,mean,na.rm=T)%>%
+ select(!!!syms(nms3))%>%
+ relocate(val_use, .before = mnVal_hind)%>%
+ relocate(var_salmon, .before = var)%>%
+ ungroup()%>%data.frame()
+
+
+ grp_names <- names(hindDat_weekly_strata_tmp)[!names(hindDat_weekly_strata_tmp)%in%c("strata","lognorm","type",
+ "strata_area_km2",
+ "mo","wk","season","basin","mday")]
+ nms <- c("val_use","jday","mnDate",
+ "mnVal_hind","sdVal_hind","nVal_hind",
+ "seVal_hind")
+
+ nms2 <- grp_names[!grp_names%in%nms]
+ nms3 <- grp_names[!grp_names%in%rm_cols]
+
+ hindDat_avg_tmp <- hindDat_weekly_strata_tmp%>%
+ select(!!!syms(grp_names))%>%group_by(!!!syms(nms2))%>%
+ summarize_at(nms,mean,na.rm=T)%>%
+ select(!!!syms(nms3))%>%
+ relocate(val_use, .before = mnVal_hind)%>%
+ relocate(var_salmon, .before = var)%>%
+ ungroup()%>%data.frame()
+
+
+ if(i==1){
+ hindDat_avg <- hindDat_avg_tmp
+ hindDat_monthly <- hindDat_monthly_tmp
+ hindDat_weekly_strata <- hindDat_weekly_strata_tmp
+
+ }else{
+ hindDat_avg <- rbind(hindDat_avg,hindDat_avg_tmp)
+ hindDat_monthly <- rbind(hindDat_monthly,hindDat_monthly_tmp)
+ hindDat_weekly_strata <- rbind(hindDat_weekly_strata,hindDat_weekly_strata_tmp)
+ }
+ rm(hindDat_avg_tmp)
+ rm(hindDat_monthly_tmp)
+ rm(hindDat_weekly_strata_tmp)
+
+}
+
+
+
+ # output the data as Rdata and CSV
+ write_csv(hindDat_weekly_strata, file.path(outfldr,"salmon_hindDat_weekly_strata.csv"))
+ save(hindDat_weekly_strata, file=file.path(outfldr,"salmon_hindDat_weekly_stratat.Rdata"))
+
+ write_csv(hindDat_monthly, file.path(outfldr,"salmon_hindDat_monthly.csv"))
+ save(hindDat_monthly, file=file.path(outfldr,"salmon_hindDat_monthly.Rdata"))
+
+ write_csv(hindDat_avg, file.path(outfldr,"salmon_hindDat.csv"))
+ save(hindDat_avg, file=file.path(outfldr,"salmon_hindDat.Rdata"))
+
+
+ #rm(hindDat_weekly_strata)
+ rm(hindDat_monthly)
+ rm(hindDat_avg)
+
+
+# ------------------------------------
+# Projections Next
+# ------------------------------------
+select_list <- unique(c(names(hindDat_weekly_strata),"long_name","qry_date","station_id",
+ "sim_type","GCM","RCP","mod", "CMIP", "mnVal_hind","sdVal_hind","val_delta","val_raw"))
+ my_vars <- function(IN) select_list[select_list%in%names(IN)]
+
+ jj <- 0
+
+CMIP <- "CMIP6"
+ii <- 1
+i <- 1
+
+for (CMIP in c("CMIP6","CMIP5")){
+
+ if(CMIP == "CMIP6"){
+
+ fl_base_srv <- "ACLIMsurveyrep_B10K-K20P19_CMIP6_"
+ fl_path_srv <- "Data/out/K20P19_CMIP6/BC_ACLIMsurveyrep"
+
+ fl_base_reg <- "ACLIMregion_B10K-K20P19_CMIP6_"
+ fl_path_reg <- "Data/out/K20P19_CMIP6/BC_ACLIMregion"
+
+
+ # CMIP6
+ ssps <- c("ssp126","ssp585")
+ gcms <- c("cesm","miroc","gfdl")
+ set <- expand.grid(gcms,ssps,stringsAsFactors = F)
+
+ }
+
+ if(CMIP == "CMIP5"){
+ #now get monthly
+ fl_base_srv <- "ACLIMsurveyrep_B10K-K20P19_CMIP5_"
+ fl_path_srv <- "Data/out/K20P19_CMIP5/BC_ACLIMsurveyrep"
+
+ fl_base_reg <- "ACLIMregion_B10K-K20P19_CMIP5_"
+ fl_path_reg <- "Data/out/K20P19_CMIP5/BC_ACLIMregion"
+
+ # now for CMIP5
+ gcms <- c("CESM","MIROC","GFDL")
+ ssps <- c("rcp45","rcp85")
+ set <- expand.grid(gcms,ssps,stringsAsFactors = F)
+
+ }
+
+
+ for(i in 1:dim(set)[1]){
+ fl <- paste0(fl_base_reg,set[i,1],"_",set[i,2],"_BC_fut.Rdata")
+ load(file.path(fl_path_reg,fl))
+
+ for(ii in 1:length(strata_set)){
+ cat(paste("running proj summary for",CMIP,": strata_set",names(var_set)[ii],set[i,1],set[i,2],"\n"))
+
+
+ fut$mo <- as.numeric(substr(fut$mnDate,6,7))
+ fut$mday <- as.numeric(substr(fut$mnDate,9,10))
+
+ if(var_set[[ii]]=="NEwinds"){
+
+ valset <- c("val_use","mnVal_hind","sdVal_hind","val_raw","val_delta")
+ # PCOD_vars_wide_op$NE_winds <-
+ # getNE_winds(vNorth=PCOD_vars_wide_op$vNorth_surface5m,
+ # uEast=PCOD_vars_wide_op$uEast_surface5m)
+
+ futDat_weekly_strata_tmp <- fut%>%
+ filter(strata%in%strata_set[[ii]],
+ mo%in%mo_set[[ii]],
+ var%in%c("vNorth_surface5m","uEast_surface5m"))%>%
+ mutate(val_use = val_delta,units = "meter second-1" )%>%
+ mutate(var_salmon = names(var_set)[ii])%>%
+ ungroup()%>%data.frame()
+
+ # remove all the extra stuff:
+ futDat_weekly_strata_tmp <- futDat_weekly_strata_tmp%>%select(my_vars(futDat_weekly_strata_tmp))
+
+ futDat_weekly_strata_tmp <- futDat_weekly_strata_tmp%>%
+ pivot_wider(names_from = var,
+ values_from = all_of(valset))%>%
+ ungroup()%>%data.frame()
+
+ for(k in 1:length(valset)){
+
+ eval(parse(text = paste0(" futDat_weekly_strata_tmp <- futDat_weekly_strata_tmp%>%
+ rowwise()%>%mutate(",valset[k], "= getNE_winds
+ (vNorth = ",valset[k],"_vNorth_surface5m,
+ uEast = ",valset[k],"_uEast_surface5m))%>%
+ select(-",valset[k],"_vNorth_surface5m,-",
+ valset[k],"_uEast_surface5m)%>%
+ ungroup()%>%data.frame()") ))
+ }
+ futDat_weekly_strata_tmp <- futDat_weekly_strata_tmp%>%mutate(var = "NEwinds")
+
+
+ }else{
+
+ futDat_weekly_strata_tmp <- fut%>%
+ filter(strata%in%strata_set[[ii]],
+ mo%in%mo_set[[ii]],
+ var%in%var_set[[ii]])%>%
+ mutate(val_use = val_delta)%>%
+ mutate(var_salmon = names(var_set)[ii])%>%
+ ungroup()%>%data.frame()
+
+ # remove all the extra stuff:
+ futDat_weekly_strata_tmp <- futDat_weekly_strata_tmp%>%select(my_vars(futDat_weekly_strata_tmp))
+
+ }
+
+
+ grp_names <- names(futDat_weekly_strata_tmp)[!names(futDat_weekly_strata_tmp)%in%c("strata","strata_area_km2","wk","season","basin","mday")]
+ nms <- c("val_use","jday","mnDate",
+ "mnVal_hind","sdVal_hind","val_raw","val_delta")
+
+ nms2 <- grp_names[!grp_names%in%nms]
+
+ futDat_monthly_tmp <- futDat_weekly_strata_tmp%>%
+ select(!!!syms(grp_names))%>%group_by(!!!syms(nms2))%>%
+ summarize_at(nms,mean,na.rm=T)%>%
+ relocate(val_use, .before = mnVal_hind)%>%
+ relocate(var_salmon, .before = var)%>%
+ ungroup()%>%data.frame()
+
+
+ grp_names <- names(futDat_weekly_strata_tmp)[!names(futDat_weekly_strata_tmp)%in%c("strata","lognorm","type",
+ "strata_area_km2",
+ "mo","wk","season","basin","mday")]
+ nms <- c("val_use","jday","mnDate",
+ "mnVal_hind","sdVal_hind","val_raw","val_delta")
+
+ nms2 <- grp_names[!grp_names%in%nms]
+
+ futDat_avg_tmp <- futDat_weekly_strata_tmp%>%
+ select(!!!syms(grp_names))%>%group_by(!!!syms(nms2))%>%
+ summarize_at(nms,mean,na.rm=T)%>%
+ relocate(val_use, .before = mnVal_hind)%>%
+ relocate(var_salmon, .before = var)%>%
+ ungroup()%>%data.frame()
+
+ jj <- jj + 1
+ if(jj==1){
+ futDat_avg <- futDat_avg_tmp
+ futDat_monthly <- futDat_monthly_tmp
+ futDat_weekly_strata <- futDat_weekly_strata_tmp
+
+ }else{
+ futDat_avg <- rbind(futDat_avg,futDat_avg_tmp)
+ futDat_monthly <- rbind(futDat_monthly,futDat_monthly_tmp)
+ futDat_weekly_strata <- rbind(futDat_weekly_strata,futDat_weekly_strata_tmp)
+ }
+
+ rm(futDat_avg_tmp)
+ rm(futDat_monthly_tmp)
+ rm(futDat_weekly_strata_tmp)
+
+
+ }
+ rm(fut)
+ }
+}
+
+ fix_caps <- function (dataset){
+ out <- dataset%>%
+ mutate(GCM = gsub("MIROC","miroc",GCM))%>%
+ mutate(GCM = gsub("GFDL","gfdl",GCM))%>%
+ mutate(GCM = gsub("CESM","cesm",GCM))%>%
+ mutate(sim = gsub("MIROC","miroc",sim))%>%
+ mutate(sim = gsub("GFDL","gfdl",sim))%>%
+ mutate(sim = gsub("CESM","cesm",sim))%>%
+ mutate(GCM = factor(GCM, levels =c("hind","gfdl","cesm","miroc")))%>%
+ mutate(GCM_scen = paste0(GCM,"_",RCP))%>%data.frame()
+ return(out)
+ }
+
+ cat("saving outputs\n")
+
+ # fix the caps issue:
+ futDat_monthly <- fix_caps(futDat_monthly)
+ futDat_weekly_strata <- fix_caps(futDat_weekly_strata)
+ futDat_avg <- fix_caps(futDat_avg)
+
+ # output the data as Rdata and CSV
+ write_csv(futDat_monthly, file.path(outfldr,"futDat_monthly.csv"))
+ save(futDat_monthly, file=file.path(outfldr,"futDat_monthly.Rdata"))
+
+ # output the data as Rdata and CSV
+ write_csv(futDat_weekly_strata, file.path(outfldr,"futDat_weekly_strata.csv"))
+ save(futDat_weekly_strata, file=file.path(outfldr,"futDat_weekly_strata.Rdata"))
+
+ # output the data as Rdata and CSV
+ write_csv(futDat_avg, file.path(outfldr,"futDat_avg.csv"))
+ save(futDat_avg, file=file.path(outfldr,"futDat_avg.Rdata"))
+
+### ----- PREV
+
+
+if (1 ==10){
+ # loads packages, data, setup, etc.
+ suppressMessages( suppressWarnings(source("R/make.R")))
+
+ # get the variable you want:
+ df <- get_var( typeIN = "monthly",
+ monthIN = 9:10,
+ plotvar = "temp_bottom5m",
+ bcIN = c("raw","bias corrected"),
+ CMIPIN = "K20P19_CMIP6",
+ plothist = T, # ignore the hist runs
+ removeyr1 = T) # "Remove first year of projection ( burn in)")
+
+ df$plot+coord_cartesian(ylim = c(0, 7))
+ head(df$dat)
+
+ # concat the hind and fut runs by removing years from projection
+ stitchDate <- "2020-12-30"
+
+ newdat <- stitchTS(dat = df$dat,
+ maxD = stitchDate)
+
+ # newdat has the full set of data
+ # select miroc_ssp126
+ head(newdat%>%dplyr::filter(GCM_scen==paste0(GCMs[1],"_",scens[1])))
+ tail(newdat%>%dplyr::filter(GCM_scen==paste0(GCMs[1],"_",scens[1])))
+
+ pp <- plotTS(newdat )
+ pp
+
+}
\ No newline at end of file
diff --git a/R/sub_scripts/older_scripts/section5_chunk4Andy.R b/R/sub_scripts/older_scripts/section5_chunk4Andy.R
old mode 100644
new mode 100755
diff --git a/R/sub_scripts/plotACLIM2_indices.R b/R/sub_scripts/plotACLIM2_indices.R
old mode 100644
new mode 100755
diff --git a/R/sub_scripts/plot_BC_station.R b/R/sub_scripts/plot_BC_station.R
old mode 100644
new mode 100755
diff --git a/R/sub_scripts/plot_BC_stratawk.R b/R/sub_scripts/plot_BC_stratawk.R
old mode 100644
new mode 100755
diff --git a/R/sub_scripts/plot_NEBSnSEBS.R b/R/sub_scripts/plot_NEBSnSEBS.R
old mode 100644
new mode 100755
index 879b6d8..dfbacf0
--- a/R/sub_scripts/plot_NEBSnSEBS.R
+++ b/R/sub_scripts/plot_NEBSnSEBS.R
@@ -2,6 +2,8 @@
#'
#'plot_NEBSnSEBS.R
#'
+#'
+
fldrout <- "Figs/prod_plots"
if(dir.exists(fldrout))
dir.remove(fldrout)
@@ -18,6 +20,7 @@
vardef <- vardef%>%filter(name%in%unique(ACLIM_weekly_hind$var))
+
i <-grep("temp_bottom5m",vardef$name)
w<-9; h<-6; dpi <-350
w2<-6; h2<-4
@@ -50,9 +53,10 @@
rm(tmp2)
rm(tmp)
- tmp3 <- suppressMessages(plotNEBS_productivity_futBC(datIN=ACLIM_weekly_hind,alphaIN = .4,
- datIN_fut=ACLIM_weekly_fut,
- varlistIN=vardef[i,]))
+ tmp3 <- suppressMessages(plotNEBS_productivity_futBC(datIN = ACLIM_weekly_hind,
+ alphaIN = .4,
+ datIN_fut = ACLIM_weekly_fut%>%mutate(scen=paste0(RCP,"_",GCM)),
+ varlistIN = vardef[i,]))
jpeg(file.path(fldrout,paste0(vv,"_futBC.jpg")),width = w*sclr, height=h*sclr, res=dpi, units="in")
print(tmp3$p2)
diff --git a/R/sub_scripts/plot_RKC_Indices_Andre.R b/R/sub_scripts/plot_RKC_Indices_Andre.R
new file mode 100644
index 0000000..1ab66dd
--- /dev/null
+++ b/R/sub_scripts/plot_RKC_Indices_Andre.R
@@ -0,0 +1,269 @@
+#'
+#'
+#'
+#'plot_RKC_Indices_Andre
+#'ALCIM3 K. Holsman
+#'2024
+#'
+
+# rmlist=ls()
+message("Set outfldr to your local folder with the RKC indiceies \n e.g., /Users/KKH/Documents/GitHub_mac/ACLIM2/Data/out/RKC_indices")
+outfldr <- "Data/out/RKC_indices"
+
+# ------------------------------------
+# Load packages
+# ------------------------------------
+load(file.path(outfldr,"RKC_indices_setup.Rdata"))
+library(ggplot2)
+library(dplyr)
+library(RColorBrewer)
+#
+# # load ACLIM packages and functions
+# suppressMessages(source("R/make.R"))
+
+# ------------------------------------
+# Load Data
+# ------------------------------------
+
+
+
+load(file.path(outfldr,"AVG_hind_TC_small_males.Rdata"))
+load(file.path(outfldr,"AVG_hind_RKC_immature_males.Rdata"))
+load(file.path(outfldr,"AVG_fut_TC_small_males.Rdata"))
+load(file.path(outfldr,"AVG_fut_RKC_immature_males.Rdata"))
+
+load(file.path(outfldr,"fut_TC_small_males_monthly.Rdata"))
+load(file.path(outfldr,"fut_TC_small_males_monthly_avg.Rdata"))
+load(file.path(outfldr,"fut_RKC_immature_males_monthly.Rdata"))
+load(file.path(outfldr,"fut_RKC_immature_males_monthly_avg.Rdata"))
+
+load(file.path(outfldr,"hind_TC_small_males_monthly.Rdata"))
+load(file.path(outfldr,"hind_TC_small_males_monthly_avg.Rdata"))
+load(file.path(outfldr,"hind_RKC_immature_males_monthly.Rdata"))
+load(file.path(outfldr,"hind_RKC_immature_males_monthly_avg.Rdata"))
+
+load(file.path(outfldr,"hind_TC_small_males.Rdata"))
+load(file.path(outfldr,"hind_RKC_immature_males.Rdata"))
+load(file.path(outfldr,"fut_TC_small_males.Rdata"))
+load(file.path(outfldr,"fut_RKC_immature_males.Rdata"))
+
+# ------------------------------------
+# Plot results
+# ------------------------------------
+
+ # Set color palette
+ # ------------------------------------
+ display.brewer.all()
+
+ GCM_scens <- c("hind",unique(AVG_fut_RKC_immature_males$GCM_scen))
+ GCM_scens <- GCM_scens[
+ c(1,
+ grep("ssp126",GCM_scens),
+ grep("rcp45",GCM_scens),
+ grep("rcp85",GCM_scens),
+ grep("ssp585",GCM_scens))]
+
+ cc <- brewer.pal(n = 11, name ="Spectral") [c(10,9,4,3)]
+
+ nn <- length(grep("ssp126",GCM_scens))
+ c1 <- rep(cc[1],nn)
+ nn <- length(grep("rcp45",GCM_scens))
+ c2 <- rep(cc[2],nn)
+ nn <- length(grep("rcp85",GCM_scens))
+ c3 <- rep(cc[3],nn)
+ nn <- length(grep("ssp585",GCM_scens))
+ c4 <- rep(cc[4],nn)
+
+ cols <- c("gray",c1,c2,c3,c4)
+ valuesIN <- cols
+ names(valuesIN) <- GCM_scens
+
+ # Plot AVG across subset stations:
+ # --------------------------------------
+ AVGfut_RKC <- AVG_fut_RKC_immature_males
+ plot_stationAVG <- plotTS(AVGfut_RKC%>%
+ mutate(mn_val = mn_val_use,
+ scen = GCM,
+ GCM_scen = factor(GCM_scen,levels = GCM_scens)))+
+ labs(title="average variables (RKC_immature_males) ")+
+ scale_color_manual(name = "Scenario",
+ values = valuesIN)+
+ scale_fill_manual(name = "Scenario",
+ values = valuesIN)
+ plot_stationAVG
+
+ # Plot Indiv station:
+ # --------------------------------------
+ stN <- 3
+ fut_RKC <- fut_RKC_immature_males
+ station_set <- unique(fut_RKC$station_id)
+
+ plot_station <- plotTS(fut_RKC%>%mutate(mn_val=val_use,
+ scen=GCM,
+ GCM_scen = factor(GCM_scen,levels = GCM_scens))%>%
+ filter(station_id==station_set[stN]) )+
+ labs(title=paste("Station: ", station_set[stN]))+
+ scale_color_manual(name = "Scenario",
+ values = valuesIN)+
+ scale_fill_manual(name = "Scenario",
+ values = valuesIN)
+ plot_station
+
+
+# ------------------------------------
+# Plot restults with hindcast
+# ------------------------------------
+
+ # Plot AVG across subset stations:
+ # --------------------------------------
+ stN <- 3 # for further down
+ vlN <- 1:length(varlist)
+
+
+ AVGfut_RKC <- AVG_fut_RKC_immature_males%>%mutate(mn_val = mn_val_use,
+ scen = GCM,
+ GCM_scen = factor(GCM_scen,levels = GCM_scens))
+ AVGhind_RKC <- AVG_hind_RKC_immature_males%>%mutate(mn_val = mn_val_use,
+ scen = "hind",
+ GCM_scen = factor("hind",levels = GCM_scens))
+
+ hindfut_avg <- ggplot( )+
+ geom_line(data=AVGfut_RKC%>%mutate(scen=GCM)%>%filter(var%in%varlist[vlN]),
+ aes(x=mnDate,y=mn_val_use,color= GCM_scen,linetype = basin),
+ alpha = 0.2,show.legend = FALSE)+
+ geom_smooth(data=AVGfut_RKC%>%mutate(scen=GCM)%>%filter(var%in%varlist[vlN]),
+ aes(x=mnDate,y=mn_val_use,color= GCM_scen,
+ fill=GCM_scen,linetype = basin),
+ alpha=0.1,method="loess",formula='y ~ x',span = .5)+facet_grid(var~GCM,scales="free_y")+
+
+ geom_line(data=AVGhind_RKC%>%filter(var%in%varlist[vlN]),
+ aes(x=mnDate,y=mn_val_use,color="hind",linetype = basin),
+ alpha = 0.2,show.legend = FALSE)+
+ geom_smooth(data=AVGhind_RKC%>%filter(var%in%varlist[vlN]),
+ aes(x=mnDate,y=mn_val_use,color= "hind",
+ fill="hind",linetype = basin),
+ alpha=0.1,method="loess",formula='y ~ x',span = .5)+
+ theme_minimal() +
+ labs(x="Year",
+ subtitle = "",
+ title = "cross-station average (RKC_immature_males)",
+ legend = "")+
+ scale_color_manual(name = "Scenario",
+ values = valuesIN)+
+ scale_fill_manual(name = "Scenario",
+ values = valuesIN)
+
+ hindfut_avg
+
+
+
+ # Plot Indiv station:
+ # --------------------------------------
+ hind_RKC <- hind_RKC_immature_males%>%mutate(mn_val = val_use,
+ scen = "hind",
+ GCM_scen = factor("hind",levels = GCM_scens))
+ fut_RKC <- fut_RKC_immature_males%>%mutate(mn_val = val_use,
+ scen = GCM,
+ GCM_scen = factor(GCM_scen,levels = GCM_scens))
+
+
+ hindfut_station <- ggplot( )+
+ geom_line(data=fut_RKC%>%mutate(mn_val=val_use,scen=GCM)%>%filter(station_id==station_set[stN],var%in%varlist[vlN]),
+ aes(x=mnDate,y=val_use,color= GCM_scen,linetype = basin),
+ alpha = 0.2,show.legend = FALSE)+
+ geom_smooth(data=fut_RKC%>%mutate(mn_val=val_use,scen=GCM)%>%filter(station_id==station_set[stN],var%in%varlist[vlN]),
+ aes(x=mnDate,y=val_use,color= GCM_scen,
+ fill=GCM_scen,linetype = basin),
+ alpha=0.1,method="loess",formula='y ~ x',span = .5)+facet_grid(var~GCM,scales="free_y")+
+
+ geom_line(data=hind_RKC%>%filter(station_id%in%station_set[stN],var%in%varlist[vlN]),
+ aes(x=mnDate,y=val_use,color="hind",linetype = basin),
+ alpha = 0.2,show.legend = FALSE)+
+ geom_smooth(data=hind_RKC%>%filter(station_id==station_set[stN],var%in%varlist[vlN]),
+ aes(x=mnDate,y=val_use,color= "hind",
+ fill="hind",linetype = basin),
+ alpha=0.1,method="loess",formula='y ~ x',span = .5)+
+ theme_minimal() +
+ labs(x="Year",
+ subtitle = "",
+ title = paste("Station:", station_set[stN]),
+ legend = "")+
+ scale_color_manual(name = "Scenario",
+ values = valuesIN)+
+ scale_fill_manual(name = "Scenario",
+ values = valuesIN)
+
+ hindfut_station
+
+
+ # plot monthly indices
+ # ------------------------------------
+
+ # ggplot(hind_TC_small_males_monthly)+
+ # geom_line(aes(x= mnDate,y = val_use,color=var))+facet_grid(var~strata,scales="free_y")
+ #
+ stratalist<-unique( fut_TC_small_males_monthly$strata)
+ sin <- 1
+ hindfut_station_monthly <- ggplot()+
+ geom_line(data=fut_TC_small_males_monthly%>%mutate(mn_val=val_use,scen=GCM)%>%
+ filter(strata==stratalist[sin],var%in%varlist[vlN]),
+ aes(x=mnDate,y=val_use,color= GCM_scen,linetype = basin),
+ alpha = 0.2,show.legend = FALSE)+
+ geom_smooth(data=fut_TC_small_males_monthly%>%mutate(mn_val=val_use,scen=GCM)%>%
+ filter(strata==stratalist[sin],var%in%varlist[vlN]),
+ aes(x=mnDate,y=val_use,color= GCM_scen,
+ fill=GCM_scen,linetype = basin),
+ alpha=0.1,method="loess",formula='y ~ x',span = .5)+
+ facet_grid(var~GCM,scales="free_y")+
+ geom_line(data=hind_TC_small_males_monthly%>%filter(strata==stratalist[sin],var%in%varlist[vlN]),
+ aes(x=mnDate,y=val_use,color="hind",linetype = basin),
+ alpha = 0.2,show.legend = FALSE)+
+ geom_smooth(data=hind_TC_small_males_monthly%>%filter(strata==stratalist[sin],var%in%varlist[vlN]),
+ aes(x=mnDate,y=val_use,color= "hind",
+ fill="hind",linetype = basin),
+ alpha=0.1,method="loess",formula='y ~ x',span = .5)+
+ theme_minimal() +
+ labs(x="Year",
+ subtitle = "",
+ title = paste("Monthly by strata:",stratalist[sin] ),
+ legend = "")+
+ scale_color_manual(name = "Scenario",
+ values = valuesIN)+
+ scale_fill_manual(name = "Scenario",
+ values = valuesIN)
+
+ hindfut_station_monthly
+
+
+ # Save plots
+ # ------------------------------------
+
+ sclr <- 1.2
+ jpeg(filename = file.path(outfldr,"Station_Averaged_RKC_vars.jpg"),
+ width=8*sclr, height=7*sclr, units="in", res = 350)
+ print(plot_stationAVG)
+ dev.off()
+
+ jpeg(filename = file.path(outfldr,paste0("RKC_vars_",station_set[stN],".jpg")),
+ width=8*sclr, height=7*sclr, units="in", res = 350)
+ print(plot_station)
+ dev.off()
+
+
+ sclr <- 1.2
+ jpeg(filename = file.path(outfldr,"RKC_indices_station_W_hind.jpg"),
+ width=8*sclr,height=7*sclr,units="in",res=350)
+ print(hindfut_station)
+ dev.off()
+
+ jpeg(filename = file.path(outfldr,"RKC_indices_avg_W_hind.jpg"),
+ width=8*sclr,height=7*sclr,units="in",res=350)
+ print(hindfut_avg)
+ dev.off()
+
+ jpeg(filename = file.path(outfldr,"TC_small_males_monthly.jpg"),
+ width=8*sclr,height=7*sclr,units="in",res=350)
+ print(hindfut_station_monthly)
+ dev.off()
+
+
\ No newline at end of file
diff --git a/R/sub_scripts/plot_salmon_indices.R b/R/sub_scripts/plot_salmon_indices.R
new file mode 100644
index 0000000..00d6821
--- /dev/null
+++ b/R/sub_scripts/plot_salmon_indices.R
@@ -0,0 +1,213 @@
+#'
+#'
+#'
+#'plot_salmon_Indices
+#'ALCIM3 K. Holsman
+#'2024
+#'
+
+# rmlist=ls()
+message("Set outfldr to your local folder with the salmon indiceies \n e.g., /Users/KKH/Documents/GitHub_mac/ACLIM2/Data/out/salmon_indices")
+outfldr <- "Data/out/salmon_indices"
+
+# ------------------------------------
+# Load packages
+# ------------------------------------
+load(file.path(outfldr,"salmon_indices_setup.Rdata"))
+library(ggplot2)
+library(dplyr)
+library(RColorBrewer)
+#
+# # load ACLIM packages and functions
+# suppressMessages(source("R/make.R"))
+
+# ------------------------------------
+# Load Data
+# ------------------------------------
+
+dirlist <- dir(outfldr)
+dirlist <- dirlist[grep("Dat",dirlist)]
+dirlist <- dirlist[grep("Rdata",dirlist)]
+
+for(d in dirlist)
+ load(file.path(outfldr,d))
+
+# ------------------------------------
+# Plot results
+# ------------------------------------
+
+# Set color palette
+# ------------------------------------
+display.brewer.all()
+
+GCM_scens <- c("hind",unique(futDat_weekly_strata$GCM_scen))
+GCM_scens <- GCM_scens[
+ c(1,
+ grep("ssp126",GCM_scens),
+ grep("rcp45",GCM_scens),
+ grep("rcp85",GCM_scens),
+ grep("ssp585",GCM_scens))]
+
+cc <- brewer.pal(n = 11, name ="Spectral") [c(10,9,4,3)]
+
+nn <- length(grep("ssp126",GCM_scens))
+c1 <- rep(cc[1],nn)
+nn <- length(grep("rcp45",GCM_scens))
+c2 <- rep(cc[2],nn)
+nn <- length(grep("rcp85",GCM_scens))
+c3 <- rep(cc[3],nn)
+nn <- length(grep("ssp585",GCM_scens))
+c4 <- rep(cc[4],nn)
+
+cols <- c("gray",c1,c2,c3,c4)
+valuesIN <- cols
+names(valuesIN) <- GCM_scens
+
+# Plot AVG across subset stations:
+# --------------------------------------
+var_salmons <- unique(futDat_avg$var_salmon)
+
+AVGfut <- futDat_avg
+plot_stationAVG <- plotTS(AVGfut%>%
+ mutate(mn_val = val_use,
+ scen = RCP,
+ basin = var_salmon,
+ GCM_scen = factor(GCM_scen,levels = GCM_scens)))+
+ labs(title="average variables (RKC_immature_males) ")+
+ scale_color_manual(name = "Scenario",
+ values = valuesIN)+
+ scale_fill_manual(name = "Scenario",
+ values = valuesIN)+
+ facet_grid(var~scen,scales='free_y')
+plot_stationAVG
+
+
+plot_stationAVG2 <- plotTS(AVGfut%>%
+ mutate(mn_val = val_use,
+ scen = GCM,
+ basin = var_salmon,
+ GCM_scen = factor(GCM_scen,levels = GCM_scens)))+
+ labs(title="average variables (RKC_immature_males) ")+
+ scale_color_manual(name = "Scenario",
+ values = valuesIN)+
+ scale_fill_manual(name = "Scenario",
+ values = valuesIN)+
+ facet_grid(var~scen,scales='free_y')
+plot_stationAVG2
+
+
+# ------------------------------------
+# Plot restults with hindcast
+# ------------------------------------
+
+# Plot AVG across subset stations:
+# --------------------------------------
+varlist <- unique(futDat_avg$var)
+var_sal_set <- var_salmons
+
+
+vlN <- 1:length(varlist)
+
+
+AVGfut <- futDat_avg%>%mutate(mn_val = val_use,
+ scen = GCM,
+ GCM_scen = factor(GCM_scen,levels = GCM_scens))
+AVGhind_salmon <- hindDat_avg%>%mutate(mn_val = val_use,
+ scen = "hind",
+ GCM_scen = factor("hind",levels = GCM_scens))
+
+hindfut_avg <- ggplot( )+
+ geom_line(data=AVGfut%>%mutate(scen=GCM)%>%filter(var%in%varlist[vlN]),
+ aes(x=mnDate,y=val_use,color= GCM_scen,linetype = var_salmon),
+ alpha = 0.2,show.legend = FALSE)+
+ geom_smooth(data=AVGfut%>%mutate(scen=GCM)%>%filter(var%in%varlist[vlN]),
+ aes(x=mnDate,y=val_use,color= GCM_scen,
+ fill=GCM_scen,linetype = var_salmon),
+ alpha=0.1,method="loess",formula='y ~ x',span = .5)+facet_grid(var~GCM,scales="free_y")+
+
+ geom_line(data=AVGhind_salmon%>%filter(var%in%varlist[vlN]),
+ aes(x=mnDate,y=val_use,color="hind",linetype = var_salmon),
+ alpha = 0.2,show.legend = FALSE)+
+ geom_smooth(data=AVGhind_salmon%>%filter(var%in%varlist[vlN]),
+ aes(x=mnDate,y=val_use,color= "hind",
+ fill="hind",linetype = var_salmon),
+ alpha=0.1,method="loess",formula='y ~ x',span = .5)+
+ theme_minimal() +
+ labs(x="Year",
+ subtitle = "",
+ title = "cross-station average (RKC_immature_males)",
+ legend = "")+
+ scale_color_manual(name = "Scenario",
+ values = valuesIN)+
+ scale_fill_manual(name = "Scenario",
+ values = valuesIN)
+
+hindfut_avg
+
+
+
+# Plot Indiv station:
+# --------------------------------------
+hind_salmon <- hindDat_avg%>%mutate(mn_val = val_use,
+ scen = "hind",
+ GCM_scen = factor("hind",levels = GCM_scens))
+fut_salmon <- futDat_avg%>%mutate(mn_val = val_use,
+ scen = GCM,
+ GCM_scen = factor(GCM_scen,levels = GCM_scens))
+
+
+stN <- 3
+for(stN in 1:length(var_sal_set)){
+ var_sal_set[stN]
+
+ hindfut_station <- ggplot( )+
+ geom_line(data=fut_salmon%>%mutate(mn_val=val_use,scen=GCM)%>%filter(var_salmon==var_sal_set[stN],var%in%varlist[vlN]),
+ aes(x=mnDate,y=val_use,color= GCM_scen,linetype = var_salmon),
+ alpha = 0.2,show.legend = FALSE)+
+ geom_smooth(data=fut_salmon%>%mutate(mn_val=val_use,scen=GCM)%>%filter(var_salmon==var_sal_set[stN],var%in%varlist[vlN]),
+ aes(x=mnDate,y=val_use,color= GCM_scen,
+ fill=GCM_scen,linetype = var_salmon),
+ alpha=0.1,method="loess",formula='y ~ x',span = .5)+facet_grid(var~GCM,scales="free_y")+
+
+ geom_line(data=hind_salmon%>%filter(var_salmon%in%var_sal_set[stN],var%in%varlist[vlN]),
+ aes(x=mnDate,y=val_use,color="hind",linetype = var_salmon),
+ alpha = 0.2,show.legend = FALSE)+
+ geom_smooth(data=hind_salmon%>%filter(var_salmon==var_sal_set[stN],var%in%varlist[vlN]),
+ aes(x=mnDate,y=val_use,color= "hind",
+ fill="hind",linetype = var_salmon),
+ alpha=0.1,method="loess",formula='y ~ x',span = .5)+
+ theme_minimal() +
+ labs(x="Year",
+ subtitle = "",
+ title = paste("Salmon strata/mo:", var_sal_set[stN]),
+ legend = "")+
+ scale_color_manual(name = "Scenario",
+ values = valuesIN)+
+ scale_fill_manual(name = "Scenario",
+ values = valuesIN)
+
+ hindfut_station
+
+ jpeg(filename = file.path(outfldr,paste0("salmon_indices_indiv_hind",stN,".jpg")),
+ width=8*sclr,height=7*sclr,units="in",res=350)
+ print(hindfut_avg)
+ dev.off()
+}
+
+
+sclr <- 1.2
+jpeg(filename = file.path(outfldr,"salmon_vars.jpg"),
+ width=8*sclr, height=7*sclr, units="in", res = 350)
+print(plot_stationAVG)
+dev.off()
+
+
+sclr <- 1.2
+jpeg(filename = file.path(outfldr,"salmon_indices_fut_hind.jpg"),
+ width=8*sclr,height=7*sclr,units="in",res=350)
+print(hindfut_station)
+dev.off()
+
+
+
+
diff --git a/R/sub_scripts/profile_attach.R b/R/sub_scripts/profile_attach.R
old mode 100644
new mode 100755
diff --git a/R/sub_scripts/updateREADME.R b/R/sub_scripts/updateREADME.R
old mode 100644
new mode 100755
diff --git a/R/sub_scripts/update_base_data.R b/R/sub_scripts/update_base_data.R
old mode 100644
new mode 100755
diff --git a/R/sub_scripts/xplot_weekly.R b/R/sub_scripts/xplot_weekly.R
old mode 100644
new mode 100755
diff --git a/README.md b/README.md
old mode 100644
new mode 100755
diff --git a/Rplot.png b/Rplot.png
old mode 100644
new mode 100755
diff --git a/Vignettes/Ellen_regional_temperture_index.Rmd b/Vignettes/Ellen_regional_temperture_index.Rmd
old mode 100644
new mode 100755
diff --git a/Vignettes/Ellen_regional_temperture_index.pdf b/Vignettes/Ellen_regional_temperture_index.pdf
old mode 100644
new mode 100755
diff --git a/Vignettes/Ellen_regional_temperture_index.tex b/Vignettes/Ellen_regional_temperture_index.tex
old mode 100644
new mode 100755
diff --git a/_config.yml b/_config.yml
old mode 100644
new mode 100755
diff --git a/aclim.png b/aclim.png
old mode 100644
new mode 100755
diff --git a/index.html b/index.html
old mode 100644
new mode 100755
diff --git a/nul b/nul
new file mode 100755
index 0000000..ac1d590
--- /dev/null
+++ b/nul
@@ -0,0 +1,325 @@
+
+
+
+
+
+ The CTAN archive
+
+
+
+ Comprehensive TeX Archive Network
+ The CTAN root directory
+
+
+ This is the root of the Comprehensive
+ TeX Archive Network directory tree.
+ Most information in CTAN is contained in subdirectories:
+
+
+
+ - biblio
+ -
+ Systems for maintaining and presenting bibliographies within documents
+ typeset using TeX
+
+
+ - digests
+ - Collections of
+ TeX mailing list digests,
+ TeX-related electronic magazines, and indexes,
+ etc., of printed publications
+
+
+ - dviware
+ -
+ Printer drivers and previewers, etc., for .dvi files
+
+
+ - fonts
+ -
+ Fonts written in Metafont, and support for using fonts from other
+ sources (e.g., those in Adobe Type 1 format)
+
+
+ - graphics
+ -
+ Systems and TeX macros for producing graphics
+
+
+ - help
+ -
+ FAQs
+ (de
+ fr
+ uk)
+ and similar direct assistance
+
+
+ - indexing
+ -
+ Systems for maintaining and presenting indexes of documents typeset
+ using TeX
+
+
+ - info
+ -
+ Manuals and extended how-to information; errata for
+ TeX-related publications, collections of project
+ (e.g., LaTeX and NTS)
+ documents, etc.
+
+
+ - install
+ -
+ ZIP files of ready-to-install CTAN packages. The package contained in
+ CTAN directory a/b/c/d/ appears as
+ install/a/b/c/d.tds.zip – each such file may be unzipped
+ directly into
+ a TDS
+ tree, and should contain no “tree prefix”
+
+
+ - language
+ -
+ Support for various languages; hyphenation patterns
+
+
+ - macros
+ - TeX
+ macros; several directories have significant sub-trees:
+
+ - context
+ -
+ The ConTeXt distribution
+ and contributed material
+
+
+ - generic
+ - Macros that work in several environments
+
+ - latex
+ -
+ The LaTeX
+ distribution and contributed material
+
+
+ - luatex
+ - Macros that require the LuaTeX engine
+
+ - plain
+ - Donald Knuth's example macro set and contributed material
+
+ - xetex
+ - Macros that require the XƎTeX engine
+
+
+
+ - obsolete
+ -
+ Material which is now obsolete, but which needs to remain
+ on the archive for some reason (such as the 1992
+ LaTeX 2.09
+ distribution). The tree maintains a “shadow hierarchy”
+ – for example obsolete/macros
+ has the same basic structures macros.
+
+
+ - support
+ -
+ TeX support environments and the like
+
+
+ - systems
+ -
+ TeX systems;
+ organised by operating environment, but also including:
+
+
+ - knuth
+ - Donald Knuth's current distribution.
+
+ - generic
+ -
+ Complete systems that can potentially operate in more than one
+ operating environment.
+
+
+ - MacTeX
+ - The distribution of MacTeX for Mac OS X
+
+ -
+ TeX Live
+
+ -
+ The Live collection of freely distributable TeX
+ material, including the TeX
+ Live on-line installer.
+
+
+ - MiKTeX
+ -
+ The distribution of MiKTeX; CTAN provides a
+ home for MiKTeX's package repository, part
+ of its online installer.
+
+
+
+
+ - tds
+ -
+ The TeX Directory Structure standard
+ (the output of the TUG TDS working group)
+
+
+ - usergrps
+ - Information supplied by TeX User Groups
+
+ - web
+ -
+ Literate Programming tools and systems
+
+
+
+ CTAN on the Web
+
+ The home page for CTAN is at
+ https://ctan.org,
+ which provides several features:
+
+
+
+ Searching CTAN
+
+ You can search the content of the Catalogue with the
+ following form on the CTAN site:
+
+
+
+
+
+ In addition, CTAN archives maintain lists of their contents
+ as plain text files, in a number of ways:
+
+
+ -
+ Sorted by name
+ (This is a big file…around 20M bytes).
+
+ -
+ New or changed in the last 7 days
+ (This is usually a more manageable size…)
+
+
+
+
+
+ This is an experimental page. Comments, to
+ the CTAN team,
+ will be most welcome.
+
+ Last updated 2020-03-31
+
+
+
diff --git a/plot_ts.R b/plot_ts.R
old mode 100644
new mode 100755
diff --git a/shiny/._NPSSP b/shiny/._NPSSP
deleted file mode 100644
index 169e65d..0000000
Binary files a/shiny/._NPSSP and /dev/null differ
diff --git a/shiny/NPSSP/NPSSP.R b/shiny/NPSSP/NPSSP.R
old mode 100644
new mode 100755
diff --git a/~$LIM2_quickStart.docx b/~$LIM2_quickStart.docx
deleted file mode 100644
index a837b48..0000000
Binary files a/~$LIM2_quickStart.docx and /dev/null differ