## Bigheaded Carp Growth Rate Potential (GRP) model code: Foraging and Bioenergetic functions ## Function of this script: Calculate suitable habitat for each month and mean GRP in suitable habitat for the model year 2010 ## Creator: Peter Alsip, UM School for Environment and Sustainability ## Submission date: 9/7/2018 ## Set species # species<-"bighead" #species<-"silver" ###################### Read in Biophysical Model outputs ###################### ## note: This requires access to FVCOM-GEM outputs for 2010 in order for this script to be functional. ## As it is currently written, this script is not functional without FVCOM-GEM outputs and therefore ## only contains model equations that were developed independently from the Biophysical model. ## Define path to read in input data from biophysical model #input_path<- ## Open all files ncfiles <- list.files(input_path,pattern='*.nc',full.names=TRUE) ncfiles <- ncfiles[which(nchar(basename(ncfiles)) == 10)] ## Load grid data load("E:/Polygons_Grid/shoreline_latlon.Rdata") load("E:/Polygons_Grid/Mich03_grid_polys.Rdata") load("E:/Polygons_Grid/Mich03_grid_coords.Rdata") load("E:/Polygons_Grid/Mich03_grid_coords_gb.Rdata")# green bay nodes (1 = in GB, 0= not in GB) ######## Simulate Growth on monthly, non-dynamic time step (growth.fun) and on a daily, dynamic time step (yearGrowth.fun) ## Calculate growth rate potential (species == "bighead" or "silver", month_num [1...12], file_list == FVCOM-GEM output netcdf files) ## note: growth.fun is not dynamic (i.e. if you run it for months 1-12, the start weight for each time step does not change) growth.fun = function(species, month_num, file_list){ print(paste("Preparing dataframe w/ GRP values for", species, "in", month.name[month_num])) require(ncdf4) load("E:/Polygons_Grid/shoreline_latlon.Rdata") load("E:/Polygons_Grid/Mich03_grid_polys.Rdata") load("E:/Polygons_Grid/Mich03_grid_coords.Rdata") load("E:/Polygons_Grid/Mich03_grid_coords_gb.Rdata") # Assign ncfile from month of interest to nc nc <- file_list[month_num] # Open ncfile nc <- nc_open(nc, readunlim = FALSE) PP_carbon <- ncvar_get(nc, varid = "Phytoplankton", start = c(1,1,371), count = c(-1,-1,1)) ZP_carbon <- ncvar_get(nc, varid = "Zooplankton", start = c(1,1,371), count = c(-1,-1,1)) Detritus <- ncvar_get(nc, varid = "Detritus", start = c(1,1,371), count = c(-1,-1,1)) temp <- ncvar_get(nc, varid = "temp", start = c(1,1,371), count = c(-1, -1, 1)) nc_close(nc) temp[which(temp == 0)]<-NA # Areas with Ice are not possible habitat thus the final Growth output should be NA ED_PP <- 2600 #From Anderson et al. 2015; ED of Green algae (J/g) ED_ZP <- 2512 #Cooke and Hill 2010 use 2512 J/g of WET biomass for ED_ZP (from Cummins and Wuycheck 1971) ED_Det<- 979*.13 #Value and Adjustment factor derived from Anderson et al. 2016 #Convert PP_carbon (mgC/L) to wet phytoplankton biomass (g/L) PP <- (PP_carbon*36)/(1000) #1/.0278 = 36 Chla<-PP_carbon*.036*1000 #Convert ZP_carbon (mgC/L) to Wet zooplankton biomass (g/L) ZP <- (ZP_carbon * 2.5/.1)/1000 # 0.40 is C mass:dry mass (1/0.40=2.5), Det<- (Detritus/1000)/.044/.16 #(g/L) #.044 from Ozersky et al.; .16 dry:wet from Anderson et al. 2016 # Shared parameters # Cons eqn 2 CA <- .369 # Anderson et al. 2017 CB <- -.225 #Anderson et al. 2015, 2017 (taken from Wang 1989) CQ <- 2.5 #Value for bighead and silver from Kolar et al. 2007 # Resp eqn 1 SDA <- 0.1 #Tilapia value from Cooke and Hill 2010 RTO <- 0 ACT <- 1 # Waste eqn 2 FA <- .212 FB <- -.222 FG <- .631 UA <- .031 UB <- .58 UG <- -.299 if (species == "bighead"){ CTO <- 26 #Bighead, from Kolar et al. 2007 CTM <- 38 #Bighead, Kolar et al. 2007 W <- 5480 #Bighead wt. from Anderson et al. 2015 ED_carp <- 3500 #J/g Anderson et al. 2015 RA <- .00528 #Bighead, from Cooke and Hill 2010 RB <- -.299 #Bighead, from Cooke and Hill 2010 RQ <- .048 } if (species == "silver"){ CTO <- 29 #Silver, from Kolar et al. 2007 CTM <- 43 #Silver, Kolar et al. 2007 W <- 4350 #Silver wt. from Anderson et al. 2015 ED_carp <- 5200 #J/g Silver Anderson et al. 2015 RA <- .00279 #Silver, from Cooke and Hill 2010 RB <- -.239 #Silver, from Cooke and Hill 2010 RQ <- .076 #Silver, from Cooke and Hill 2010 } #Temperature function (f_t) V <- (CTM - temp)/(CTM - CTO) Z <- log(CQ) * (CTM-CTO) Y <- log(CQ) * (CTM-CTO+2) X <- ((Z^2) * (1+(1+40/Y)^0.5)^2)/400 # Calculate Cmax and f_t Cmax <- CA*W^CB f_t <- (V^X)*exp(X*(1-V)) # Filtration rate equation from Smith 1989 FR <- 1.54 * (W^0.713) * 24 # Calculate energy density ratio between carp and prey for ZP+PP diet and ZP+PP+Det diet ED_prey <- (ED_PP * PP) + (ED_ZP * ZP) + (ED_Det*Det) # J/L ED_prey2<- (ED_PP * PP) + (ED_ZP *ZP) #J/L Prey_ed <-ED_prey/(PP+ZP+Det) #Gives wtd. average j/g of prey Prey_ed2<-ED_prey2/(PP+ZP) Prey_ed[is.nan(Prey_ed)]<-0 #Set all NaNs to 0. Prey_ed2[is.nan(Prey_ed2)]<-0 ED_ratio<- Prey_ed/ED_carp #Converts gPrey to gCarp ED_ratio2<- Prey_ed2/ED_carp # Calculate Maximum Consumption rate (gPrey/gCarp/day) Cons <- Cmax * f_t #Calculate foraging model consumption rate (gPrey/gCarp/day) C <- (FR *(PP+ZP+Det)/W)*f_t #w/Det C2<- (FR *(PP+ZP)/W)*f_t #w/o Det # Pairwise min of each consumption output, assign minimum to Cnew and multiply by ED_ratio to convert gPrey/gCarp/day to gCarp/gCarp/day Cnew1 <- pmin(C, Cons)*ED_ratio #w/ det Cnew2<- pmin(C2, Cons)*ED_ratio2 #w/o det # Calculate proportion of maximum consumption (P) P1 <- pmin(C, Cons)/Cons #w/ det P1[is.nan(P1)]<-0 P2<- pmin(C2, Cons)/Cons #w/o det P2[is.nan(P2)]<-0 #Take pairwise MAX of Cnew1, Cnew2. Assign higher value to Cnew Cnew<- pmax(Cnew1, Cnew2) #Define function that allows fish to supplement their ZP+PP diet with detritus det.ration<-function(Cnew){ #Subtract foraging-based consumption rate from maximum cons to determine how much det the carp needs to eat to maximize consumption Det_ration <- Cons-C2 #gPrey/gcarp/d, amount of detritus required to reach maximum Cons #Calculate consumption based on total available detritus in the habitat (gDet/gCarp/d) C_det<- ((FR*Det)/W)*f_t #Take pmin; the amount of det eaten can never be greater than Det_ration CnewDet<-pmin(C_det, Det_ration) *(ED_Det/ED_carp) #puts CnewDet in gCarp/gCarp/d #Add CnewDet to Cnew Cnew<- Cnew+CnewDet return(Cnew) } #If Cnew2 (w/o det) defines Cnew, and Cnew is less than Cons. Then add det mass of difference b/w Cnew and Cons Cnew<-ifelse(Cnew == Cnew2 & Cnew < Cons, det.ration(Cnew), Cnew) #determine P_det, then add to P2 to define P when Cnew == Cnew2 Det_ration <- Cons-C2 C_det<- ((FR*Det)/W)*f_t P_det<-pmin(C_det, Det_ration)/Cons P2<- P2 + P_det # Define P P<- ifelse(Cnew == Cnew1, P1, P2) # Calculate temperature function for respiration Resp_f_T<- exp(RQ*temp) # Calculate respiration (g/g/day) R <- RA * (W^RB) * Resp_f_T * ACT # Calculate Egestion F1 <- FA * (temp^FB)*(exp(FG*P) * Cnew) # Calculate Excretion U <- UA * (temp^UB)*(exp(UG*P)*(Cnew-F1)) # Calculate SDA S <- SDA * (Cnew-F1) # Calculate growth (g/g/day) G <- (Cnew-(R+S+F1+U)) # Convert to dataframe Gdf<- as.data.frame(G) return(Gdf) } ## Embedded function that calculates the daily growth and outputs the final weight at the end of a given day. ## This function is embedded in the yearGrowth.fun, which allows the model to be run continuously Annualgrowth.fun = function(species,startWeight, month_num, day_num, file_list){ print(paste("Preparing matrix w/ end weight values for", species, "in", month.name[month_num],day_num)) require(ncdf4) load("E:/Polygons_Grid/shoreline_latlon.Rdata") load("E:/Polygons_Grid/Mich03_grid_polys.Rdata") load("E:/Polygons_Grid/Mich03_grid_coords.Rdata") load("E:/Polygons_Grid/Mich03_grid_coords_gb.Rdata") # Assign ncfile from month of interest to nc nc <- file_list[month_num] # Create vectors that grab daily data from 10 AM # Open ncfile nc <- nc_open(nc, readunlim = FALSE) #Time step indices for 10 AM on each day of the month dayTimeIndex<-c(11,35,59,83,107,131,155,179,203,227,251,275,299,323,347,371,395,419,443,467, 491,515,539,563,587,611,635,659,683,707) dayTimeNum<-dayTimeIndex[day_num] #extract environmental data PP_carbon <- ncvar_get(nc, varid = "Phytoplankton", start = c(1,1,dayTimeNum), count = c(-1,-1,1)) #Node, siglay, time ZP_carbon <- ncvar_get(nc, varid = "Zooplankton", start = c(1,1,dayTimeNum), count = c(-1,-1,1)) Detritus <- ncvar_get(nc, varid = "Detritus", start = c(1,1,dayTimeNum), count = c(-1,-1,1)) temp <- ncvar_get(nc, varid = "temp", start = c(1,1,dayTimeNum), count = c(-1, -1, 1)) nc_close(nc) #set temp == 0 to NA for earlier winter months to identify not available habitat in those months. #Don't do this for December. if(month_num < 11){ # I don't expect ice onset for any scenario to be earlier than november temp[which(temp == 0)]<-NA # Areas with Ice are not possible habitat early in year thus the final Growth output should be NA } ##If weights in previous habitats were NA due to ice cover ##then change those start weights to the mean seen across that water column (nodal mean or rowMean) ##Rationale: If surface cell is frozen then the fish is at the bottom experiencing weight change in accordance with the habitat quality at the bottom ##, when the surface cell thaws then that fish has access to the water surface cell. So the starting weight for the surface cell ## should be based on the weight of the fish that uses that water column. Ice cover only occurs in shallow areas, so water column mean ## should be a reasonable value to use ## if whole water column is frozen then assign a start weight to the once frozen cell by using Inverse Distance Weighting of non-frozen cells nearby if((any(is.na(startWeight)))){ NAnodes<-coordsn[which(is.na((rowMeans(startWeight)))),]#c("node","lon", "lat")] #loop through nodes that have at least one NA and assign it a startWeight based on the mean end weight of its wc or mean end weight of nearby cells NAnodes<-as.data.frame(NAnodes) for(i in 1:nrow(NAnodes)){ #loop through each frozen node frznNode<-NAnodes[i,] #NAnodes[i,c("node", "lon", "lat")] frznNode$node<-as.numeric(frznNode[,c("node")]) #if(startWeight[which(is.nan(rowMeans(startWeight, na.rm = T)))]){ if(is.nan(mean(startWeight[frznNode$node,], na.rm=TRUE)) == FALSE){ #if the frozen cell is in a water column that is not completely frozen over than assign the water column mean endweight as the #frozen cell's start weight startWeight[frznNode$node,which(is.na(startWeight[frznNode$node,]))]<- mean(startWeight[frznNode$node,], na.rm = T) } else if(is.nan(mean(startWeight[frznNode$node,], na.rm=TRUE)) == TRUE){ #if frozen cell is in a water column that is completely frozen over, #assign it a start weight by averaging nearby non-frozen nodes #neighborCells = search neighborhood neighborCells<- coordsn[which(coordsn$lon <= (frznNode$lon + .02) & coordsn$lon >= (frznNode$lon - .02) & coordsn$lat <= (frznNode$lat + .05) & coordsn$lat >= (frznNode$lat - .05)),c("x","y","node","lon", "lat")] neighborCells<-neighborCells[-c(which(neighborCells$node == frznNode$node)),] # removes frzn node from search neighborhood # calculate distance (m) between frozen cell and neighbor cells and assign start weight using inverse distance weighting nearnodes<-data.frame(matrix(nrow = 0, ncol = 9)) for(q in 1:nrow(neighborCells)){ #loop through each neighboring node of each frozen node nearn<-neighborCells[q,] nearn$weight<-mean(startWeight[nearn$node,], na.rm= TRUE) #average water column weight of nearby nodes x1<-nearn$x y1<-nearn$y x2<-frznNode$x y2<-frznNode$y nearn$dist<-sqrt((x1-x2)^2 + (y1-y2)^2) #calculate non-frozen nodes distance from frozen node nearn$wtd<-nearn$weight/nearn$dist ##wtd == distance weighted contribution nearn$InvD<- (1/nearn$dist) ##inverse distance to be used in IDW interpolated value colnames(nearnodes)<-colnames(nearn) nearnodes<-rbind(nearnodes, nearn) rm(nearn) #neighborCells<-cbind(neighborCells,nearnodes$dist,nearnodes$wtd) } IDW_weight<-sum(nearnodes$wtd, na.rm = T)/sum(nearnodes$InvD, na.rm = T) #if other nearby nodes haven't are also frozen and haven't been interpolated yet, then remove nas and only average non-frozen nearby cells startWeight[frznNode$node,]<-IDW_weight } } } ED_PP <- 2600 #From Anderson et al. 2015; ED of Green algae (J/g) ED_ZP <- 2512 #Cooke and Hill 2010 use 2512 J/g of WET biomass for ED_ZP (from Cummins and Wuycheck 1971) ED_Det<- 979*.13 #Value and Adjustment factor derived from Anderson et al. 2016 #Convert PP_carbon (mg/L) to wet phytoplankton biomass (g/L) PP <- (PP_carbon*36)/(1000) #1/.0278 = 36 #Convert to Chla To compare Chla levels to Andersons (ug/L) Chla<-PP_carbon*.036*1000 #Convert ZP_carbon (mgC/L) to Wet zooplankton biomass (g/L) ZP <- (ZP_carbon * 2.5/.1)/1000 # 0.40 is C mass:dry mass (1/0.40=2.5), Det<- (Detritus/1000)/.044/.16 #(g/L) #.044 from Ozersky et al.; .16 dry:wet from Anderson et al. 2016 # Shared parameters # Cons eqn 2 CA <- .369 # Anderson et al. 2017 CB <- -.225 #Anderson et al. 2015, 2017 (taken from Wang 1989) CQ <- 2.5 #Value for bighead and silver from Kolar et al. 2007 # Resp eqn 1 SDA <- 0.1 #Tilapia value from Cooke and Hill 2010 RTO <- 0 ACT <- 1 # Waste eqn 2 FA <- .212 FB <- -.222 FG <- .631 UA <- .031 UB <- .58 UG <- -.299 if (species == "bighead"){ CTO <- 26 #Bighead, from Kolar et al. 2007 CTM <- 38 #Bighead, Kolar et al. 2007 W <- startWeight#5480 #Bighead wt. from Anderson et al. 2015 ED_carp <- 3500 #J/g Anderson et al. 2015 RA <- .00528 #Bighead, from Cooke and Hill 2010 RB <- -.299 #Bighead, from Cooke and Hill 2010 RQ <- .048 } if (species == "silver"){ CTO <- 29 #Silver, from Kolar et al. 2007 CTM <- 43 #Silver, Kolar et al. 2007 W <- startWeight#4350 #Silver wt. from Anderson et al. 2015 ED_carp <- 5200 #J/g Silver Anderson et al. 2015 RA <- .00279 #Silver, from Cooke and Hill 2010 RB <- -.239 #Silver, from Cooke and Hill 2010 RQ <- .076 #Silver, from Cooke and Hill 2010 } #Temperature function (f_t) V <- (CTM - temp)/(CTM - CTO) Z <- log(CQ) * (CTM-CTO) Y <- log(CQ) * (CTM-CTO+2) X <- ((Z^2) * (1+(1+40/Y)^0.5)^2)/400 # Calculate Cmax and f_t Cmax <- CA*W^CB f_t <- (V^X)*exp(X*(1-V)) # Filtration rate equation from Smith 1989 FR <- 1.54 * (W^0.713) * 24 # Calculate energy density ratio between carp and prey for ZP+PP diet and ZP+PP+Det diet ED_prey <- (ED_PP * PP) + (ED_ZP * ZP) + (ED_Det*Det) # J/L ED_prey2<- (ED_PP * PP) + (ED_ZP *ZP) #J/L Prey_ed <-ED_prey/(PP+ZP+Det) #Gives wtd. average j/g of prey Prey_ed2<-ED_prey2/(PP+ZP) Prey_ed[is.nan(Prey_ed)]<-0 #Set all NaNs to 0 Prey_ed2[is.nan(Prey_ed2)]<-0 ED_ratio<- Prey_ed/ED_carp #Converts gPrey to gCarp ED_ratio2<- Prey_ed2/ED_carp # Calculate Maximum Consumption rate (gPrey/gCarp/day) Cons <- Cmax * f_t #Calculate foraging model consumption rate (gPrey/gCarp/day) C <- (FR *(PP+ZP+Det)/W)*f_t #w/Det C2<- (FR *(PP+ZP)/W)*f_t #w/o Det # Pairwise min of each consumption output, assign minimum to Cnew and multiply by ED_ratio to convert gPrey/gCarp/day to gCarp/gCarp/day Cnew1 <- pmin(C, Cons)*ED_ratio #w/ det Cnew2<- pmin(C2, Cons)*ED_ratio2 #w/o det # Calculate proportion of maximum consumption (P) P1 <- pmin(C, Cons)/Cons #w/ det P1[is.nan(P1)]<-0 P2<- pmin(C2, Cons)/Cons #w/o det P2[is.nan(P2)]<-0 #Take pairwise MAX of Cnew1, Cnew2. Assign higher value to Cnew Cnew<- pmax(Cnew1, Cnew2) #Define function that allows fish to supplement their ZP+PP diet with detritus det.ration<-function(Cnew){ #Subtract foraging-based consumption rate from maximum cons to determine how much det the carp needs to eat to maximize consumption Det_ration <- Cons-C2 #gPrey/gcarp/d, amount of detritus required to reach maximum Cons #Calculate consumption based on total available detritus in the habitat (gDet/gCarp/d) C_det<- ((FR*Det)/W)*f_t #Take pmin; the amount of det eaten can never be greater than Det_ration CnewDet<-pmin(C_det, Det_ration) *(ED_Det/ED_carp) #puts CnewDet in gCarp/gCarp/d #Add CnewDet to Cnew Cnew<- Cnew+CnewDet return(Cnew) } #If Cnew2 (w/o det) defines Cnew, and Cnew is less than Cons. Then add det mass of difference b/w Cnew and Cons Cnew<-ifelse(Cnew == Cnew2 & Cnew < Cons, det.ration(Cnew), Cnew) #determine P_det, then add to P2 to define P when Cnew == Cnew2 Det_ration <- Cons-C2 C_det<- ((FR*Det)/W)*f_t P_det<-pmin(C_det, Det_ration)/Cons P2<- P2 + P_det # Define P P<- ifelse(Cnew == Cnew1, P1, P2) # Calculate temperature function for respiration Resp_f_T<- exp(RQ*temp) # Calculate respiration (g/g/day) R <- RA * (W^RB) * Resp_f_T * ACT # Calculate Egestion F1 <- FA * (temp^FB)*(exp(FG*P) * Cnew) # Calculate Excretion U <- UA * (temp^UB)*(exp(UG*P)*(Cnew-F1)) # Calculate SDA S <- SDA * (Cnew-F1) # Calculate growth (g/g/day) G <- (Cnew-(R+S+F1+U)) # Multiply growth rate by weight and add to initial weight to get final weight final_wt<- ((G*W)+W) #needs to be a matrix in order for it to work in yearGrowth.fun() # If final_wt is NA because temp in that cell = 0, then assign the last viable startWeight (non-NA value) as the final weight. #The habitat is not unavailable due to ice cover and the final weight should not change until it thaws final_wt[which(is.na(final_wt))]<-startWeight[which(is.na(final_wt))] return(final_wt) } ## Function that loops through specified time period (startMonth thru endMonth) on a daily time step and outputs final weight after that time period yearGrowth.fun<-function(species, startMonth,endMonth, file_list){ mo<- c(startMonth:endMonth) for(i in seq_along(ncfiles[startMonth:endMonth])){ #loop through months month_num<-mo[i] days<-c(1:30) for(q in seq_along(days)){ #nested loop that loops through days day_num<-days[q] if(i == 1 & day_num == 1){ #ensures that only the first day of the first month will start at initial weights if(species == "bighead"){ # looping through the 1st month the start weight should be 5480 or 4350 startWeight<- 5480 } if(species == "silver"){ startWeight<- 4350 } } else{ startWeight<-endWeight } endWeight<-Annualgrowth.fun(species, startWeight, month_num, day_num, ncfiles) } } return(endWeight)# returns endWeight at the final month (endMonth) } ### General scripts for Analyses ## Function that selects GRP values within specified depth_range ("NS" = 0-10m, "DCL"= 10-50m, "WC", "S") and averages across that depth range growth.depths<-function(GRP_df, depth_range){ #Load grid data and spatial plots load("E:/Polygons_Grid/shoreline_latlon.Rdata") load("E:/Polygons_Grid/Mich03_grid_polys.Rdata") load("E:/Polygons_Grid//Mich03_grid_coords.Rdata") load("E:/Polygons_Grid/Mich03_grid_coords_gb.Rdata")# green bay nodes ncfile <- "E:/FVCOM-GEM outputs/1M_L_2010(baseline)/mi_0008.nc" #Month doesn't matter here, only grabbing grid dimension vars nc <- nc_open(ncfile, readunlim = FALSE) siglay <- ncvar_get(nc, "siglay", start=c(1,1), count=c(1,1)) h <- ncvar_get(nc, "h", start=c(1), count=c(-1)) #meters siglev <- ncvar_get(nc, "siglev", start=c(1,1), count=c(1,-1)) nc_close(nc) rm(ncfile) #Make siglev and h into conformable arrays siglevr<-t(as.data.frame(siglev)) siglevr<-siglevr[rep(row.names(siglevr),5795),] siglevr<- drop(siglevr[,2:21]) #removes the surface level # store depth of each node in a data frame df.depths<- as.data.frame(h) df.depths<- as.data.frame(rep(df.depths,20)) #calculate depth of each cell (in a water column) depths<- -siglevr*df.depths colnames(depths)<- c(1:20) #print(max(depths)) #View(depths) if(depth_range == "S"){ depth_range_max = 1} #surface, 1m if(depth_range == "NS"){ depth_range_min = 0 depth_range_max = 10 #Near Surface 0-10m (similar to range of depths of BHC in IR) or 0-20 depending on DCL depth_range_maxB = max(depths[which(depths[,1] >10),1])} # for nodes where 1st siglay is >10m deep. 13.5585==max(depths[which(depths[,1] >10),1]) if(depth_range == "DCL"){ depth_range_min = 10 #20-60 m is historical depth (Fahnenstiel & Scavia 1987) depth_range_minB = max(depths[which(depths[,1] >10),1]) ## For nodes where 1st siglay is >10m deep. This ensures that the 1st siglay doesn't get counted twice in NS and DCL evals depth_range_max = 50} #10-50 m is post quagga depth (Bramburger & Reavie 2016) if(depth_range == "WC"){ depth_range_min = 0 depth_range_max = 272} #whole wc 0-272 m ## Identify nodes where 1st siglay is deeper than 10 m deepNodes<-coordsn[which(depths[,1] >10),c("node")] #Create empty df to store new values in GRPdepth_df<- data.frame(matrix(nrow= 5795, ncol = 20)) #loop through the whole dataset and select GRP value from cells within depth range for(i in 1:nrow(depths)){ if(depth_range == "S"){ if(depths[i,1]>depth_range_max){ #for surface evals, if surface layer is deeper than 1m; or for NS evals, if 1st siglay is deeper than 10m GRPdepth_df[i,1]<-GRP_df[i,1] #then assign only the first layer to output df for that node }else{ i2<-1 for(i2 in 1:ncol(depths)){ if(depths[i,i2]<=depth_range_max){ GRPdepth_df[i,i2]<-GRP_df[i,i2] }else{ GRPdepth_df[i,i2]<- NA } i2+1 } i+1 } } # only for Surface evals else{ if(depths[i,20]<=depth_range_max & depths[i,1]>= depth_range_min){ GRPdepth_df[i,]<-GRP_df[i,] } else{ i2<-1 for(i2 in 1:ncol(depths)){ if(any(deepNodes == i) & i2 == 1){ ## Use different depths for nodes where 1st siglay >10m if(depth_range == "NS"){ if(depths[i,1]<=depth_range_maxB & depths[i,1]>= depth_range_min){ GRPdepth_df[i,1]<-GRP_df[i,i2] } } if(depth_range == "DCL"){ if(depths[i,1]<=depth_range_max & depths[i,1]>= depth_range_minB){ GRPdepth_df[i,1]<-GRP_df[i,1] } } } else{ if(depths[i,i2]<=depth_range_max & depths[i,i2]>= depth_range_min){ GRPdepth_df[i,i2]<-GRP_df[i,i2] } else{ GRPdepth_df[i,i2]<- NA } } i2+1 } i+1 } #depths other than surface } }#end of loop # Append nodal mean GRPdepth_df$meanG<-rowMeans(GRPdepth_df,na.rm = TRUE) GRPdepth_df$Gmax<-apply(GRPdepth_df, 1, function(n) max(n,na.rm = TRUE)) #this produces warnings but they are harmless return(GRPdepth_df) } ## Calculate suitable habitat extent in surface area (km2) and volume (km3) given a df of GRP suit.hab = function(df){ ncfile <- ncfiles[5] #Month doesn't matter here, only grabbing grid dimension vars nc <- nc_open(ncfile, readunlim = FALSE) siglay <- ncvar_get(nc, "siglay", start=c(1,1), count=c(1,1)) h <- ncvar_get(nc, "h", start=c(1), count=c(-1)) #meters siglev <- ncvar_get(nc, "siglev", start=c(1,1), count=c(1,-1)) nc_close(nc) #Calculate total volume of the lake, use coordsn$depth it is more accurate than 'h' art1 <- calcArea(psTces) #area of each TCE in km2 areas <- art1$area vol<- coordsn$depth*(areas*1000*1000) #mult by 1000 twice to put in m3 volumes<-cbind(art1,vol) totalArea <- sum(areas) totalVolume <- sum(vol)# total volume of lake in m3 #Make siglev and vol into conformable arrays siglevr<-t(as.data.frame(siglev)) #Put siglevr in terms of individual proportions instead of cumulative proportions for(i in 1:ncol(siglevr)){ if(i == 1){ siglevr [i] <- -siglev[i] } else if(i > 1){ siglevr[i] <- -siglev[i] - -siglev[i-1] } } #sum(siglevr) == 1 # Check to make sure each proportion adds up to 1 siglevr<-siglevr[rep(row.names(siglevr),5795),] siglevr<- drop(siglevr[,2:21]) #removes the surface level # store volumes of each node (whole water column volume) in a data frame df.vol<- as.data.frame(vol) df.vol<- as.data.frame(rep(df.vol,20)) # store surface areas of each node in a data frame df.areas<- as.data.frame(areas) #calculate vol of each cell (in a water column) cellvols<- siglevr*df.vol #each cells volume in m3 colnames(cellvols)<- c(1:20) if (sum(cellvols) != totalVolume){ break print("sum(cellvols != totalVolume)") } #check to make sure volume is correct # Append node grp maxima to the grp data frame dfmax<- df dfmax <- apply(dfmax, 1, function(n) max(n,na.rm = TRUE)) dfmax<- as.data.frame(dfmax) # If no nodes are suitable than store 0's for area, vol, and meangrp = NA if(sum(dfmax>0)==0){ #suithab_df <- matrix(nrow = 1, ncol = 7) suithab_df <- data.frame(matrix(nrow = 1, ncol = 0)) suithab_df$Area <- 0 suithab_df$Volume<- 0 suithab_df$MeanGRP<- NA suithab_df$MinGRP<-NA suithab_df$GRP_Q1<-NA suithab_df$MedGRP<-NA suithab_df$GRP_Q3<-NA suithab_df$MaxGRP<-NA #suithab_df<-suithab_df[,4:6] return(suithab_df) } else{ #Define function to only select nodes with at least one positive growth cell suit.nodes<-function(df){ suitnodes<-data.frame(matrix(nrow = sum(dfmax>=0), ncol = 20)) for(i in 1:nrow(dfmax)){ if(any(dfmax[i,]>=0)){ suitnodes[i,]<-(df[i,]) } i+1} suitnodes } suitnodes<-suit.nodes(df) # remove rows w/o any positive growth suitcells<-suitnodes[complete.cases(suitnodes),] colnames(suitcells)<-c(1:20) # Select nodes with at least one positive sigma layer in their water column suitareas <- as.data.frame(df.areas[rownames(suitcells),1]) # Make negative cells == NA suitcells[suitcells < 0] <- NA # Select volumes for all cells in positive nodes only suitvols<- cellvols[rownames(suitcells), ] #Set any negative grp cell volumes and areas to NA suitvols[is.na(suitcells[,1:20])]<-NA #Sum total volume of suitable habitat, m3 SuitHabVol<- sum(suitvols, na.rm = TRUE) SuitHabVol_km3<- SuitHabVol/(1000*1000*1000) #put in km3 #Sum total surface area of suitable habitat SuitHabArea_km2<- sum(suitareas, na.rm = TRUE) #km2 #Mean growth rate in suitable habitat cells meangrowth<- sum(rowSums((suitcells*(suitvols/SuitHabVol)), na.rm = T)) #Cell Volume Weighted GRP average #Min, Q1, Median, Q3, and Max growth rates in suitable cells q_suit<-quantile(suitcells, names = FALSE, na.rm = TRUE) MinGRP<-q_suit[1] GRP_Q1<-q_suit[2] MedGRP<-q_suit[3] GRP_Q3<-q_suit[4] MaxGRP<-q_suit[5] #suithab_df <- matrix(nrow = 1, ncol = 0) suithab_df <- data.frame(matrix(nrow = 1, ncol = 0)) suithab_df$Area <- SuitHabArea_km2 suithab_df$Volume<- SuitHabVol_km3 suithab_df$MeanGRP<- meangrowth suithab_df$MinGRP<-MinGRP suithab_df$GRP_Q1<-GRP_Q1 suithab_df$MedGRP<-MedGRP suithab_df$GRP_Q3<-GRP_Q3 suithab_df$MaxGRP<-MaxGRP #suithab_df<-suithab_df[,4:6] return(suithab_df) } #end of if statement }