BruchPrg.R 15.1 KB
Newer Older
1
2
3
4
library(openxlsx)
# ============================================================

dataBruch = read.xlsx("BDalosesBruch.xlsx")
5
#dataBruch = read.xlsx("/home/patrick/Documents/AA/CC et migrateur/thèse Poulet/BDalosesBruch.xlsx")
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
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){
51
52
53
54
55
56
57
58
  #  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)
}
59
60
61
62

# ------------------------------------------------
# temperature effect on growth
# ------------------------------------------------
Lambert Patrick's avatar
Lambert Patrick committed
63
64
65
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")])
Lambert Patrick's avatar
Lambert Patrick committed
66
67
TrefAtSea = (12 + Tref) / 2

Lambert Patrick's avatar
Lambert Patrick committed
68
temperature = seq(0,30,.1)
Lambert Patrick's avatar
Lambert Patrick committed
69
70
# temperature effect on growth
tempEffects = temperatureEffect(temperature,  3, 17, 26)
Lambert Patrick's avatar
Lambert Patrick committed
71
plot(temperature,tempEffects, type = 'l')
Lambert Patrick's avatar
Lambert Patrick committed
72

Lambert Patrick's avatar
Lambert Patrick committed
73
74
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)
Lambert Patrick's avatar
Lambert Patrick committed
75

Lambert Patrick's avatar
Lambert Patrick committed
76
77
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)
78

Lambert Patrick's avatar
Lambert Patrick committed
79
plot(tempData$Year[sel], tempData$Winter[sel], type = 'l')
80
81
82
83
84

# ----------------------------------------------
# growth simulation
# ----------------------------------------------
vonBertalaffyGrowth = function(age, L0, Linf, K){
Lambert Patrick's avatar
Lambert Patrick committed
85
86
  t0 = log(1 - L0 / Linf) / K
  return(Linf * (1 - exp(-K * (age - t0))))
87
88
}

Lambert Patrick's avatar
Lambert Patrick committed
89
vonBertalaffyInverse = function(L, L0, Linf, K){
Lambert Patrick's avatar
Lambert Patrick committed
90
91
  t0 = log(1 - L0/Linf)/K
  return(t0 - log(1 - L/Linf)/K)
Lambert Patrick's avatar
Lambert Patrick committed
92
93
}

Lambert Patrick's avatar
Lambert Patrick committed
94
95
KvonBertalaffy = function(age, L, L0, Linf) {
  return(-log((Linf - L ) / (Linf - L0)) / age)
Lambert Patrick's avatar
Lambert Patrick committed
96
97
98
99
}

#vonBertalaffyInverse(L=40, L0, Linf, koptMale)

Lambert Patrick's avatar
Lambert Patrick committed
100
101
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))
Lambert Patrick's avatar
Lambert Patrick committed
102
103
}

Lambert Patrick's avatar
Lambert Patrick committed
104
105
Pauly = function(age, t0, Linf, K, D){
  return(Linf / 10 * ((1 - exp(-K * D * (age - t0)))^(1 / D)))
106
107
}

Lambert Patrick's avatar
Lambert Patrick committed
108
109
110
# 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)
Lambert Patrick's avatar
Lambert Patrick committed
111
112
  L = matrix(nrow = nStep + 1)
  L[1] = L0
113
  for (i in 1:nStep) {
Lambert Patrick's avatar
Lambert Patrick committed
114
    mu = log((Linf - L[i]) * (1 - exp(-K * deltaT))) - sigma * sigma / 2
115
    increment = exp(rnorm(1, mu, sigma))
Lambert Patrick's avatar
Lambert Patrick committed
116
117
    if (withTempEffect) {
      increment = increment * tempEffect[((i - 1) %% 4) + 1]
118
    }
Lambert Patrick's avatar
Lambert Patrick committed
119
    L[i + 1] = L[i] + increment
120
121
122
123
  }
  return(L)
}

Lambert Patrick's avatar
Lambert Patrick committed
124
125
126
127
# ---------------------------------------------------
# parametres
# ----------------------------------------------------
L0 = 2
128
Linf = 70
Lambert Patrick's avatar
Lambert Patrick committed
129
130
131
#kopt = 0.3900707
#koptMale = 0.2
#koptFemale = 0.3
132
timestep = 0.25 # time step of the simumation
Lambert Patrick's avatar
Lambert Patrick committed
133
134
sigma = 0.2

Lambert Patrick's avatar
Lambert Patrick committed
135
tempEffect = mean(temperatureEffect( TrefAtSea, 3, 17, 26))
Lambert Patrick's avatar
Lambert Patrick committed
136

Lambert Patrick's avatar
Lambert Patrick committed
137
#c(mean(temperatureEffect( Tref, 3, 17, 26)), mean(temperatureEffect( TrefAtSea, 3, 17, 26)))
Lambert Patrick's avatar
Lambert Patrick committed
138
age = seq(0,8,.25)
Lambert Patrick's avatar
Lambert Patrick committed
139

140

Lambert Patrick's avatar
Lambert Patrick committed
141
142
143
144
# ###########################################################################"
# comparaison male femelle
# ###########################################################################
Lfemale = 55
Lambert Patrick's avatar
Lambert Patrick committed
145
(koptFemale = KvonBertalaffy(5.5, L = Lfemale, L0=L0, Linf=Linf)/tempEffect)
Lambert Patrick's avatar
Lambert Patrick committed
146
147
148

Lmale = 40
(koptMale = KvonBertalaffy(4.5, L = Lmale, L0, Linf)/tempEffect)
149

Lambert Patrick's avatar
Lambert Patrick committed
150

Lambert Patrick's avatar
Lambert Patrick committed
151
plot(x=NaN, xlim=range(age), ylim=c(0, Linf), xlab ='age (year)', ylab ='length (cm)')
152
for (i in 1:100) {
Lambert Patrick's avatar
Lambert Patrick committed
153
  lines(age, vonBertalaffyIncrement(max(age)/timestep, L0, Linf, koptMale, timestep, sigma, withTempEffect = TRUE), col='blue')
154
}
Lambert Patrick's avatar
Lambert Patrick committed
155
156
lines(age, vonBertalaffyGrowth(age, L0, Linf, koptMale * tempEffect), lty=2, lwd = 2, col = 'black')

Lambert Patrick's avatar
Lambert Patrick committed
157
158
res=matrix(nrow = max(age)/timestep + 1, ncol= 100 )
# verification par rapport à la médiane et la moyenne
159
for (i in 1:100) {
Lambert Patrick's avatar
Lambert Patrick committed
160
  res[,i]=vonBertalaffyIncrement(max(age)/timestep, L0, Linf, koptMale, timestep, sigma, withTempEffect = TRUE)
161
}
Lambert Patrick's avatar
Lambert Patrick committed
162
163
lines(age, apply(res, 1, quantile, probs =0.5 ),  lty=2, lwd = 2, col = 'red')
lines( age, rowMeans(res))
164
165

for (i in 1:100) {
Lambert Patrick's avatar
Lambert Patrick committed
166
  lines(age, vonBertalaffyIncrement(max(age)/timestep, L0, Linf, koptFemale, timestep, sigma, withTempEffect = TRUE), col='pink')
167
}
Lambert Patrick's avatar
Lambert Patrick committed
168
169
170
lines(age, vonBertalaffyGrowth(age, L0, Linf, koptFemale * tempEffect), lty=2, lwd = 2)
abline(h = Lmale)
abline(h = Lfemale)
171

Lambert Patrick's avatar
Lambert Patrick committed
172
173
174
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))
175

Lambert Patrick's avatar
Lambert Patrick committed
176
ageFemale=vonBertalaffyInverse(L=Lfemale, L0, Linf, koptFemale*tempEffect)
Lambert Patrick's avatar
Lambert Patrick committed
177
178
179
180
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 =" ")
181

Lambert Patrick's avatar
Lambert Patrick committed
182
183
184
185
186
187
188
189
190
# --------------------------------------------------------------
# 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))
191

Lambert Patrick's avatar
Lambert Patrick committed
192
193
194
195
# ======================================================================
#  temp en fonction de la latitude
# ====================================================================
library(dplyr)
196

Lambert Patrick's avatar
Lambert Patrick committed
197
198
199
200
201
202
extractTemp<-as.data.frame(tempData %>%
                             filter(Year>1901 & Year <= 1910) %>%
                             group_by(NOM,Lat)%>%
                             summarize(meanSpring=mean(Spring)) 
                           %>% arrange (Lat)
)
203

Lambert Patrick's avatar
Lambert Patrick committed
204
plot(extractTemp$Lat, extractTemp$meanSpring, pch=20)
205
abline(v=45.25) # Garonne
Lambert Patrick's avatar
Lambert Patrick committed
206
abline(v=49.25) # vire
207
abline(v=52.25) # Rhine
208
209


Lambert Patrick's avatar
Lambert Patrick committed
210
211
abline(h=11)
abline(h=12.5)
212
213
214
# ======================================================================
# exploration of growth for male and female
# ======================================================================
215
correction=mean(temperatureEffect(TrefAtSea, 3, 17, 26))
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240

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
# =====================================================================
Lambert Patrick's avatar
Lambert Patrick committed
241
simData=read.csv("/home/patrick/Documents/workspace/GR3D/data/output/lengthAgeDistribution_1-RCP85.csv", sep=";", row.names = NULL)
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262

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)
263
264
265
266
267
268
269

# ===================================================================
# 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. 

270
271
#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° 
272
#Ratio = 0.2 
273
#n= paramètre simulant l'effet Allee 
274
275


276
#-----------On cherche a reproduire la relation SR telle que modélisée dans GR3D-------------- 
277
278
279
280
281
282
283
284
285
286
287
288

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)
  
}

289
#Relation SR telle qu'elle est modélisée dans GR3D
290
291
292
293
294
295
296
297
298
299
300
301

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
  
302
  #parametre c de la RS de BH intégrant un effet du BV considéré 
303
304
  cj = lambda/surfaceWatershed
  
305
  #parametre b représentant la mortalité densité dépendante de la RS de BH intégrant un effet de la temperature
306
307
308
309
310
  # 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
  
311
  #parametre a (fécondité de l'espèce) de la RS de BH intégrant un effet de la temperature
312
313
  alphaj = (bj * exp(-bj * deltaTrecruitement)) / (cj * (1-exp(-bj * deltaTrecruitement)))
  
314
  #Bj = paramètre de la relation SR intégrant l'effet de la température 
315
316
  betaj = bj/(a*cj*(1-exp(-bj*deltaTrecruitement)))
  
317
  #p = proportion de géniteurs participant à la reproduction en focntion de la quantité de géniteur total
318
319
320
321
322
323
324
  #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)))
  
325
  #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 
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
  
  #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")
  
345
#-----------On cherche à déterminer le numbre de juvéniles générés par S = 100000 géniteurs en fonction de la T° -------------- 
346
347
348
349
350
351
352
353
354
355
356
357
358

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
    
359
    #parametre c de la RS de BH intégrant un effet du BV considéré 
360
361
    cj = lambda/surfaceWatershed
    
362
    #parametre b représentant la mortalité densité dépendante de la RS de BH intégrant un effet de la temperature
363
364
365
366
367
    # 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
    
368
    #parametre a (fécondité de l'espèce) de la RS de BH intégrant un effet de la temperature
369
370
    alphaj = (bj * exp(-bj * deltaTrecruitement)) / (cj * (1-exp(-bj * deltaTrecruitement)))
    
371
    #Bj = paramètre de la relation SR intégrant l'effet de la température 
372
373
    betaj = bj/(a*cj*(1-exp(-bj*deltaTrecruitement)))
    
374
    #p = proportion de géniteurs participant à la reproduction en focntion de la quantité de géniteur total
375
376
377
378
379
380
381
    #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)))
    
382
    #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 
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
    
    #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", 
400
       xlab =" Temperature (°C)",
401
402
403
       ylab = "Number Of Recruits")
  

404
#-----------On cherche à déterminer le numbre de géniteurs survivants en fonction de la T° -------------- 
405
406
407
408
409

  #Prend en compte Zsea 
  
  plot(temperature, StockRecruitementRelationship (temperature, 84810, numberOfSpawner),
       type="l", 
410
       xlab =" Temperature (°C)",
411
412
413
414
415
416
417
       ylab = "Number Of Recruits")