rm(list=ls()) ###Using theatmodel function for the ectotherm biophysical model coded by Rubalcaba et al. 2019 Global Ecology and Biogeography theatmodel <- function(Tb, # Body temperature at time t (the function gives body temperature at time t+delta) (ºC) A, # Total surface area of skin (m^2) M, # Body mass (g) Ta, # Air temperature (ºC) Tg, # Ground temperature (ºC) S, # Solar radiation (W m^-2) v, # Wind speed (m s^-1) HR, # Relative humidity (%) # SWP, # Soil water potential (millibars) # HS, #Hydration state (% of max. hydration) skin_humidity = 1, # Skin humidity (0-1) r, # Total skin resitance to water loss (s m^-1) posture = 1, # postural adjustment: proportion of dorsal skin surface exposed to the sunlight vent = 1/3, # proportion of ventral A delta = 60, # Number of time steps (seconds) C = 3.6 # Specific heat capacity (J g-1 ºC-1) ) { for(i in 1:delta){ ## Heat transfer function Q <- function(Tb){ qabs <- A1 * a * S + 0.2 * A1 * S + 0.2 * A1 * S # Absorption of solar radiation qevap <- L * EWL_estimate # Evaporative cooling qrad <- A1 * sigma * epsilon * (Tb^4 - Ta^4) # Emission of long-wave radiation qconv <- A1 * hc * (Tb - Ta) # COnvection qcond <- Ag * hg * (Tb - Tg) # Conduction to the ground qnet <- qabs - qevap - qrad - qconv - qcond # Net heat flux return(qnet) } ## Model parameters # Animal geometry and color Ag = vent * A # ventral area (m^2) A1 = (1-vent) * A * posture # exposed area (m^2) l = sqrt(A) # caracterstic length (m) C = C # average heat capacity body (J g^-1 ºC^-1) epsilon = 0.9 # Emissivity to long-wave radiation a = 0.8 # Absorptance to short-wave solar radiation # Constants L = 2257 # latent heat of vaporization of water (J g^-1) sigma = 5.670367e-8 # Stefan-Boltzmann constant (W m^-2 ºC^-4) # Estimation convection coefficient rho = 101325 / (287.04 * (Ta + 273)) nu = -1.1555e-14*(Ta+273)^3 + 9.5728e-11*(Ta+273)^2 + 3.7604e-08*(Ta+273) - 3.4484e-06 # kinematic viscosity (m^2 s^-1) kf = 1.5207e-11*(Ta+273)^3 - 4.8574e-08*(Ta+273)^2 + 1.0184e-04*(Ta+273) - 3.9333e-04 # thermal conductivity air (W m^-1 ºC^-1) Re = l * v / nu # Reynolds number c = 0.2 # Constants relating Re and Nu numbers n = 0.7 Nu = 2 + c * Re^n # Nusselt number hc = Nu * kf / l # Convection coefficient (W m-2 ºC-1) # Estimation of conduction coefficient ksub = 0.027 # thermal conductivity of substrate (W ºC-1) ts = 0.025 * (0.001 * M / (pi * 1000))^0.2 # thinkness of surface layer (Kleiber 1972) hg = ksub / ts # Conduction coeficient (W m-2 ºC-1) (Stevenson 1985) # Estimate evaporative water loss (g s-1) HR_prop = HR * 0.01 ps_a = exp(77.3450 + 0.0057 * (Ta+273) - 7235 / (Ta+273)) / (Ta+273)^8.2 # Air vapor pressure (Pa) rho_a = HR_prop * 2.2 * ps_a / (Ta+273) # Trasform into density (g m-3) ps_b = exp(77.3450 + 0.0057 * (Tb+273) - 7235 / (Tb+273)) / (Tb+273)^8.2 # Body vapor density (Pa) rho_b = skin_humidity * 2.2 * ps_b / (Tb+273) # Trasform into density (g m-3) EWL_estimate = A1 * (rho_b - rho_a) / r # Evaporative water loss (g s-1) #Estimate water flux to soil (g s-1) #Constant, hydraulic conductance of frog (g cm-2 s-1 mb-1) # K_fr = (1.3 * 10^-5)/60 #for a frog from 12-28 C # FrWP = -2.365e+05 + 7.788e+03*HS-8.609e+01*HS^2+3.185e-01*HS^3#Frog water potential, inferred from the hydration state of the frog (Tracy 1976) # K_soil = 10^-5 # D_soil = 10^-2 # s = 1 + ((K_fr*(Ag/(2*pi))^0.5)/K_soil)*(D_soil/(Ag/(2*pi)))^0.5 # M_eq = (SWP-FrWP)*(1/(Ag*K_fr)+1/((2*pi*Ag)^0.5*K_soil))) # M = M_eq*(1+(K_fr*(Ag/(2*pi)^0.5)*exp(s^2)*erfc(s))/K_soil) # M_soil = Ag*K_fr*(FrWP - 0) #Water flux (g s-1) # HS = HS + M_soil/M - EWL_estimate/M ## Model iteration # Body Tempearature at t + 1" Tb = Tb + 1 / (C * M) * Q(Tb) } # Output: return body temperature at time t + delta (ºC) and instantaneous rate of evaporative water loss (g s-1) return(c(Tb, EWL_estimate)) } #Functons for geometric mean & standard deviation gm_mean<- function(x){ exp(mean(log(x))) } ####1. Species characteristics #Adult snout-vent lengths for each species #Ascaphus truei (Matsuda et al 2006) SVL_at=c(30,40,50) #min. and max. reported #Pseudacris regilla (Matsuda et al 2006) SVL_pr=c(25,37.5,50) #min. and max. reported #Spea intermontana (Nussbaum et al. 1983; Matsuda et al 2006) SVL_si=c(40,52.5,65) #min. and max. reported ##Estimates of body mass for each species from SVL (Santini et al. 2018) #Ascaphus truei- Hylid curve M_at<-10^-4.462 * SVL_at^3.201 #Pseudacris regilla -Hylid curve M_pr<-10^-4.462 * SVL_pr^3.201 #Spea intermontana -Terrestrial anuran curve M_si<- 10^-4.298 * SVL_si^3.181 ##Estimates of skin surface area (Klein et al. 2016) A_at<- 9.1411 * M_at^(0.7735) * 1e-4 #Hylid curve A_pr<- 9.1411 * M_pr^(0.7735) * 1e-4 #Hylid curve A_si<- 10.1158 * M_si^(0.6012) * 1e-4 #Scaphiopod curve #For all species: r_at = 300 # Skin resistance to water loss (s m^-1), just boundary layer (Claussen 1972) r_pr = 500 # Boundary layer (3 s m^-1) + skin resistance (2 s m^-1) (Withers et al. 1982) r_sp = 800 # Boundary layer (5 s m^-1) + skin resistance (2 s m^-1) (Wygoda 1984) C = 3.6 # Specific heat capacity of the body (J g-1) Morph_at<- data.frame(M_at,A_at) Morph_pr<- data.frame(M_pr,A_pr) Morph_si<- data.frame(M_si,A_si) ####Set thermal limits - for temperature hours of restriction ###Ascaphus truei #Upper 80% thermal from Murray et al. (unpub.) t80_at<- 23.9 #Pseudacris regilla t80_pr<- 28.9 #Upper 80% thermal from Gerick et al. 2014 #Spea intermontana t80_sp<- 32.2 #Upper 80% thermal from Gerick et al. 2014 ##Define models and fitness functions for each species library(nlme) hyd_dat<- read.csv("jump_performance_trials.csv") hyd_dat_a<- subset(hyd_dat,Species=='A_truei') hyd_dat_a$W<- 1-hyd_dat_a$hydration hyd_dat_a$W.2<- hyd_dat_a$W^2 hyd_dat_a$W.3<- hyd_dat_a$W^3 hyd_dat_a$y<- hyd_dat_a$jump_dist quadratic_a<- lme(y~scale(W)+scale(W.2),random=~1|id,data=hyd_dat_a,method='REML') cubic_a_t3<- lme(y~scale(W)+scale(W.2)+scale(W.3)*scale(temp),random=~1|id,data=hyd_dat_a,method='REML') theoret_a<- expand.grid(temp=c(15,20,24),W=seq(0,0.3,length.out = 100)) theoret_a$W.2<- theoret_a$W^2 theoret_a$W.3<- theoret_a$W^3 theoret_a$pred.h<- predict(quadratic_a,newdata=theoret_a,level=0) theoret_a$pred.ht<- predict(cubic_a_t3,newdata=theoret_a,level=0) ##Bring in species locality data require(raster) a_loc<- shapefile("ascaphus_truei_gbif") q<- quantile(a_loc$X2,prob=seq(0,1,0.1)) # percentiles q<- q[-11] # function to retrieve closest quantile for a given value. get.quantile <- function(x)names(q)[which(abs(q-x)==min(abs(q-x)))] # apply this function for all values in df$Amount_Spent a_loc$Quantile<- sapply(a_loc$X2,get.quantile) a_loc<- as.data.frame(a_loc) #Subset 10 locations from each decile of the species' latitudinal distribution sub_locs<- list() for(i in 1:length(q)){ a_q<- subset(a_loc,a_loc$Quantile==names(q[i])) sub_locs[[i]]<- a_q[sample(nrow(a_q),10),] } sub_locs_df<- do.call("rbind", sub_locs) require(NicheMapR) ###Full sun microclimate data_sun<- list() month<- c(7) #For the month of july #get.global.climate(folder = 'global climate') for(i in 1:nrow(sub_locs_df)){ tryCatch( micro <- micro_global(loc = sub_locs_df[i,1:2], minshade = 50, maxshade = 100,runmoist = 1,runshade=0), error=function(e){cat("ERROR :",conditionMessage(e), "\n")}) metout <- as.data.frame(micro$metout) soil <- as.data.frame(micro$soil) soilp <- as.data.frame(micro$soilpot) met_month<- list() met_month_min<- list() for(z in 1:length(month)){ day = unique(metout$DOY)[month[z]] DAY <- day TIME <- subset(metout, DOY == day)$TIME RAD <- subset(metout, DOY == day)$SOLR AIRT <- subset(metout, DOY == day)$TALOC SOILT <- subset(soil, DOY == day)$D0cm RH <- subset(metout, DOY == day)$RHLOC V <- subset(metout, DOY == day)$VLOC SOILPOT<- subset(soilp, DOY == day)$PT0cm*10 met_month<- data.frame(DAY,TIME, RAD, AIRT, SOILT, RH, V, SOILPOT) modif <- data.frame(DAY=rep(DAY,1440),TIME=1:1440, RAD=NA, AIRT=NA, SOILT=NA, RH=NA, V=NA, SOILPOT=NA) for(j in 3:8){ f <- approxfun(met_month[,2], met_month[,j]) S <- f(seq(0, 1380, length = 1440)) modif[,j] <- S } met_month_min[[z]]<- modif } data_sun[[i]]<- do.call("rbind", met_month_min) } iter=1440 limits_data_frame_at<- list() for(k in 1:length(data_sun)){ # Body temperature and EWL for each month limits_data_frame_at[[k]]<- data.frame(X=NA,Y=NA,month=NA,h_t_1=NA,h_h_1=NA,h_ht_1=NA,h_t_2=NA,h_h_2=NA,h_ht_2=NA,h_t_3=NA,h_h_3=NA,h_ht_3=NA) for(q in 1:length(unique(data_sun[[k]]$DAY))){ limits_data_frame_at[[k]][q,1]<- sub_locs_df$X1[k] limits_data_frame_at[[k]][q,2]<- sub_locs_df$X2[k] limits_data_frame_at[[k]][q,3]<- unique(data_sun[[k]]$DAY)[q] data<- subset(data_sun[[k]],DAY==unique(data_sun[[k]]$DAY)[q]) A1 = Morph_at$A_at[1] M1 = Morph_at$M_at[1] A2 = Morph_at$A_at[2] M2 = Morph_at$M_at[2] A3 = Morph_at$A_at[3] M3 = Morph_at$M_at[3] Tb1<- rep(data$AIRT[1],1440) EWL1<- rep(0,1440) Tb2<- rep(data$AIRT[1],1440) EWL2<- rep(0,1440) Tb3<- rep(data$AIRT[1],1440) EWL3<- rep(0,1440) for(j in 2:length(Tb1)){ sun<- data$RAD[j] Ta<- data$AIRT[j] Tg<- data$SOILT[j] RH<- data$RH[j] V<- data$V[j] m_6_tb_m1<- theatmodel(Tb1[j-1], A1, M1, Ta, Tg, sun, V, RH, skin_humidity = 1, r = r_at, posture = 1, delta = 60) Tb1[j]<- m_6_tb_m1[1] EWL1[j]<- m_6_tb_m1[2] m_6_tb_m2<- theatmodel(Tb2[j-1], A2, M2, Ta, Tg, sun, V, RH, skin_humidity = 1, r = r_at, posture = 1, delta = 60) Tb2[j]<- m_6_tb_m2[1] EWL2[j]<- m_6_tb_m2[2] m_6_tb_m3<- theatmodel(Tb3[j-1], A3, M3, Ta, Tg, sun, V, RH, skin_humidity = 1, r = r_at, posture = 1, delta = 60) Tb3[j]<- m_6_tb_m3[1] EWL3[j]<- m_6_tb_m3[2] } REWL1<- EWL1/M1*60 REWL2<- EWL2/M2*60 REWL3<- EWL3/M3*60 ###Small freg df1<- data.frame(temp=Tb1,REWL1,data$TIME) a1<- subset(df1,Tb124){ temp[y]=24 } } dat<- data.frame(W=w,W.2=w^2,W.3=w^3,temp=temp) pred_h<- predict(quadratic_a,newdata=dat,level=0)/max(theoret_a$pred.h) pred_h_max<- which.max(pred_h) pred_h_sub<- pred_ht[pred_h_max:length(pred_h)] pred_ht<- predict(cubic_a_t3,newdata=dat,level=0)/max(theoret_a$pred.ht) pred_ht_max<- which.max(pred_ht) pred_ht_sub<- pred_ht[pred_ht_max:length(pred_ht)] length_h1[i]<-length(which(pred_h_sub>0.8))+pred_h_max-1 length_ht1[i]<-length(which(pred_ht_sub>0.8))+pred_ht_max-1 } ##medium freg df2<- data.frame(temp=Tb2,REWL2,data$TIME) a2<- subset(df2,Tb224){ temp[y]=24 } } dat<- data.frame(W=w,W.2=w^2,W.3=w^3,temp=temp) pred_h<- predict(quadratic_a,newdata=dat,level=0)/max(theoret_a$pred.h) pred_h_max<- which.max(pred_h) pred_h_sub<- pred_ht[pred_h_max:length(pred_h)] pred_ht<- predict(cubic_a_t3,newdata=dat,level=0)/max(theoret_a$pred.ht) pred_ht_max<- which.max(pred_ht) pred_ht_sub<- pred_ht[pred_ht_max:length(pred_ht)] length_h2[i]<-length(which(pred_h_sub>0.8))+pred_h_max-1 length_ht2[i]<-length(which(pred_ht_sub>0.8))+pred_ht_max-1 } ##big froggo df3<- data.frame(temp=Tb3,REWL3,data$TIME) a3<- subset(df3,Tb324){ temp[y]=24 } } dat<- data.frame(W=w,W.2=w^2,W.3=w^3,temp=temp) pred_h<- predict(quadratic_a,newdata=dat,level=0)/max(theoret_a$pred.h) pred_h_max<- which.max(pred_h) pred_h_sub<- pred_ht[pred_h_max:length(pred_h)] pred_ht<- predict(cubic_a_t3,newdata=dat,level=0)/max(theoret_a$pred.ht) pred_ht_max<- which.max(pred_ht) pred_ht_sub<- pred_ht[pred_ht_max:length(pred_ht)] length_h3[i]<-length(which(pred_h_sub>0.8))+pred_h_max-1 length_ht3[i]<-length(which(pred_ht_sub>0.8))+pred_ht_max-1 } limits_data_frame_at[[k]][q,4]<- nrow(a1) limits_data_frame_at[[k]][q,5]<- median(length_ht1) limits_data_frame_at[[k]][q,6]<- gm_mean(length_ht1) limits_data_frame_at[[k]][q,7]<- nrow(a2) limits_data_frame_at[[k]][q,8]<- median(length_ht2) limits_data_frame_at[[k]][q,9]<- gm_mean(length_ht2) limits_data_frame_at[[k]][q,10]<- nrow(a3) limits_data_frame_at[[k]][q,11]<- median(length_ht3) limits_data_frame_at[[k]][q,12]<- gm_mean(length_ht3) } } at_limits_summary<- do.call("rbind", limits_data_frame_at) write.csv(at_limits_summary,'Ascaphus_hours_restriction.csv') ###P. regilla ##Define models and fitness functions for each species hyd_dat_p<- subset(hyd_dat,Species=="P_regilla") hyd_dat_p$W<- 1-hyd_dat_p$hydration hyd_dat_p$W.2<- hyd_dat_p$W^2 hyd_dat_p$W.3<- hyd_dat_p$W^3 hyd_dat_p$y<- hyd_dat_p$jump_dist cubic_p<- lme(y~scale(W)+scale(W.2)+scale(W.3),random=~1|id,data=hyd_dat_p,method='REML') cubic_p_t3<- lme(y~scale(W)+scale(W.2)+scale(W.3)*scale(temp),random=~1|id,data=hyd_dat_p,method='REML') theoret_p<- expand.grid(temp=c(20,24,27,30),W=seq(0,0.3,length.out = 100)) theoret_p$W.2<- theoret_p$W^2 theoret_p$W.3<- theoret_p$W^3 theoret_p$pred.h<- predict(cubic_p,newdata=theoret_p,level=0) theoret_p$pred.ht<- predict(cubic_p_t3,newdata=theoret_p,level=0) ##Bring in species locality data require(raster) pr_loc<- shapefile("pseudacris_regilla") q<- quantile(pr_loc$X2,prob=seq(0,0.9,0.1)) # percentiles # function to retrieve closest quantile for a given value. get.quantile <- function(x)names(q)[which(abs(q-x)==min(abs(q-x)))] # apply this function for all values in df$Amount_Spent pr_loc$Quantile<- sapply(pr_loc$X2,get.quantile) pr_loc<- as.data.frame(pr_loc) #Subset 10 locations from each decile of the species' latitudinal distribution sub_locs<- list() for(i in 1:length(q)){ pr_q<- subset(pr_loc,pr_loc$Quantile==names(q[i])) sub_locs[[i]]<- pr_q[sample(nrow(pr_q),10),] } sub_locs_df<- do.call("rbind", sub_locs) require(NicheMapR) ###Full sun microclimate data_sun<- list() month<- c(7) #For the months july #get.global.climate(folder = 'global climate') for(i in 1:nrow(sub_locs_df)){ tryCatch( micro <- micro_global(loc = sub_locs_df[i,1:2], minshade = 50, maxshade = 100,runmoist = 1,runshade=0), error=function(e){cat("ERROR :",conditionMessage(e), "\n")}) metout <- as.data.frame(micro$metout) soil <- as.data.frame(micro$soil) soilp <- as.data.frame(micro$soilpot) met_month<- list() met_month_min<- list() for(z in 1:length(month)){ day = unique(metout$DOY)[month[z]] DAY <- day TIME <- subset(metout, DOY == day)$TIME RAD <- subset(metout, DOY == day)$SOLR AIRT <- subset(metout, DOY == day)$TALOC SOILT <- subset(soil, DOY == day)$D0cm RH <- subset(metout, DOY == day)$RHLOC V <- subset(metout, DOY == day)$VLOC SOILPOT<- subset(soilp, DOY == day)$PT0cm*10 met_month<- data.frame(DAY,TIME, RAD, AIRT, SOILT, RH, V, SOILPOT) modif <- data.frame(DAY=rep(DAY,1440),TIME=1:1440, RAD=NA, AIRT=NA, SOILT=NA, RH=NA, V=NA, SOILPOT=NA) for(j in 3:8){ f <- approxfun(met_month[,2], met_month[,j]) S <- f(seq(0, 1380, length = 1440)) modif[,j] <- S } met_month_min[[z]]<- modif } data_sun[[i]]<- do.call("rbind", met_month_min) } iter=1440 limits_data_frame_pr<- list() for(k in 1:length(data_sun)){ # Body temperature and EWL for each month limits_data_frame_pr[[k]]<- data.frame(X=NA,Y=NA,month=NA,h_t_1=NA,h_ht_1=NA,m_h_ht_1=NA,sd.h_ht_1=NA,h_t_2=NA,h_ht_2=NA,m_h_ht_2=NA,sd.h_ht_2=NA,h_t_3=NA,h_ht_3=NA,m_h_ht_3=NA,sd.h_ht_3=NA) for(q in 1:length(unique(data_sun[[k]]$DAY))){ limits_data_frame_pr[[k]][q,1]<- sub_locs_df$X1[k] limits_data_frame_pr[[k]][q,2]<- sub_locs_df$X2[k] limits_data_frame_pr[[k]][q,3]<- unique(data_sun[[k]]$DAY)[q] data<- subset(data_sun[[k]],DAY==unique(data_sun[[k]]$DAY)[q]) A1 = Morph_pr$A_pr[1] M1 = Morph_pr$M_pr[1] A2 = Morph_pr$A_pr[2] M2 = Morph_pr$M_pr[2] A3 = Morph_pr$A_pr[3] M3 = Morph_pr$M_pr[3] Tb1<- rep(data$AIRT[1],1440) EWL1<- rep(0,1440) Tb2<- rep(data$AIRT[1],1440) EWL2<- rep(0,1440) Tb3<- rep(data$AIRT[1],1440) EWL3<- rep(0,1440) for(j in 2:length(Tb1)){ sun<- data$RAD[j] Ta<- data$AIRT[j] Tg<- data$SOILT[j] RH<- data$RH[j] V<- data$V[j] # SWP = data$SOILPOT[j] m_6_tb_m1<- theatmodel(Tb1[j-1], A1, M1, Ta, Tg, sun, V, RH, skin_humidity = 1, r = r_pr, posture = 1, delta = 60) Tb1[j]<- m_6_tb_m1[1] EWL1[j]<- m_6_tb_m1[2] m_6_tb_m2<- theatmodel(Tb2[j-1], A2, M2, Ta, Tg, sun, V, RH, skin_humidity = 1, r = r_pr, posture = 1, delta = 60) Tb2[j]<- m_6_tb_m2[1] EWL2[j]<- m_6_tb_m2[2] m_6_tb_m3<- theatmodel(Tb3[j-1], A3, M3, Ta, Tg, sun, V, RH, skin_humidity = 1, r = r_pr, posture = 1, delta = 60) Tb3[j]<- m_6_tb_m3[1] EWL3[j]<- m_6_tb_m3[2] } REWL1<- EWL1/M1*60 REWL2<- EWL2/M2*60 REWL3<- EWL3/M3*60 ###Small freg df1<- data.frame(temp=Tb1,REWL1,data$TIME) a1<- subset(df1,Tb130){ temp[y]=30 } } dat<- data.frame(W=w,W.2=w^2,W.3=w^3,temp=temp) pred_h<- predict(cubic_p,newdata=dat,level=0)/max(theoret_p$pred.h) pred_h_max<- which.max(pred_h) pred_h_sub<- pred_ht[pred_h_max:length(pred_h)] pred_ht<- predict(cubic_p_t3,newdata=dat,level=0)/max(theoret_p$pred.ht) pred_ht_max<- which.max(pred_ht) pred_ht_sub<- pred_ht[pred_ht_max:length(pred_ht)] length_h1[i]<-length(which(pred_h_sub>0.8))+pred_h_max-1 length_ht1[i]<-length(which(pred_ht_sub>0.8))+pred_ht_max-1 } ##medium freg df2<- data.frame(temp=Tb2,REWL2,data$TIME) a2<- subset(df2,Tb230){ temp[y]=30 } } dat<- data.frame(W=w,W.2=w^2,W.3=w^3,temp=temp) pred_h<- predict(cubic_p,newdata=dat,level=0)/max(theoret_p$pred.h) pred_h_max<- which.max(pred_h) pred_h_sub<- pred_h[pred_h_max:length(pred_h)] pred_ht<- predict(cubic_p_t3,newdata=dat,level=0)/max(theoret_p$pred.ht) pred_ht_max<- which.max(pred_ht) pred_ht_sub<- pred_ht[pred_ht_max:length(pred_ht)] length_h2[i]<-length(which(pred_h_sub>0.8))+pred_h_max-1 length_ht2[i]<-length(which(pred_ht_sub>0.8))+pred_ht_max-1 } ##big froggo df3<- data.frame(temp=Tb3,REWL3,data$TIME) a3<- subset(df3,Tb330){ temp[y]=30 } } dat<- data.frame(W=w,W.2=w^2,W.3=w^3,temp=temp) pred_h<- predict(cubic_p,newdata=dat,level=0)/max(theoret_p$pred.h) pred_h_max<- which.max(pred_h) pred_h_sub<- pred_h[pred_h_max:length(pred_h)] pred_ht<- predict(cubic_p_t3,newdata=dat,level=0)/max(theoret_p$pred.ht) pred_ht_max<- which.max(pred_ht) pred_ht_sub<- pred_ht[pred_ht_max:length(pred_ht)] length_h3[i]<-length(which(pred_h_sub>0.8))+pred_h_max-1 length_ht3[i]<-length(which(pred_ht_sub>0.8))+pred_ht_max-1 } limits_data_frame_pr[[k]][q,4]<- nrow(a1) limits_data_frame_pr[[k]][q,5]<- median(length_ht1) limits_data_frame_pr[[k]][q,6]<- gm_mean(length_ht1) limits_data_frame_pr[[k]][q,7]<- nrow(a2) limits_data_frame_pr[[k]][q,8]<- median(length_ht2) limits_data_frame_pr[[k]][q,9]<- gm_mean(length_ht2) limits_data_frame_pr[[k]][q,10]<- nrow(a3) limits_data_frame_pr[[k]][q,11]<- median(length_ht3) limits_data_frame_pr[[k]][q,12]<- gm_mean(length_ht3) } print(k) } pr_limits_summary<- do.call("rbind", limits_data_frame_pr) write.csv(pr_limits_summary,'Pseudacris_hours_restriction.csv') ###P. regilla ##Define models and fitness functions for each species hyd_dat_s$W<- 1-hyd_dat_s$hydration hyd_dat_s$W.2<- hyd_dat_s$W^2 hyd_dat_s$W.3<- hyd_dat_s$W^3 hyd_dat_s$y<- hyd_dat_s$jump_dist cubic_s<- lme(y~scale(W)+scale(W.2)+scale(W.3),random=~1|id,data=hyd_dat_s,method='REML') cubic_s_t2<- lme(y~scale(W)+scale(W.2)*scale(temp)+scale(W.3),random=~1|id,data=hyd_dat_s,method='REML') theoret_s<- expand.grid(temp=c(24,27,30),W=seq(0,0.4,length.out = 100)) theoret_s$W.2<- theoret_s$W^2 theoret_s$W.3<- theoret_s$W^3 theoret_s$pred.h<- predict(cubic_s,newdata=theoret_s,level=0) theoret_s$pred.ht<- predict(cubic_s_t2,newdata=theoret_s,level=0) ##Bring in species locality data sp_loc<- shapefile("spea_intermontana") q<- quantile(sp_loc$X2,prob=seq(0,0.9,0.1)) # percentiles # function to retrieve closest quantile for a given value. get.quantile <- function(x)names(q)[which(abs(q-x)==min(abs(q-x)))] # apply this function for all values in df$Amount_Spent sp_loc$Quantile<- sapply(sp_loc$X2,get.quantile) sp_loc<- as.data.frame(sp_loc) #Subset 10 locations from each decile of the species' latitudinal distribution sub_locs<- list() for(i in 1:length(q)){ sp_q<- subset(sp_loc,sp_loc$Quantile==names(q[i])) sub_locs[[i]]<- sp_q[sample(nrow(sp_q),10),] } sub_locs_df<- do.call("rbind", sub_locs) require(NicheMapR) ###Full sun microclimate data_sun<- list() month<- c(7) #For the months of july #get.global.climate(folder = 'global climate') for(i in 1:nrow(sub_locs_df)){ tryCatch( micro <- micro_global(loc = sub_locs_df[i,1:2], minshade = 50, maxshade = 100,runmoist = 1,runshade=0), error=function(e){cat("ERROR :",conditionMessage(e), "\n")}) metout <- as.data.frame(micro$metout) soil <- as.data.frame(micro$soil) soilp <- as.data.frame(micro$soilpot) met_month<- list() met_month_min<- list() for(z in 1:length(month)){ day = unique(metout$DOY)[month[z]] DAY <- day TIME <- subset(metout, DOY == day)$TIME RAD <- subset(metout, DOY == day)$SOLR AIRT <- subset(metout, DOY == day)$TALOC SOILT <- subset(soil, DOY == day)$D0cm RH <- subset(metout, DOY == day)$RHLOC V <- subset(metout, DOY == day)$VLOC SOILPOT<- subset(soilp, DOY == day)$PT0cm*10 met_month<- data.frame(DAY,TIME, RAD, AIRT, SOILT, RH, V, SOILPOT) modif <- data.frame(DAY=rep(DAY,1440),TIME=1:1440, RAD=NA, AIRT=NA, SOILT=NA, RH=NA, V=NA, SOILPOT=NA) for(j in 3:8){ f <- approxfun(met_month[,2], met_month[,j]) S <- f(seq(0, 1380, length = 1440)) modif[,j] <- S } met_month_min[[z]]<- modif } data_sun[[i]]<- do.call("rbind", met_month_min) } iter=1440 limits_data_frame_sp<- list() for(k in 1:length(data_sun)){ # Body temperature and EWL for each month limits_data_frame_sp[[k]]<- data.frame(X=NA,Y=NA,month=NA,h_t_1=NA,h_ht_1=NA,m_h_ht_1=NA,sd.h_ht_1=NA,h_t_2=NA,h_ht_2=NA,m_h_ht_2=NA,sd.h_ht_2=NA,h_t_3=NA,h_ht_3=NA,m_h_ht_3=NA,sd.h_ht_3=NA) for(q in 1:length(unique(data_sun[[k]]$DAY))){ limits_data_frame_sp[[k]][q,1]<- sub_locs_df$X1[k] limits_data_frame_sp[[k]][q,2]<- sub_locs_df$X2[k] limits_data_frame_sp[[k]][q,3]<- unique(data_sun[[k]]$DAY)[q] data<- subset(data_sun[[k]],DAY==unique(data_sun[[k]]$DAY)[q]) A1 = Morph_si$A_si[1] M1 = Morph_si$M_si[1] A2 = Morph_si$A_si[2] M2 = Morph_si$M_si[2] A3 = Morph_si$A_si[3] M3 = Morph_si$M_si[3] Tb1<- rep(data$AIRT[1],1440) EWL1<- rep(0,1440) Tb2<- rep(data$AIRT[1],1440) EWL2<- rep(0,1440) Tb3<- rep(data$AIRT[1],1440) EWL3<- rep(0,1440) for(j in 2:length(Tb1)){ sun<- data$RAD[j] Ta<- data$AIRT[j] Tg<- data$SOILT[j] RH<- data$RH[j] V<- data$V[j] # SWP = data$SOILPOT[j] m_6_tb_m1<- theatmodel(Tb1[j-1], A1, M1, Ta, Tg, sun, V, RH, skin_humidity = 1, r = r_sp, posture = 1, delta = 60) Tb1[j]<- m_6_tb_m1[1] EWL1[j]<- m_6_tb_m1[2] m_6_tb_m2<- theatmodel(Tb2[j-1], A2, M2, Ta, Tg, sun, V, RH, skin_humidity = 1, r = r_sp, posture = 1, delta = 60) Tb2[j]<- m_6_tb_m2[1] EWL2[j]<- m_6_tb_m2[2] m_6_tb_m3<- theatmodel(Tb3[j-1], A3, M3, Ta, Tg, sun, V, RH, skin_humidity = 1, r = r_sp, posture = 1, delta = 60) Tb3[j]<- m_6_tb_m3[1] EWL3[j]<- m_6_tb_m3[2] } REWL1<- (EWL1/M1)*60 REWL2<- (EWL2/M2)*60 REWL3<- (EWL3/M3)*60 ###Small freg df1<- data.frame(temp=Tb1,REWL1,data$TIME) a1<- subset(df1,Tb130){ temp[y]=30 } } dat<- data.frame(W=w,W.2=w^2,W.3=w^3,temp=temp) pred_h<- predict(cubic_s,newdata=dat,level=0)/max(theoret_s$pred.h) pred_h_max<- which.max(pred_h) pred_h_sub<- pred_ht[pred_h_max:length(pred_h)] pred_ht<- predict(cubic_s_t2,newdata=dat,level=0)/max(theoret_s$pred.ht) pred_ht_max<- which.max(pred_ht) pred_ht_sub<- pred_ht[pred_ht_max:length(pred_ht)] length_h1[i]<-length(which(pred_h_sub>0.8))+pred_h_max-1 length_ht1[i]<-length(which(pred_ht_sub>0.8))+pred_ht_max-1 } ##medium freg df2<- data.frame(temp=Tb2,REWL2,data$TIME) a2<- subset(df2,Tb230){ temp[y]=30 } } dat<- data.frame(W=w,W.2=w^2,W.3=w^3,temp=temp) pred_h<- predict(cubic_s,newdata=dat,level=0)/max(theoret_s$pred.h) pred_h_max<- which.max(pred_h) pred_h_sub<- pred_ht[pred_h_max:length(pred_h)] pred_ht<- predict(cubic_s_t2,newdata=dat,level=0)/max(theoret_s$pred.ht) pred_ht_max<- which.max(pred_ht) pred_ht_sub<- pred_ht[pred_ht_max:length(pred_ht)] length_h2[i]<-length(which(pred_h_sub>0.8))+pred_h_max-1 length_ht2[i]<-length(which(pred_ht_sub>0.8))+pred_ht_max-1 } ##big froggo df3<- data.frame(temp=Tb3,REWL3,data$TIME) a3<- subset(df3,Tb330){ temp[y]=30 } } dat<- data.frame(W=w,W.2=w^2,W.3=w^3,temp=temp) pred_h<- predict(cubic_s,newdata=dat,level=0)/max(theoret_s$pred.h) pred_h_max<- which.max(pred_h) pred_h_sub<- pred_ht[pred_h_max:length(pred_h)] pred_ht<- predict(cubic_s_t2,newdata=dat,level=0)/max(theoret_s$pred.ht) pred_ht_max<- which.max(pred_ht) pred_ht_sub<- pred_ht[pred_ht_max:length(pred_ht)] length_h3[i]<-length(which(pred_h_sub>0.8))+pred_h_max-1 length_ht3[i]<-length(which(pred_ht_sub>0.8))+pred_ht_max-1 } limits_data_frame_sp[[k]][q,4]<- nrow(a1) limits_data_frame_sp[[k]][q,5]<- median(length_ht1) limits_data_frame_sp[[k]][q,6]<- gm_mean(length_ht1) limits_data_frame_sp[[k]][q,7]<- nrow(a2) limits_data_frame_sp[[k]][q,8]<- median(length_ht2) limits_data_frame_sp[[k]][q,9]<- gm_mean(length_ht2) limits_data_frame_sp[[k]][q,10]<- nrow(a3) limits_data_frame_sp[[k]][q,11]<- median(length_ht3) limits_data_frame_sp[[k]][q,12]<- gm_mean(length_ht3) } } sp_limits_summary<- do.call("rbind", limits_data_frame_sp) write.csv(sp_limits_summary,"Spea_hours_restriction.csv")