library(openxlsx) # ============================================================ dataBruch = read.xlsx("BDalosesBruch.xlsx") #dataBruch = read.xlsx("/home/patrick/Documents/AA/CC et migrateur/thèse Poulet/BDalosesBruch.xlsx") dataBruch$`M.gonades.(g)` = as.numeric(dataBruch$`M.gonades.(g)`) head(dataBruch) tapply(dataBruch$`Lf.(cm)`, dataBruch[,c('Année', 'Sexe')],min, na.rm = TRUE) tapply(dataBruch$`Lf.(cm)`, dataBruch[,c('Année', 'Sexe')],min, na.rm = TRUE) tapply(dataBruch$`Lf.(cm)`, dataBruch[,c('Année', 'Sexe')],quantile, na.rm = TRUE, probs=.05) tapply(dataBruch$`Lf.(cm)`, dataBruch[,c('Année', 'Sexe')],max, na.rm = TRUE) tapply(dataBruch$`Lf.(cm)`, dataBruch[,c('Année', 'Sexe')],quantile, na.rm = TRUE, probs=.95) tapply(dataBruch$`Lf.(cm)`, dataBruch[,c('Année', 'Sexe')],quantile, na.rm = TRUE, probs=.5) sel = dataBruch$Année==2013 & dataBruch$Sexe =='F' hist(dataBruch$`Lt.(cm)`[sel]) abline(v=quantile(dataBruch$`Lt.(cm)`[sel], probs = 0.05)) sel = dataBruch$Sexe=='M' lm (dataBruch$`Lt.(cm)`[sel]~dataBruch$`Lf.(cm)`[sel]) summary(lm (dataBruch$`Lt.(cm)`~dataBruch$`Lf.(cm)`)) summary(lm (dataBruch$`Lt.(cm)`~dataBruch$`Lf.(cm)` * dataBruch$Sexe)) # ==================================================== # fecundity # Taverny 1991 # ==================================================== (41*172895+33*202902+74*186424)/(41+33+74) (41*98390+33*110386+74*104325)/(41+33+74) # ============================================================ # maximal production of recruit in GR3D for the Garonne basin # ============================================================ bj=-log(1.7e-3) /.33 cj= 4.1e-4 / (84810*.5356) (alphaj = bj*exp(-bj*.33)/(cj*(1-exp(-bj*.33)))) # ================================================ # growth in GR3D # ================================================ temperatureEffect= function(temp, Tmin, Topt, Tmax){ # if (temp<=Tmin | temp >= Tmax) # return(0) # else response=(temp-Tmin)*(temp-Tmax)/((temp-Tmin)*(temp-Tmax)-(temp-Topt)^2) response[temp<=Tmin | temp >= Tmax] = 0 return(response) } # ------------------------------------------------ # temperature effect on growth # ------------------------------------------------ tempData = read.csv("/home/patrick/Documents/workspace/GR3D/data/input/reality/SeasonTempBVFacAtlant1801_2100_newCRU_RCP85.csv", sep =";") sel = tempData$NOM == "Garonne" & tempData$Yea >= 2008 & tempData$Year <= 2018 Tref = colMeans(tempData[sel, c("Winter", "Spring", "summer", "Autumn")]) TrefAtSea = (12 + Tref) / 2 temperature = seq(0,30,.1) # temperature effect on growth tempEffects = temperatureEffect(temperature, 3, 17, 26) plot(temperature,tempEffects, type = 'l') points(TrefAtSea, temperatureEffect(TrefAtSea, 3, 17, 26), col = "red", pch = 20) text(TrefAtSea, temperatureEffect(TrefAtSea, 3, 17, 26), c("Winter", "Spring", "Summer", "Autumn"), pos = 4) points(Tref, temperatureEffect(Tref, 3, 17, 26), col = "red") text(Tref, temperatureEffect(Tref, 3, 17, 26), c("Winter", "Spring", "Summer", "Autumn"), pos = 2, cex = .5) plot(tempData$Year[sel], tempData$Winter[sel], type = 'l') # ---------------------------------------------- # growth simulation # ---------------------------------------------- vonBertalaffyGrowth = function(age, L0, Linf, K){ t0 = log(1 - L0 / Linf) / K return(Linf * (1 - exp(-K * (age - t0)))) } vonBertalaffyInverse = function(L, L0, Linf, K){ t0 = log(1 - L0/Linf)/K return(t0 - log(1 - L/Linf)/K) } KvonBertalaffy = function(age, L, L0, Linf) { return(-log((Linf - L ) / (Linf - L0)) / age) } #vonBertalaffyInverse(L=40, L0, Linf, koptMale) Zfct = function(age, Ltarget, L0, Linf, K, Tref ) { return(sapply(age,function(a) (vonBertalaffyGrowth(a, L0, Linf, K) * mean(temperatureEffect(Tref, 3, 17, 26)) - Ltarget)^2)) } Pauly = function(age, t0, Linf, K, D){ return(Linf / 10 * ((1 - exp(-K * D * (age - t0)))^(1 / D))) } # pas cohérent avec la temperature effet sur le coeff de croissance mais ca marche vonBertalaffyIncrement = function(nStep, L0, Linf, K, deltaT, sigma, withTempEffect=FALSE, TrefAtSea = c(9.876946, 13.489854, 15.891487, 11.554104) ){ tempEffect = temperatureEffect( TrefAtSea , 3, 17, 26) L = matrix(nrow = nStep + 1) L[1] = L0 for (i in 1:nStep) { mu = log((Linf - L[i]) * (1 - exp(-K * deltaT))) - sigma * sigma / 2 increment = exp(rnorm(1, mu, sigma)) if (withTempEffect) { increment = increment * tempEffect[((i - 1) %% 4) + 1] } L[i + 1] = L[i] + increment } return(L) } # --------------------------------------------------- # parametres # ---------------------------------------------------- L0 = 2 Linf = 70 #kopt = 0.3900707 #koptMale = 0.2 #koptFemale = 0.3 timestep = 0.25 # time step of the simumation sigma = 0.2 tempEffect = mean(temperatureEffect( TrefAtSea, 3, 17, 26)) #c(mean(temperatureEffect( Tref, 3, 17, 26)), mean(temperatureEffect( TrefAtSea, 3, 17, 26))) age = seq(0,8,.25) # ###########################################################################" # comparaison male femelle # ########################################################################### Lfemale = 55 (koptFemale = KvonBertalaffy(5.5, L = Lfemale, L0=L0, Linf=Linf)/tempEffect) Lmale = 40 (koptMale = KvonBertalaffy(4.5, L = Lmale, L0, Linf)/tempEffect) plot(x=NaN, xlim=range(age), ylim=c(0, Linf), xlab ='age (year)', ylab ='length (cm)') for (i in 1:100) { lines(age, vonBertalaffyIncrement(max(age)/timestep, L0, Linf, koptMale, timestep, sigma, withTempEffect = TRUE), col='blue') } lines(age, vonBertalaffyGrowth(age, L0, Linf, koptMale * tempEffect), lty=2, lwd = 2, col = 'black') res=matrix(nrow = max(age)/timestep + 1, ncol= 100 ) # verification par rapport à la médiane et la moyenne for (i in 1:100) { res[,i]=vonBertalaffyIncrement(max(age)/timestep, L0, Linf, koptMale, timestep, sigma, withTempEffect = TRUE) } lines(age, apply(res, 1, quantile, probs =0.5 ), lty=2, lwd = 2, col = 'red') lines( age, rowMeans(res)) for (i in 1:100) { lines(age, vonBertalaffyIncrement(max(age)/timestep, L0, Linf, koptFemale, timestep, sigma, withTempEffect = TRUE), col='pink') } lines(age, vonBertalaffyGrowth(age, L0, Linf, koptFemale * tempEffect), lty=2, lwd = 2) abline(h = Lmale) abline(h = Lfemale) ageMale=vonBertalaffyInverse(L=Lmale, L0, Linf, koptMale*tempEffect) points(ageMale, vonBertalaffyGrowth(ageMale, L0, Linf, koptMale * tempEffect), col = 'red') segments(x0 = ageMale, y0=0, y1=vonBertalaffyGrowth(ageMale, L0, Linf, koptMale * tempEffect)) ageFemale=vonBertalaffyInverse(L=Lfemale, L0, Linf, koptFemale*tempEffect) points(ageFemale, vonBertalaffyGrowth(ageFemale, L0, Linf, koptFemale * tempEffect), col = 'red') segments(x0 = ageFemale, y0=0, y1=vonBertalaffyGrowth(ageFemale, L0, Linf, koptFemale * tempEffect)) cat(koptFemale, koptMale, sep =" ") cat(ageFemale, ageMale, sep =" ") # -------------------------------------------------------------- # croissance sur le Rhin # -------------------------------------------------------------- sel = tempData$NOM=="Rhine" & tempData$Year>=2008 & tempData$Year<=2018 TrefRhine = (12 + colMeans(tempData[sel, c("Winter", "Spring", "summer", "Autumn")])) / 2 rbind(TrefAtSea, TrefRhine) tempEffectRhine = mean(temperatureEffect( TrefRhine, 3, 17, 26)) c(vonBertalaffyInverse(L=Lfemale, L0, Linf, koptFemale*tempEffectRhine), vonBertalaffyInverse(L=Lmale, L0, Linf, koptMale*tempEffectRhine)) # ====================================================================== # temp en fonction de la latitude # ==================================================================== library(dplyr) extractTemp<-as.data.frame(tempData %>% filter(Year>1901 & Year <= 1910) %>% group_by(NOM,Lat)%>% summarize(meanSpring=mean(Spring)) %>% arrange (Lat) ) plot(extractTemp$Lat, extractTemp$meanSpring, pch=20) abline(v=45.25) # Garonne abline(v=49.25) # vire abline(v=52.25) # Rhine abline(h=11) abline(h=12.5) # ====================================================================== # exploration of growth for male and female # ====================================================================== correction=mean(temperatureEffect(TrefAtSea, 3, 17, 26)) age=seq(0,10,.25) present = vonBertalaffyGrowth(age, 2, 60, 0.3900707 * correction) plot(age, present, type='l', lwd=3, ylim =c(0,80)) present[age == 5] abline(v=5) male = vonBertalaffyGrowth(age, 2, 65, 0.3900707 * correction) lines(age, male, type='l', lwd=3, ylim =c(0,80), col='blue') abline(h=40, col='blue', lwd=2, lty=2) male[age == 5] female = vonBertalaffyGrowth(age, 2, 75, 0.3900707*55/40 * correction) lines(age, female, lwd=3, col='red') abline(h=55, col='red', lwd=2, lty=2) female[age == 5] ## a partir d'individus en mer donc à croissance de plus lente à mesure qu'ils sont agées (taverny = Pauly(age,t0=-0.7294, Linf=701.59, K=0.4491, D=.5912)) lines (age, taverny, lwd=2, col ='green') taverny[age == 5] # =================================================================== # GR3D outputs # ===================================================================== simData=read.csv("/home/patrick/Documents/workspace/GR3D/data/output/lengthAgeDistribution_1-RCP85.csv", sep=";", row.names = NULL) simGaronne= simData[simData$basin =="Garonne",] sel=simGaronne$nbSpawn == 0 tapply(simGaronne$length[sel], simGaronne[sel,c('year')],quantile, na.rm = TRUE, probs=.5) # masse des gonades avant sel = (dataBruch$LOT =='Tuilières' | dataBruch$LOT =='Golfech') & !is.na(dataBruch$`M.gonades.(g)`) & dataBruch$Sexe =='F' mean(dataBruch$`M.gonades.(g)`[sel]/dataBruch$`M.tot.(g)`[sel]) sel = (dataBruch$LOT =='Tuilières' | dataBruch$LOT =='Golfech') & dataBruch$Sexe =='F' sum(sel) Wpre = mean(dataBruch$`M.tot.(g)`[sel]) Wgonad =mean(dataBruch$`M.gonades.(g)`[sel], na.rm = TRUE) sel = (! (dataBruch$LOT =='Tuilières' | dataBruch$LOT =='Golfech')) & dataBruch$Sexe =='F' Wpost= mean(dataBruch$`M.tot.(g)`[sel]) WgonadSpent =mean(dataBruch$`M.gonades.(g)`[sel], na.rm = TRUE) (Wloss=(Wpre - Wpost)/Wpre) # =================================================================== # Exploration of Stock recruitement-relationship for GR3D calibration # ===================================================================== #Use to improve the likelihood between observations and GR3D outputs in terms of abudances and North limit colonization. #a = fécondité de l'espèce, a = 135000 #S = quantité de géniteurs: ici on veut la quantité R0 produite par 1000 géniteurs en fonction de la T° #Ratio = 0.2 #n= paramètre simulant l'effet Allee #-----------On cherche a reproduire la relation SR telle que modélisée dans GR3D-------------- temperatureEffect= function(tempWater, TminRep, ToptRep, TmaxRep){ # if (tempWater<=TminRep | tempWater >= TmaxRep) # return(0) # else response=(tempWater-TminRep)*(tempWater-TmaxRep)/((tempWater-TminRep)*(tempWater-TmaxRep)-(tempWater-ToptRep)^2) response[tempWater<=TminRep | tempWater >= TmaxRep] = 0 return(response) } #Relation SR telle qu'elle est modélisée dans GR3D numberOfSpawner<- seq(0:400000) StockRecruitementRelationship <-function (temp, surfaceWatershed, S) { lambda = 4.1E-4 deltaTrecruitement = 0.33 survOptRep = 0.0017 n= 2.4 ratioTeta = 1.9 a = 135000 #parametre c de la RS de BH intégrant un effet du BV considéré cj = lambda/surfaceWatershed #parametre b représentant la mortalité densité dépendante de la RS de BH intégrant un effet de la temperature # bj = (-(1/deltaTrecruitement))* # log(survOptRep * temperatureEffect(temp, 9.8, 20.0, 26.0)) bj = - log(survOptRep * temperatureEffect(temp, 9.8, 20.0, 26.0)) / deltaTrecruitement #parametre a (fécondité de l'espèce) de la RS de BH intégrant un effet de la temperature alphaj = (bj * exp(-bj * deltaTrecruitement)) / (cj * (1-exp(-bj * deltaTrecruitement))) #Bj = paramètre de la relation SR intégrant l'effet de la température betaj = bj/(a*cj*(1-exp(-bj*deltaTrecruitement))) #p = proportion de géniteurs participant à la reproduction en focntion de la quantité de géniteur total #p = 1/(1+exp(-log(19)*(S-n)/(Ratio*surfaceWatershed))) S95 = n * surfaceWatershed S50 = S95/ratioTeta p= 1/(1+exp(-log(19)*(S-S50)/(S95-S50))) #relation Stock Recrutement ie calcul le nombre de recrues en fonction du nombre de géniteurs et de la T en intégrant l'effet Allee #R0 = aj * S * p AlleeEffect = 1/ (1+exp(-log(19)*(S -n/ratioTeta*surfaceWatershed)/(n*surfaceWatershed -n/ratioTeta*surfaceWatershed))) Rj = (alphaj * S * AlleeEffect)/(betaj +S * AlleeEffect) #Rj = ((aj * S) * p)/(Bj +S * p) StockRecruitement = as.data.frame(Rj) return (Rj) } StockRecruitement<-StockRecruitementRelationship (18, 84810, numberOfSpawner) plot(numberOfSpawner, StockRecruitement,type = 'l', xlab= "Number of spawners", ylab = "Number of recruits") #-----------On cherche à déterminer le numbre de juvéniles générés par S = 100000 géniteurs en fonction de la T° -------------- temperature <- seq (8,30,.1) numberOfSpawner=100000 StockRecruitementRelationship <-function (temp, surfaceWatershed, S) { lambda = 4.1E-4 deltaTrecruitement = 0.33 survOptRep = 0.0017 n= 2.4 ratioTeta = 1.9 a = 135000 #parametre c de la RS de BH intégrant un effet du BV considéré cj = lambda/surfaceWatershed #parametre b représentant la mortalité densité dépendante de la RS de BH intégrant un effet de la temperature # bj = (-(1/deltaTrecruitement))* # log(survOptRep * temperatureEffect(temp, 9.8, 20.0, 26.0)) bj = - log(survOptRep * temperatureEffect(temp, 9.8, 20.0, 26.0)) / deltaTrecruitement #parametre a (fécondité de l'espèce) de la RS de BH intégrant un effet de la temperature alphaj = (bj * exp(-bj * deltaTrecruitement)) / (cj * (1-exp(-bj * deltaTrecruitement))) #Bj = paramètre de la relation SR intégrant l'effet de la température betaj = bj/(a*cj*(1-exp(-bj*deltaTrecruitement))) #p = proportion de géniteurs participant à la reproduction en focntion de la quantité de géniteur total #p = 1/(1+exp(-log(19)*(S-n)/(Ratio*surfaceWatershed))) S95 = n * surfaceWatershed S50 = S95/ratioTeta p= 1/(1+exp(-log(19)*(S-S50)/(S95-S50))) #relation Stock Recrutement ie calcul le nombre de recrues en fonction du nombre de géniteurs et de la T en intégrant l'effet Allee #R0 = aj * S * p AlleeEffect = 1/ (1+exp(-log(19)*(S -n/ratioTeta*surfaceWatershed)/(n*surfaceWatershed -n/ratioTeta*surfaceWatershed))) Rj = (alphaj * S * AlleeEffect)/(betaj +S * AlleeEffect) #Rj = ((aj * S) * p)/(Bj +S * p) StockRecruitement = as.data.frame(Rj) return (Rj) } plot(temperature, StockRecruitementRelationship (temperature, 84810, numberOfSpawner), type="l", xlab =" Temperature (°C)", ylab = "Number Of Recruits") #-----------On cherche à déterminer le numbre de géniteurs survivants en fonction de la T° -------------- #Prend en compte Zsea plot(temperature, StockRecruitementRelationship (temperature, 84810, numberOfSpawner), type="l", xlab =" Temperature (°C)", ylab = "Number Of Recruits")