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)