-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy path7d_all_met_sim.R
304 lines (256 loc) · 11.3 KB
/
7d_all_met_sim.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
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
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
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
### inserire dati completi Hg
setwd('C:\\Users\\gi\\Dropbox\\BlackSea2\\nuovi_script')
dati<-read.table('Dataset_Finale_new_no_bosforo.txt')
dati_mehg<-data.frame(dati$MeHg, dati$sigma)
str(dati_mehg)
medie_hg_pM<-c(1.86,2.127058824,1.7675,1.902,2.056666667,3.0875,2.810909091,3.714,3.713684211)
medie_mehg_pM<-c(0.116285714,0.1408125,0.120611111,0.108857143,0.157888889,0.767916667,0.551913043,0.751090909,0.737428571)
x1<-medie_hg_pM;x2<-medie_mehg_pM
sd1<-c(0.4911890,0.6179782,0.5060026,0.5832838,0.6269617,0.7147727,0.5266774,0.3722066,0.4352804)
sd2<-c(0.036545471,0.034533498,0.034428338,0.035267684,0.083861152,0.170347104,0.149375644,0.130289388,0.095633974)
setwd("C:/Users/gi/Desktop/")
tiff('LACombo17.tiff',
height=25, width=45, units='cm',
compression="lzw", res=300)
par(mfrow=c(1,4), oma = c(.1,9,3,9)+0.1,
mgp=c(4.2, 2, 0),mar = c(5.6,0.5,1,0)+ 0.1,
cex.axis=2.5, cex.lab=2.5,
cex.main=2.5)
names(dati_mehg)<-c('mehg','sigma')
### ALL WATER COL
## _______ MEHG ________________________##
plot(dati_mehg$mehg, dati_mehg$sigma,
bty='n',lwd=1.5,
cex=5,
ylim=c(17.24,10.5), lty=1, yaxt='n',
xlim=c(0, 1.3), pch=19, xaxt='n',
ylab=" ", xlab=" ", type="b",
col="#26261900", bg='#26261900',main=" ")
abline(h=15.64, lty=2, col='grey60')
abline(h=16.2, lty=2, col='grey60')
y2<-c(10.5,11.5,12.5,13.5,14.5,15.5,
16.5,17, 17.5)
y1<-c(11,12,13,14,15,16,17)
axis(2, at=y1,line=0, col="black", cex.axis=2.8 , las=2)
axis(2, at=y2,line=0,labels=F, col="black" , tck =-.02, cex.axis=2.8)
axx1<-c(0.1,0.3,0.5,0.7,0.9, 1.1, 1.3)
axx<-c(0,0.2,0.4,0.6,0.8,1,1.2)
text(1.2,10.7,'a', col='grey50', cex=5.3)
mtext(text="Modeled and observed MeHg concentrations",
cex=2.5, outer = T,
side=3,line=0)
mtext(expression(paste(sigma[theta]*' (kg m'^-3*')')),
side= 2, line=5.5,at=13.5, col="black", cex=2.5)
polygon(xx,yy,col='#b3db2544',border = '#b3db2544')
str(xx)
par(new=T)
plot(dati_mehg$mehg, dati_mehg$sigma, yaxt='n',
bty='n',lwd=1.8, cex=2.5, ylab='', xlab=' pM',
ylim=c(17.24,10.5), lty=1,
xlim=c(0, 1.3), pch=21,
type="p", xaxt='n',
col="#4d4d4d66", bg='#4d4d4d55',main=" ")
axis(1, at=aty2, labels =F,tick=T, las=1, cex.axis=2.8)
axis(1, at=aty2c, labels =aty2c,tick=T, las=1, cex.axis=2.8)
axis(1, at=maty2, labels =F,tick=T, las=1, cex.axis=2.8, tcl=-.2)
par(new=T)
plot(F_base$dissMehg_pM,type='l', prof, yaxt='n',
bty='n',lwd=1.5, cex=2, col='grey50',xaxt='n', #bg='#b3db2577',
ylim=c(17.24,10.5), lty=3, pch=23, ylab='',xlab='',
xlim=c(0, 1.3))
par(new=T)
plot(F_met_min$dissMehg_pM,type='l', prof, ylab='',
bty='n',lwd=1.5, cex=2,col='grey50',yaxt='n',
ylim=c(17.24,10.5), lty=3,xlab='',xaxt='n',
xlim=c(0, 1.3))
par(new=T)
plot(F_met_max$dissMehg_pM,type='l', prof, ylab='',
bty='n',lwd=1.5, cex=5,col='grey50',yaxt='n', xaxt='n',
ylim=c(17.24,10.5), lty=3, yaxt='n', xlab='',
xlim=c(0, 1.3))
par(new=T)
plot(F_demet_max$dissMehg_pM,type='l', prof, ylab='',
bty='n',lwd=1.5, cex=5,col='grey50',yaxt='n',xaxt='n',
ylim=c(17.24,10.5), lty=3, xlab='',
xlim=c(0, 1.3))
par(new=T)
plot(new1$dissMehg_pM,type='l', prof, ylab='',
bty='n',lwd=1.5, cex=5,col='grey50',yaxt='n',
ylim=c(17.24,10.5), lty=3, xlab='',xaxt='n',
xlim=c(0, 1.3))
text(0.58,12.3,expression(paste('k'[m])),col='grey20',font= 2,cex=4.3)
text(0.6,12.8,expression(paste('k'[dm])),col='grey20',cex=4.3)
text(0.6,15.9,expression(paste('k'[m])),col='grey20',cex=4.3)
text(1.1,15.9,expression(paste('k'[dm])),col='grey20',cex=4.3)
text(0.17,16.7,expression(paste('k'[m])),col='grey20',cex=4.3)
text(0.19,17.2,expression(paste('k'[dm])),col='grey20',cex=4.3)
##### SOL met + ubiquitus dem
## _______ MEHG ________________________## mar=c(4,5.5,2,0)
plot(medie_mehg_pM, prof, bty='n',lwd=1.5, cex=4,
ylim=c(17.24,10.5), lty=1, yaxt='n',
xlim=c(0, 1.3), pch=21, xaxt='n',
ylab=" ", xlab=" ", type="b",
col="#26261900", bg='#26261900',main=" ")
abline(h=15.64, lty=2, col='grey60')
abline(h=16.2, lty=2, col='grey60')
axis(1, at=aty2, labels =F,tick=T, las=1, cex.axis=2.8)
axis(1, at=aty2c, labels =aty2c,tick=T, las=1, cex.axis=2.8)
axis(1, at=maty2, labels =F,tick=T, las=1, cex.axis=2.8, tcl=-.2)
par(new=T)
plot(dati_mehg$mehg, dati_mehg$sigma, yaxt='n',
bty='n',lwd=1.8, cex=2.5, ylab='',xlab=' pM',
ylim=c(17.24,10.5), lty=1, xaxt='n',
xlim=c(0, 1.3), pch=21,
type="p",
col="#4d4d4d66", bg='#4d4d4d55',main=" ")
axis(1, at=aty2, labels =F,tick=T, las=1, cex.axis=2.8)
axis(1, at=aty2c, labels =aty2c,tick=T, las=1, cex.axis=2.8)
axis(1, at=maty2, labels =F,tick=T, las=1, cex.axis=2.8, tcl=-.2)
polygon(xx1,yy,col='#25b3db44',border ='#25b3db44')
par(new=T)
plot(L_base$dissMehg_pM,type='l', prof, yaxt='n',
bty='n',lwd=1.5, cex=2,col='grey50', xaxt='n', #bg='#b3db2577',
ylim=c(17.24,10.5), lty=3, pch=23, ylab='',xlab='',
xlim=c(0, 1.3))
par(new=T)
plot(L_met_min$dissMehg_pM,type='l', prof, ylab='',
bty='n',lwd=1.5, cex=2,col='grey50',yaxt='n',
ylim=c(17.24,10.5), lty=3,xlab='',xaxt='n',
xlim=c(0, 1.3))
par(new=T)
plot(L_met_max$dissMehg_pM,type='l', prof, ylab='',
bty='n',lwd=1.5, cex=4,col='grey50',yaxt='n', xaxt='n',
ylim=c(17.24,10.5), lty=3, yaxt='n', xlab='',
xlim=c(0, 1.3))
par(new=T)
plot(L_demet_min$dissMehg_pM,type='l', prof, ylab='',
bty='n',lwd=1.5, cex=4,col='grey50',yaxt='n', xaxt='n',
ylim=c(17.24,10.5), lty=3, xlab='',
xlim=c(0, 1.3))
par(new=T)
plot(L_demet_max$dissMehg_pM,type='l', prof, ylab='',
bty='n',lwd=1.5, cex=4,col='grey50',yaxt='n', xaxt='n',
ylim=c(17.24,10.5), lty=3, xlab='',
xlim=c(0, 1.3))
text(0.58,12.3,expression(paste('k'[m])),col='grey20',cex=4.3)
text(0.58,12.3,expression(paste('X')),col='#cc282899',cex=4.6)
text(0.58,12.3,expression(paste('X')),col='#cc282899',cex=4.6)
text(0.6,12.8,expression(paste('k'[dm])),col='grey20',cex=4.3)
text(0.55,15.9,expression(paste('k'[m])),col='grey20',cex=4.3)
text(1,15.9,expression(paste('k'[dm])),col='grey20',cex=4.3)
text(0.18,16.7,expression(paste('k'[m])),col='grey20',cex=4.3)
text(0.18,16.7,expression(paste('X')),col='#cc282899',cex=4.6)
text(0.18,16.7,expression(paste('X')),col='#cc282899',cex=4.6)
text(0.21,17.2,expression(paste('k'[dm])),col='grey20',cex=4.3)
why<-c(0,10, 25, 50, 75, 100, 250, 500, 2000)
text(1.2,10.7,'b', col='grey50',cex=5.3, outer=F)
###L demet mar=c(4,5.5,2,0)
plot(medie_mehg_pM, prof, bty='n',lwd=1.5, cex=4,
ylim=c(17.24,10.5), lty=1, yaxt='n',
xlim=c(0, 1.3), pch=21,xaxt='n',
ylab=" ", xlab=" ", type="b",
col="#26261900", bg='#26261900',main=" ")
axis(1, at=aty2, labels =F,tick=T, las=1, cex.axis=2.8)
axis(1, at=aty2c, labels =aty2c,tick=T, las=1, cex.axis=2.8)
axis(1, at=maty2, labels =F,tick=T, las=1, cex.axis=2.8, tcl=-.2)
abline(h=15.64, lty=2, col='grey60')
abline(h=16.2, lty=2, col='grey60')
par(new=T)
plot(dati_mehg$mehg, dati_mehg$sigma, yaxt='n',
bty='n',lwd=1.8, cex=2.5, ylab='',xlab=' pM',
ylim=c(17.24,10.5), lty=1, xaxt='n',
xlim=c(0, 1.3), pch=21,
type="p",
col="#4d4d4d66", bg='#4d4d4d55',main=" ")
polygon(xx2,yy,col='#db25b344',border = '#db25b344')
par(new=T)
plot(L1_base$dissMehg_pM,type='l', prof, yaxt='n',
bty='n',lwd=1.5, cex=2,col='grey50', xaxt='n', #bg='#b3db2577',
ylim=c(17.24,10.5), lty=3, pch=23, ylab='',xlab='',
xlim=c(0, 1.3))
par(new=T)
plot(L1_met_min$dissMehg_pM,type='l', prof, ylab='',
bty='n',lwd=1.5, cex=2,col='grey50',yaxt='n',
ylim=c(17.24,10.5), lty=3,xlab='',xaxt='n',
xlim=c(0, 1.3))
par(new=T)
plot(L1_met_max$dissMehg_pM,type='l', prof, ylab='',
bty='n',lwd=1.5, cex=4,col='grey50',yaxt='n', xaxt='n',
ylim=c(17.24,10.5), lty=3, yaxt='n', xlab='',
xlim=c(0, 1.3))
par(new=T)
plot(L1_demet_min$dissMehg_pM,type='l', prof, ylab='',
bty='n',lwd=1.5, cex=4,col='grey50',yaxt='n', xaxt='n',
ylim=c(17.24,10.5), lty=3, xlab='',
xlim=c(0, 1.3))
par(new=T)
plot(L1_demet_max$dissMehg_pM,type='l', prof, ylab='',
bty='n',lwd=1.5, cex=4,col='grey50',yaxt='n', xaxt='n',
ylim=c(17.24,10.5), lty=3, xlab='',
xlim=c(0, 1.3))
text(0.58,12.3,expression(paste('k'[m])),cex=4.3, col='grey20')
text(0.58,12.3,expression(paste('X')),col='#cc282899',cex=4.6)
text(0.58,12.3,expression(paste('X')),col='#cc282899',cex=4.6)
text(0.6,12.8,expression(paste('k'[dm])),col='grey20',cex=4.3)
text(0.6,12.8,expression(paste('X')),col='#cc282899',cex=4.6)
text(0.6,12.8,expression(paste('X')),col='#cc282899',cex=4.6)
text(0.55,15.9,expression(paste('k'[m])),col='grey20',cex=4.3)
text(0.98,15.9,expression(paste('k'[dm])),col='grey20',cex=4.3)
text(1.1,16.8,expression(paste('k'[m])),col='grey20',cex=4.3)
text(1.1,16.8,expression(paste('X')),col='#cc282899',cex=4.6)
text(1.1,16.8,expression(paste('X')),col='#cc282899',cex=4.6)
text(1.15,17.2,expression(paste('k'[dm])),col='grey20',cex=4.3)
text(1.15,17.2,expression(paste('X')),col='#cc282899',cex=4.6)
text(1.15,17.2,expression(paste('X')),col='#cc282899',cex=4.6)
y2<-c(10.5,11.5,12.5,13.5,14.5,15.5,
16.5,17, 17.5)
y1<-c(11,12,13,14,15,16,17)
text(1.2,10.7,'c', col='grey50',cex=5.3, outer=F)
####_____ COMBO mar=c(4,5.5,2,0) mar=c(4,5.5,2,0),
plot(medie_mehg_pM, prof, bty='n',lwd=1.5, cex=4,
ylim=c(17.24,10.5), lty=1, yaxt='n',
xlim=c(0, 1.3), pch=21, xaxt='n',
ylab=" ", xlab=" ", type="b",
col="#26261900", bg='#26261900',main=" ")
abline(h=15.64, lty=2, col='grey60')
abline(h=16.2, lty=2, col='grey60')
segments(x2-sd2,y,x2+sd2,y, col= '#262619', lwd=1.6)
epsilon <- 0.12
segments(x2-sd2,y-epsilon,x2-sd2,y+epsilon, lwd=1.6, col='#262619')
segments(x2+sd2,y-epsilon,x2+sd2,y+epsilon, lwd=1.6, col='#262619')
#title(expression(' All hypoteses \n modeled VS observed MeHg profiles'), line=1, cex=1.1)
text(1.2,10.7,'d', col='grey50',cex=5.3)
#lines(x, y, lwd = 2, col='#748929')
#points(x, y, lwd = 2, col='#748929', cex=4, pch=23)
par(new=T)
plot(medie_mehg_pM, prof, yaxt='n', xaxt='n',
bty='n',lwd=1.8, cex=3.9, ylab='',xlab=' pM',
ylim=c(17.24,10.5), lty=1,
xlim=c(0, 1.3), pch=21,
type="b",
col="#262619", bg='#26261944',main=" ")
polygon(xx,yy,col='#b3db2544',border = '#b3db25')
polygon(xx1,yy,col='#25b3db44',border = '#25b3db')
polygon(xx2,yy,col='#db25b344',border = '#db25b3')
axis(1, at=aty2, labels =F,tick=T, las=1, cex.axis=2.8)
axis(1, at=aty2c, labels =aty2c,tick=T, las=1, cex.axis=2.8)
axis(1, at=maty2, labels =F,tick=T, las=1, cex.axis=2.8, tcl=-.2)
polygon(xx,yy,col='#b3db2544',border = '#b3db25')
polygon(xx1,yy,col='#25b3db44',border = '#25b3db')
polygon(xx2,yy,col='#db25b344',border = '#db25b3')
text(0.9,12.5,'OL', cex=3.3, font=2, col='grey30')
text(1.1,15.9,'SOL', cex=3.3, font=2, col='grey30')
text(1,16.8,'AOL', cex=3.3, font=2, col='grey30')
y1<-c( 10.5,11.5, 13, 14.25, 15.64, 16.2, 16.55, 16.9, 17.24)
why<-c(0,10, 25, 50, 75, 100, 250, 500, 2000)
axis(4, at = y1, labels = why, tick = TRUE, cex.axis=2.8, las=2)
mtext('depth (m)', at=13.5, side=4,line=7.4,
cex=2.4, las=3)
par(new=T)
plot(medie_mehg_pM, prof, yaxt='n', xaxt='n',
bty='n',lwd=1.8, cex=3.9, ylab='',xlab=' ',
ylim=c(17.24,10.5), lty=1,
xlim=c(0, 1.3), pch=21,
type="b",
col="#262619", bg='#26261944',main=" ")
dev.off()