-
Notifications
You must be signed in to change notification settings - Fork 1
/
GOM.R
198 lines (170 loc) · 11.2 KB
/
GOM.R
1
2
3
4
5
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
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
#EMAX GOM model as modified by Gaichas for the herring assessment
#and by Lucey for QNM comparison
#modified 8 July 2019 from https://github.com/raw/NOAA-EDAB/QNM/master/QNMcompare/GOM.R
library(Rpath); library(data.table)
gom.groups <- c("Phytoplankton- Primary Producers", "Bacteria", "Microzooplankton",
"Small copepods", "Large Copepods", "Gelatinous Zooplankton", "Micronekton",
"Macrobenthos- polychaetes", "Macrobenthos- crustaceans", "Macrobenthos- molluscs",
"Macrobenthos- other", "Megabenthos- filterers", "Megabenthos- other",
"Shrimp et al.", "Larval-juv fish- all", "Small Pelagics- commercial",
"Small Pelagics- other", "Small Pelagics- squid", "Small Pelagics- anadromous",
"Medium Pelagics- (piscivores & other)", "Demersals- benthivores",
"Demersals- omnivores", "Demersals- piscivores", "Sharks- pelagics",
"HMS", "Pinnipeds", "Baleen Whales", "Odontocetes", "Sea Birds",
"Discard", "Detritus-POC", "Fishery")
gom.type <- c(1, rep(0, 28), rep(2, 2), 3)
gom.par <- create.rpath.params(gom.groups, gom.type)
gom.par$model[, Biomass := c(22.126, 5.484, 4.885, 10.403, 11.955, 1.283, 4.874,
18.942, 4.04, 9.866, 24.936, 2.879, 3.505, 0.396, 0.207,
5.714, 1.275, 0.29, 0.153, 0.0229, 2.981, 0.4, 4.006,
0.00296, 0.00587, 0.063, 0.602, 0.0336, 0.0035, NA,
NA, NA)]
gom.par$model[, PB := c(163.143, 91.25, 72, 30.918, 35, 35, 14.25, 2.55, 3.3, 2.24,
2.04, 0.864, 1.68, 2, 15, 0.52, 0.44, 1.4, 0.437, 0.649,
0.459, 0.54, 0.55, 0.15, 0.5, 0.0673, 0.042, 0.04, 0.275,
rep(NA, 3))]#55.84456, 29.90675, NA)]
gom.par$model[, QB := c(0, 380.208, 242.424, 127.75, 109.5, 146, 36.5, 17.5, 21,
13.72, 11.777, 10, 11.03, 5, 45, 1.882, 2, 2, 2, 1.428,
0.9, 0.9, 1.014, 0.623, 2.362, 4.85, 2.3, 8.5, 5.362,
rep(NA, 3))]
#float path_EE[NUM_GROUPS+1] = { 1, 0.8798624, 0.8832949, 0.9413487, 0.9117724, 0.7440153, 0.9114986, 0.8848614, 0.8999999, 0.8819566, 0.8633214, 0.8862541, 0.8823805, 0.8879386, 0.8046172, 0.8990406, 0.8963319, 0.9790097, 0.9941966, 0.5794412, 0.8889518, 0.9287031, 0.910743, 0.878319, 0.8430014, 0.8114244, 0.2865187, 0.01683634, 0.1337793, 0.1249854, 0.005459445, 0.9885055, 0};
gom.par$model[, BioAcc := c(rep(0, 31), NA)]
gom.par$model[, Unassim := c(0, 0.2, 0.1, 0.25, 0.25, 0.35, 0.25, 0.5, 0.5, 0.6,
0.5, 0.7, 0.3, 0.3, 0.15, 0.15, 0.35, 0.15, 0.15,
0.15, 0.3, 0.35, 0.15, 0.15, 0.15, 0.2, 0.2, 0.2,
0.15, 0, 0, NA)]
#float path_DtImp[NUM_GROUPS+1] = { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
gom.par$model[, Discard := c(rep(0, 31), 1)]
gom.par$model[, 'Detritus-POC' := c(rep(1, 29), rep(0, 3))]
gom.par$model[, Fishery := c(rep(0, 11), 0.0917, 0.345, 0.0622, 0, 0.874, 0.0000151,
0.0138, 0.00561, 0.0101, 0.142, 0.00718, 0.301, 0.00024,
0.00182, rep(0, 6), NA)]
gom.par$model[, Fishery.disc := c(rep(0, 5), 6.36E-07, 0, 0.000135, 0.0000183, 0.000189,
0.000724, 0.0302, 0.104, 0.0195, 2.65E-09, 0.131,
0.000156, 0.000562, 0.0167, 0.00303, 0.043, 0.00228,
0.0904, 0.0000721, 0.000545, 0.00115, 0.000405,
0.000139, 0.0000588, rep(0, 2), NA)]
DC.groups <- gom.groups[which(!gom.groups %in% c('Discard', 'Detritus-POC', 'Fishery'))]
DC.c.code <-as.data.table(matrix(c(
0,0,0.4305193,0.2117425,0.7240587,0.6317893,0.08602893,0.1824,0.1241,0.1351713,0.4000163,0.2008283,0.7909901,0,0.07548419,0.0589119,0.01121766,0.1656693,0,0.0117771,0,0,0,0,0,0,0,0,0,0,0,0,0,
0,0,0,0.1971663,0,0,0.02000673,0,0.2954761,0.1374624,0.1951299,0.2008283,0.1188099,0.1418395,0.4473137,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
0,0,0,0.09702543,0.1133734,0.04287973,0.05001682,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
0,0,0,0,0.032,0.1097328,0.3301111,0.2208001,0,0.005453273,0,0,0,0,0,0.4305102,0.1193151,0.08079781,0,0.05613042,0,0,0,0,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,0.1276366,0.3601212,0.3664,0,0.02749248,0,0,0,0,0,0.2492427,0.4344292,0.676395,0.1006037,0.9020878,0,0,0,0,0.03340645,0,0,0.153928,0,0.03427451,0,0,0,
0,0,0,0,0,0.02392604,0.04867304,0,0,0,0,0,0,0,0,0,0.01923198,0.02177701,0,0,0.01960942,0.008961361,0.08314439,0.04308448,0,0.09437622,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0.0464,0.05909523,0.05383944,0.03414774,0.06733654,0,0,0.1509684,0.06797528,0.2172146,0.04431918,0.5030184,0.02034226,0,0.06147739,0.0328972,0.02203885,0,0,0.00313569,0.5024855,0.000274618,0.1405256,0,0,0,
0,0,0,0,0,0,0,0,0.02672455,0.1351713,0,0.06561421,0,0.07919377,0,0.02096442,0.01121766,0,0,0.001553198,0,0.1471426,0.1663658,0.01101942,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,0.000327961,0,0,0.004904852,0.03310977,0,0.009938519,0,0.066255,0.008152554,0.01274636,0.0545647,0.000721625,0.1509055,0.002170789,0.01226363,0.1844322,0.1663658,0.01101942,0.01113548,0,0.005412052,0,0.002591398,0,0,0,0,
0,0,0,0,0,0,0,0,0.000639306,0.02196813,0.004286418,0.02239663,0,0.2375813,0,0.005,0.00803651,0.00072132,0,0,0,0.1229548,0.1663658,0.03305828,0,0,0.005412052,0,0.002591398,0,0,0,0,
0,0,0,0,0,0.001748159,0,0,0.02106318,0.1346288,0.022337,0.0363167,0,0.2375813,0.1048391,0.02096442,0.01121766,0.002212085,0.01676727,0,0.01226363,0.1844322,0.1109106,0.1091925,0.01113548,0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,0.001053159,0.004799868,0.002984548,0.001145771,0,0.009741648,0,0,0,0,0,0,0,0.05240696,0.04041656,0.01101942,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,0.001404212,0.004421631,0.000735662,0.006267365,0,0.04154278,0,0,0,0,0,0,0.1941741,0.09574349,0.1109106,0.02203885,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.001258382,0,0,0,0.02361985,0,0.01720621,0.01696089,0.003024275,0.1207828,0,0,0.00313569,0,0.000274618,0.02157306,0,0,0,
0,0,0,0,0,0,0.004059562,0,0,0,0,0,0,0,0,0.0747728,0.1080974,0.00738653,0.1676727,0.005938409,0.01226363,0,0.01127904,0.01101942,0,0,0,0,0,0.10,0,0,0,
0,0,0,0,0,0,0.000358531,0,0,0,0,0,0,0,0,0,0,0,0.0075355,0,0.3638209,0.04829028,0.04365378,0.2704768,0.2338451,0.1801901,0.1359761,0.1550972,0.2288119,0.3118983,0,0,0,
0,0,0,0,0,0,0.000454902,0,0,0,0,0,0,0,0,0,0,0,0.007557068,0,0.03679087,0.002015652,0.001879841,0.0560989,0.05567741,0.7006319,0.256553,0.08101361,0.06583849,0.2684838,0,0,0,
0,0,0,0,0,0,9.18E-05,0,0,0,0,0,0,0,0,0,0,0,0.02131235,0,0.02591145,0.003626702,0.001691515,0.04413503,0.1781677,0.02480174,0.06566829,0.03142953,0.3652223,0.03182262,0,0,0,
0,0,0,0,0,0,4.37E-05,0,0,0,0,0,0,0,0,0,0,0,0.001007609,0,0.01925174,0,0,0.001604244,0.02227096,0,0,0,0,0.02017987,0,0,0,
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.000721882,0,0,0,0.031213,0,0,0,0,0.000353094,0,0,0,
0,0,0,0,0,0,0,0,0,0.001607866,0,0.00181822,0,0.002898938,0,0,0.001819155,0,0,0,0.1216143,0.03728956,0.01221896,0.01101942,0.05567741,0,0.251599,0.01532696,0.1114401,0.02,0,0,0,
0,0,0,0,0,0,0,0,0,0.000314047,0,0.000115925,0,0.000863633,0,0,0.001819155,0,0,0,0.04249413,0.006046957,0.0093992,0.004007064,0.08908386,0,0.02449874,0.0110239,0.04565654,0.01322618,0,0,0,
0,0,0,0,0,0,0,0,0,0,0,0.001225253,0,0.001656536,0,0,0.001819155,0,0,0,0.1216143,0.004031305,0.005639521,0.218385,0.07794838,0,0.2486094,0.04969525,0.1772986,0.02856211,0,0,0,
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.03340645,0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.01113548,0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.03340645,0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.01113548,0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.02227096,0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.03340645,0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.0091007,0,0,0,
0,0,0.5694807,0.4940658,0.130568,0.06195949,0.1000336,0.184,0.465539226,0.304559621,0.340362397,0.386167839,0.0902,0.18084539,0.21198369,0.0589119,0,0,0,0,0,0.02418782,0.03383712,0,0.05567741,0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0), 32, 33, byrow = T))
DC.c.code[, V1 := NULL]
DC.c.code[, V31 := NULL]
DC.c.code[, V32 := NULL]
DC.c.code[, V33 := NULL]
gom.par$diet <- cbind(gom.par$diet[, Group], DC.c.code)
setnames(gom.par$diet, c(paste('V', 1:30, sep = '')), c('Group', DC.groups))
check.rpath.params(gom.par)
#Run model
GOM <- rpath(gom.par, 'Gulf of Maine')
webplot(GOM, highlight = "Sharks- pelagics", labels = T)
#Check sim
gom.base <- rsim.scenario(GOM, gom.par, 1:100)
gom.run <- rsim.run(gom.base)
rsim.plot(gom.run, gom.groups)
#scenario plots
GOM.b2 <- adjust.fishing(gom.base, parameter = 'EFFORT', group = 'Fishery',
value =2, sim.year = 25:100)
GOM.run1 <- rsim.run(GOM.b2)
rsim.plot(GOM.run1, gom.groups[1:30])
# #Create community matrices for Qpress
# GOM.10 <- Rpath2Qpress(GOM.params, .1, .1)
# GOM.20 <- Rpath2Qpress(GOM.params, .2, .2)
# GOM.30 <- Rpath2Qpress(GOM.params, .3, .3)
# GOM.40 <- Rpath2Qpress(GOM.params, .4, .4)
# GOM.50 <- Rpath2Qpress(GOM.params, .5, .5)
#
# write.csv(GOM.10, file = file.path(out.dir, 'GOM_Qpress_10.csv'), row.names = F)
# write.csv(GOM.30, file = file.path(out.dir, 'GOM_Qpress_30.csv'), row.names = F)
#
# #MTI
# #Mixed trophic impact
# MTI <- function(Rpath, Rpath.params){
# x <- copy(Rpath.params)
# f.dc <- data.table(Fishing = rep(0, nrow(x$diet)))
# DC <- as.data.table(cbind(f.dc, x$diet))
# DC <- DC[Group != 'Import']
# DC[, Group := NULL]
#
# #Need to make write functions able to output to RData
# y <- copy(Rpath)
# ngroup <- y$NUM_LIVING + y$NUM_DEAD
#
# #Calculate F mortality
# totcatch <- y$Catch + y$Discards
# Fmort <- as.data.frame(totcatch / y$BB[row(as.matrix(totcatch))])
# setnames(Fmort, paste0('V', seq_along(y$NUM_GEARS)),
# paste0('F.', y$Group[(ngroup +1):y$NUM_GROUPS]))
# Fmort <- Fmort[1:ngroup, ]
#
# #Calculate M2
# bio <- y$BB[1:y$NUM_LIVING]
# BQB <- bio * y$QB[1:y$NUM_LIVING]
# diet <- as.data.frame(y$DC)
# nodetrdiet <- diet[1:y$NUM_LIVING, ]
# detrdiet <- diet[(y$NUM_LIVING +1):ngroup, ]
# newcons <- nodetrdiet * BQB[col(as.matrix(nodetrdiet))]
# predM <- newcons / bio[row(as.matrix(newcons))]
# detcons <- detrdiet * BQB[col(as.matrix(detrdiet))]
# predM <- rbind(predM, detcons)
# setnames(predM, paste('V', 1:y$NUM_LIVING, sep = ''),
# paste('M2.', y$Group[1:y$NUM_LIVING], sep = ''))
#
# FC <- as.data.table(cbind(Fmort, predM))
#
# FC[, mort.sum := rowSums(.SD), .SDcols = names(FC)]
# FC <- FC[, .SD / mort.sum, .SDcols = names(FC)]
# FC[, mort.sum := NULL]
#
# #Merge pred and prey
# MTI <- as.matrix(DC + FC)
#
# #Add small increase
# MTI.diag <- diag(MTI) + 1
#
# diag(MTI) <- MTI.diag
#
# #Inverse
# out <- MASS::ginv(MTI)
#
# return(out)
# }
#
# MTI.GOM <- MTI(GOM, GOM.params)
#
# n.prey <- nrow(MTI.GOM)
# opar <- par()
# par(mfcol = c(n.prey, 1), mar = c(0,0,0,0))
# for(i in 1:n.prey){
# barplot(MTI.GOM[i,])
# }