-
Notifications
You must be signed in to change notification settings - Fork 5
/
Copy pathJARAplotting.R
174 lines (146 loc) · 6.29 KB
/
JARAplotting.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
#------------------------------------
# Examples for easy plotting in JARA
# Henning Winker ([email protected])
# JRC - European Commission
#------------------------------------
#Test
library(JARA)
data("jaradata")
dat = jrdat
# Set working directory
setwd("C:/Work/Research/GitHub/JARA.Rcheck")
# create subfolder
output.dir = "JARAplotting"
dir.create(file.path(getwd(),output.dir),showWarnings = F)
# Example Smoothhound Shark
shsh.input = build_jara(I=jrdat$SmoothhoundShark$I,se=jrdat$SmoothhoundShark$SE,model.type = "relative",assessment = "Smoothhound",GL=13.7)
# check trawl survey indices
jrplot_indices(shsh.input)
# Now save as png to output dir
jrplot_indices(shsh.input,as.png=T,output.dir = output.dir)
# fit jata
shsh.fit = fit_jara(shsh.input,quickmcmc=F,do.ppc = T) # quick test run
# now plot overall fit
jrplot_trjfit(shsh.fit)
# you can save this also as jpg
dev.print(jpeg,paste0(output.dir,"/shsh.trj.jpg"), width = 5, height = 4.5, res = 300, units = "in")
# or use the inbuilt save option
jrplot_trjfit(shsh.fit,as.png = T,output.dir = output.dir)
# check out the plotting options
?jrplot_trjfit
# another useful function is jrpar() to setup the graphic parameters
jrpar(mfrow=c(1,1),labs=T,plot.cex = 0.8)
jrplot_trjfit(shsh.fit)
dev.print(jpeg,paste0(output.dir,"/shsh.trj.jpg"), width = 5, height = 4.5, res = 300, units = "in")
# Lets look at Striped Marlin (relative)
# create JARA input
mls.input = build_jara(I=dat$StripedMarlin_IO_CPUE$I,se=dat$StripedMarlin_IO_CPUE$SE,model.type = "relative",assessment = "MLS",GL=5.5,variance.weighting = "model")
mls.fix = build_jara(I=dat$StripedMarlin_IO_CPUE$I,se=dat$StripedMarlin_IO_CPUE$SE,model.type = "relative",assessment = "MLS",GL=5.5,variance.weighting = "equal")
# Check input
jrplot_indices(mls.input)
# fit JARA model
mls.fit = jara = fit_jara(mls.input,quickmcmc=T,do.ppc=T) # quick test run
fit.fix = fit_jara(mls.fix,quickmcmc=T,do.ppc=T)
load("MLS_s1_posteriors.rdata")
jrplot_fits(mls.fit)
jrplot_PPC(mls.fit,joint.ppc = F)
jrplot_PPC(fit.fix,joint.ppc = F)
ppc = ap$PPC
jrpar(mfrow=c(1,1))
plot(ppc[,2:3])
abline(0,1,col=2)
# Test Colour option
jrplot_trjfit(mls.fit,cols=rainbow(6))
# now if we want to show two examples in the same plot we use the option add=T
jrpar(mfrow = c(1,2),plot.cex = 0.8)
jrplot_indices(shsh.input,add=T)
mtext(c("a)"),side=3,outer=F,cex=1.1,line=0.5,adj=-0.1)
jrplot_indices(mls.input,add=T)
mtext(c("b)"),side=3,outer=F,cex=1.1,line=0.5,adj=-0.1)
# save to file
dev.print(png,paste0(output.dir,"/Indices2x1.png"), width = 8, height = 4, res = 300, units = "in")
# note the legend may distorted if using R studio to avoid this use of
#e.g. windows(width=8,height=4) or quartz()
# We can do the same with fits
jrpar(mfrow = c(1,2),plot.cex = 0.8)
jrplot_trjfit(shsh.fit,add=T)
mtext(c("a)"),side=3,outer=F,cex=1.1,line=0.5,adj=-0.1)
jrplot_trjfit(mls.fit,add=T)
mtext(c("b)"),side=3,outer=F,cex=1.1,line=0.5,adj=-0.1)
dev.print(png,paste0(output.dir,"/trj2x1.png"), width = 8, height = 4, res = 300, units = "in")
# Save Fits for Readme file illustration
jrplot_fits(mls.fit,as.png = T,output.dir = output.dir)
jrplot_logfits(mls.fit,as.png = T,output.dir = output.dir)
# Now take a look at the census data for Penguin
# African Penguin
# create JARA input
ap.input = build_jara(I=dat$Afr_penguin$I,assessment="AfrPenguin",model.type = "census",GL=9.9,fixed.obsE = 0.15)
# Check input indices
jrplot_indices(ap.input)
# fit JARA model
ap.fit = fit_jara(ap.input,quickmcmc=F,do.ppc=T)
# and check the fits and population trajectory
jrpar(mfrow = c(1,2),plot.cex = 0.7)
jrplot_trjfit(ap.fit,add=T)
mtext(c("a)"),side=3,outer=F,cex=1.1,line=0.5,adj=-0.1)
jrplot_poptrj(ap.fit,add=T)
mtext(c("b)"),side=3,outer=F,cex=1.1,line=0.5,adj=-0.1)
dev.print(png,paste0(output.dir,"/AP2x1.png"), width = 8, height = 4, res = 300, units = "in")
# we can quickly compare to a relative abudance example
jrpar(mfrow = c(1,2),plot.cex = 0.7)
jrplot_trjfit(shsh.fit,add=T)
mtext(c("a)"),side=3,outer=F,cex=1.1,line=0.5,adj=-0.1)
jrplot_poptrj(shsh.fit,add=T)
mtext(c("b)"),side=3,outer=F,cex=1.1,line=0.5,adj=-0.1)
dev.print(png,paste0(output.dir,"/SHSH_trend.png"), width = 8, height = 4, res = 300, units = "in")
# Probably the most important JARA result is
# the change over time
jrpar(mfrow = c(1,2),plot.cex = 0.7)
jrplot_changes(ap.fit,add=T)
mtext(c("a)"),side=3,outer=F,cex=1.1,line=0.5,adj=-0.1)
jrplot_iucn(ap.fit,add=T)
mtext(c("b)"),side=3,outer=F,cex=1.1,line=0.5,adj=-0.1)
dev.print(png,paste0(output.dir,"/APstatus2x1.png"), width = 8, height = 4, res = 300, units = "in")
#-----------------------------------------------------
# An alternative option is to read in the JARA plots from a completed assessment folder
# Examples of joining JARA outputs into multiplots
#-----------------------------------------------------
# All plot will be save in the respective species assessment folders
# For example run Yellowspotted Skate
# create a folder
assessment = "YellowspottedSkate"
newoutput.dir = assessment
dir.create("YellowspottedSkate",showWarnings = F)
ysk=build_jara(I=dat$Yellowspot_skate$I,se=dat$Yellowspot_skate$SE,model.type = "relative",assessment = assessment,GL=12)
# fit
fit.ysk = fit_jara(ysk)
# Plot all output in newoutput.dir
jara_plots(fit.ysk,as.png = T,output.dir = newoutput.dir)
#----------------------------------------------------------------
# Produce IUCN Supplement Information "style" Plot
#----------------------------------------------------------------
DIMs=c(5,4) # Plot dimension
sp=8 # Select species
Plot = c("IUCN_SoupfinShark")
DIMs=c(2*5,2*4)/1.5
# Choos plot types
plots = c("TrjFits_","PopTrend_","AnnualRate_","IUCNplot_")
j=1
runs = 1 # run name
# layout the plots into a matrix columns, by row
library(png)
par(mar=rep(0,4),omi= c(0, 0, 0, 0),cex=1.3) # no margins
layout(matrix(1:4, ncol=2, byrow=TRUE))
for(j in 1:length(plots)){
run = runs[1]
# plot path
get_plot = file.path(getwd(),assessment,paste0(plots[j],assessment,"_s1.png"))
# load image
img <- readPNG(paste0(get_plot))
# do the plotting
plot(NA,xlim=0:1,ylim=0:1,xaxt="n",yaxt="n",bty="n")
rasterImage(img,0,0,1,1)
legend("topleft",paste0("(",letters[j],")"),bty="n",cex=1.2,x.intersp = -0.55,y.intersp=-0.2)
}
# Write plot
dev.print(png,paste0(output.dir,"/",assessment,"_status.png"), width = DIMs[1], height = DIMs[2], res = 300, units = "in")