-
Notifications
You must be signed in to change notification settings - Fork 1
/
analysis-martin_curve_calc.R
executable file
·109 lines (83 loc) · 3.27 KB
/
analysis-martin_curve_calc.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
#!/usr/bin/env Rscript
library("plyr")
library("reshape2")
library("ggplot2")
library("gridExtra")
library("stringr")
library("viridis")
library("RColorBrewer")
source("functions.R")
Spectral <- brewer.pal(8, "Spectral")
# load biomes
load('data/biomes_chl_mixedlayer.Rdata')
N_BIOME=4
seafloor <- read.csv("data/woa_seafloor_depth.csv")
# Import Export to depth data for all three biomes
exportres = data.frame()
uores = data.frame()
for (C_ID in 1:3){
CASES = c("0-baseline","1-high_biomass","2-low_biomass")
CASE = CASES[C_ID]
load(file=str_c("data/",CASE,"/baseline/allexport.Rdata")) # units are all g C m^-3 y^-1
load(file=str_c("data/",CASE,"/baseline/upperocean_ens_mean.Rdata"))
allexport$case <- CASE
ens_mean$case <- CASE
exportres <- rbind(exportres, allexport)
uores <- rbind(uores, ens_mean)
}
exportres$case = factor(exportres$case, levels=c("2-low_biomass", "0-baseline", "1-high_biomass"))
exportres <- join(exportres, biomes_df, by=c("lat","lon"))
sumres=exportres[exportres$taxon=="sum",]
exprt=sumres[,c("lat",'lon','export','case','biome')]
seq_d=sumres[,c("lat",'lon','seq_d','case','biome')]
seaflr=sumres[,c("lat",'lon','seafloor','case','biome')]
exprt <- join(exprt, seafloor, by=c("lat","lon"))
exprt$depth = ifelse(exprt$depth > 100, 100, exprt$depth)
seq_d <- join(seq_d, seafloor, by=c("lat","lon"))
seq_d$depth = ifelse(seq_d$depth > 1000, 1000, seq_d$depth)
seaflr = join(seaflr, seafloor, by=c("lat","lon"))
names(seaflr) = names(seq_d) = names(exprt)
db=rbind(exprt, seq_d, seaflr)
db$log10depth=log10(db$depth)
db$log10export=log10(db$export)
db = db[is.finite(db$log10export) & is.finite(db$log10depth), ]
m = lm(-log10depth~log10export, data=db)
print('Results from all sinking POC:')
print(summary(m))
# b = 0.18 +/- 0.003
###
exprt=sumres[,c("lat",'lon','export_cf','case','biome')]
seq_d=sumres[,c("lat",'lon','seq_d_cf','case','biome')]
seaflr=sumres[,c("lat",'lon','seafloor_cf','case','biome')]
exprt <- join(exprt, seafloor, by=c("lat","lon"))
exprt$depth = ifelse(exprt$depth > 100, 100, exprt$depth)
seq_d <- join(seq_d, seafloor, by=c("lat","lon"))
seq_d$depth = ifelse(seq_d$depth > 1000, 1000, seq_d$depth)
seaflr = join(seaflr, seafloor, by=c("lat","lon"))
names(seaflr) = names(seq_d) = names(exprt)
db=rbind(exprt, seq_d, seaflr)
db$log10depth=log10(db$depth)
db$log10export=log10(db$export)
db = db[is.finite(db$log10export) & is.finite(db$log10depth), ]
m = lm(-log10depth~log10export, data=db)
print('Results from carcass flux only:')
print(summary(m))
# b = 0.15 +/- 0.003
###
exprt=sumres[,c("lat",'lon','export_eg','case','biome')]
seq_d=sumres[,c("lat",'lon','seq_d_eg','case','biome')]
seaflr=sumres[,c("lat",'lon','seafloor_eg','case','biome')]
exprt <- join(exprt, seafloor, by=c("lat","lon"))
exprt$depth = ifelse(exprt$depth > 100, 100, exprt$depth)
seq_d <- join(seq_d, seafloor, by=c("lat","lon"))
seq_d$depth = ifelse(seq_d$depth > 1000, 1000, seq_d$depth)
seaflr = join(seaflr, seafloor, by=c("lat","lon"))
names(seaflr) = names(seq_d) = names(exprt)
db=rbind(exprt, seq_d, seaflr)
db$log10depth=log10(db$depth)
db$log10export=log10(db$export)
db = db[is.finite(db$log10export) & is.finite(db$log10depth), ]
m = lm(-log10depth~log10export, data=db)
print('Results from fecal matter only:')
print(summary(m))
# b = 0.12 +/- 0.002