64dbe0e6c71787177bc46114306d6dffbfd57b2b
[covid19.git] / euromomo / Output_IA_ET_v42.R
1 ### Ambient temperatures ###
2 ET <- read.table(file=paste0(indir,"/temperature.txt"), header = T, sep =";", dec=".")
3 ET <- ET[with(ET, order(year, week)),]
4 ET$t = ifelse((ET$week %% 2) == 0, min(ET$temp, ET$ptmin, na.rm = TRUE), max(ET$temp, ET$ptmax, na.rm = TRUE))
5 ET$temp = ifelse(is.na(ET$temp), ET$ptemp, ET$temp)
6 ET$Uexecc <- with(ET, pmax(ptmax,temp))
7 ET$Lexecc <- with(ET, pmin(ptmin,temp))
8 ET$wk <- 1:nrow(ET)
9 ET$yw <- sprintf("%04iw%02d", ET$year, ET$week)
10
11 png(paste0(outdir,"/temperature_v4R.png"), width=1200, height=800, res=130)
12 par(mar=c(9,4,3,2))
13 with(ET, {
14 plot(x=wk, y=t, type="n", xlab=NA, ylab="Degree Celsius",
15 main="Ambient temperatures", bty="l", xaxt="n")
16 axis(1, at=wk[as.integer(substr(yw,6,7)) %in% c(1,27)],
17 labels=yw[as.integer(substr(yw,6,7)) %in% c(1,27)], las=2, cex.axis=0.9)
18 polygon(x=c(wk,rev(wk)), y=c(Uexecc, rev(ptmax)),
19 border=do.call(rgb, c(as.list(col2rgb("red")/255), alpha=0.5)),
20 col=do.call(rgb, c(as.list(col2rgb("red")/255), alpha=0.5)))
21 polygon(x=c(wk,rev(wk)), y=c(Lexecc, rev(ptmin)),
22 border=do.call(rgb, c(as.list(col2rgb("blue")/255), alpha=0.5)),
23 col=do.call(rgb, c(as.list(col2rgb("blue")/255), alpha=0.5)))
24 points(x=wk, y=ptemp, type="l", col="black")
25 points(x=wk, y=ptmin, type="l", col="blue")
26 points(x=wk, y=ptmax, type="l", col="red")
27 points(x=wk, y=temp, pch=19, cex=0.7)
28 mtext("Year-week", side=1, line=4.5)
29 legend("bottom", c("expected temperature", "max expected temperature", "min expected temperature"),
30 lty="solid", col=c("black", "red", "blue"), lwd=2,
31 horiz=TRUE, inset=c(0,-0.39), xpd=NA, bty="n", cex=1)
32 legend("bottom", c("observed temperature", "extreme heat", "extreme cold"),
33 col=c("black", do.call(rgb, c(as.list(col2rgb("red")/255), alpha=0.5)), do.call(rgb, c(as.list(col2rgb("blue")/255), alpha=0.5))), pch=c(19,15,15),
34 horiz=TRUE, inset=c(0,-0.46), xpd=NA, bty="n", cex=1, pt.cex=c(1,2,2))
35 })
36 dev.off()
37 rm(ET)
38
39
40 ### Influenza Activity ###
41
42 results4 <- read.table(
43 file = paste0(outdir, "/", country, "_output_v4", ifelse(IArest, "_IArestricted", ""), ".txt"),
44 header = TRUE, sep =";", dec=".")
45 results4 <- results4[with(results4, order(agegrp, year, week)),]
46 results4$wk <- unlist(lapply(table(results4$agegrp), seq))
47 results4$agegrp.labels <- ordered(results4$agegrp,
48 labels = c("0-4 years", "5-14 years", "15-64 years", "Aged 65", "Total"))
49 results4$yw <- sprintf("%04iw%02d", results4$year, results4$week)
50
51
52 png(paste0(outdir,"/IA_agegroups_v4R.png"), width=1200, height=1600, res=130)
53 par(mar=c(3,4,2,2), oma=c(5,2,2,0), mfrow=c(5,1))
54 for (a in 0:4)
55 with(results4[(results4$agegrp==a),], {
56 plot(x=wk, y=IA, xaxt="n", ylab=NA, xlab=NA, type="n", bty="l",
57 main=unique(results4$agegrp.labels)[a+1])
58 if (a==4) {
59 axis(1, at=wk[as.integer(substr(yw,6,7)) %in% c(1,27)],
60 labels=yw[as.integer(substr(yw,6,7)) %in% c(1,27)], las=2, cex.axis=0.9)
61 }
62 grid(nx=NA, ny=NULL)
63 abline(v=wk[as.integer(substr(yw,6,7)) %in% c(1,27)], col="lightgray", lty="dotted")
64 points(x=wk, y=IA, col="red", lwd=2, type="l")
65 })
66 mtext("Influenza Activity per week", side=2, outer=TRUE)
67 mtext("Influenza Activity by age group", side=3, outer=TRUE, adj=0, cex=1.2)
68 mtext("year-week", side=1, line=5)
69 dev.off()
70
71
72 for (a in 0:4) {
73 png(paste0(outdir,"/IA_agegroup_",a,"_v4R.png"), width=1200, height=800, res=130)
74 par(mar=c(6,4,4,2))
75 with(results4[(results4$agegrp==a),], {
76 plot(x=wk, y=IA, xaxt="n", ylab="Influenza Activity per week", xlab=NA, type="n", bty="l")
77 axis(1, at=wk[as.integer(substr(yw,6,7)) %in% c(1,27)],
78 labels=yw[as.integer(substr(yw,6,7)) %in% c(1,27)], las=2, cex.axis=0.9)
79 grid(nx=NA, ny=NULL)
80 abline(v=wk[as.integer(substr(yw,6,7)) %in% c(1,27)], col="lightgray", lty="dotted")
81 mtext(paste("Influenza Activity, age group:",unique(results4$agegrp.labels)[a+1]),
82 side=3, line=1, adj=0)
83 points(x=wk, y=IA, col="red", lwd=2, type="l")
84 mtext("year-week", side=1, line=4.5)
85 })
86 dev.off()
87 }
88 rm(a)