Fukushima region radiation decrease

[This article was first published on Wiekvoet, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

Fukushima is still a topic which gets headlines, and somewhere in a comment was a link to actual and historical radiation data: Fukushima prefecture radioactivity measurement map.  It takes a bit of time to load but then you have a map of the Fukushima region with current radiation levels at measurement stations. Somewhere on that website there is data from the past year, this post examines the decrease at four of the stations.

Data

If you look at the map, there are loads of blue dots, indicating low radiation, and a line of yellow, orange and red dots, striking north east from the reactor. It is an interactive map. You can click on any dot and see some details for that region. And at the bottom of that is a link to a detail page. On that page there are, among other things, graphs such as this one.
直近90日間の空間線量の推移
One of these plots contains ‘transit of air dose rates in the most recent 365 days. It is a dynamic plot, the link to it contains the data, just like the plot above. It is pretty easy to extract the data from the link. Since getting the links was manual work, I decided to restrict my data to four dots. To discriminate between locations I copied some header and a number. For example 【相双地方】室原公民館の測定値 (Google translate: Measurements [phase] of bi-district community center Murohara) with number 704. The complete data with location information are at the bottom of the page.
There may be a more simple way to get the data, but since I don’t speak Japanese, I cannot say either way.

Analysis

Since the data looks pretty choppy, I wanted a robust model. To quote the robust task view:  “robust depends on robustbase, the former providing convenient routines for the casual user where the latter will contain the underlying functionality, and provide the more advanced statistician with a large range of options for robust modeling.” Regarding robust I am certainly the casual user. Proportional decrease by using a linear model on log transformed data seemed an appropriate modelling strategy. 
The plot below contains the data and the predictions. For me, the predictions seem just like I would make them when ignoring some of the crazy jumps.
The data themselves, show a clear correlation. The same drop mid January, a jump in September, though I do not know an explanation.
Half lives can be estimated from the models. Two of them are around 2. Cesium-134 has a half live of 2.04 years and has apparently been released, it is easy to speculate Cesium-134 is part of the source of this radiation.
                                  Lower CI     mean Upper CI
community center Murohara         3.171937 4.463582 7.529786
meeting place to the power before 1.922496 2.074042 2.251525
Sugiuchi meeting place            1.764994 1.817172 1.872530
township village center           2.602469 2.854075 3.159538

code

#【相双地方】室原公民館の測定値       
#704 
# Measurements [phase] of bi-district community center Murohara
ys704 <- "http://chart.apis.google.com/chart?chd=t:433,432,433,410,431,432,426,430,441,438,435,432,435,416,422,431,422,422,407,425,405,421,419,424,421,418,415,428,412,402,414,404,413,412,401,402,406,413,424,276,310,334,352,357,370,389,396,398,402,405,404,396,397,405,404,404,407,417,419,411,413,411,332,409,398,400,409,388,402,402,403,394,387,400,408,394,394,389,405,384,391,387,399,404,407,416,407,403,409,409,412,423,428,415,411,408,416,431,403,411,415,420,436,426,420,407,416,414,411,411,404,402,415,406,405,393,420,403,389,401,403,390,397,403,410,388,393,388,395,412,403,408,413,391,395,383,377,392,401,399,402,386,394,392,405,403,387,392,394,404,406,411,391,406,401,408,387,389,392,401,399,391,390,405,402,411,408,409,408,407,404,411,413,417,405,390,399,393,392,401,402,405,411,410,405,414,414,411,405,404,417,395,396,401,403,396,406,398,401,403,407,408,388,390,393,396,398,402,405,392,404,409,406,410,405,419,423,420,417,400,402,398,400,387,395,389,397,399,386,395,386,396,396,386,392,383,392,390,380,379,391,392,385,390,400,405,408,416,417,415,415,415,417,417,419,419,420,405,411,405,401,406,408,411,403,403,408,417,402,401,395,405,395,389,394,395,381,389,392,395,399,390,396,388,370,380,383,383,405,391,387,384,389,382,369,374,377,384,379,382,375,374,372,365,366,374,379,388,378,378,373,366,368,343,355,353,358,361,347,351,334,346,340,334,327,332,334,343,331,337,335,338,344,331,333,336,341,331,336,338,322,320,322,331,330,330,331,335,326,331,325,325,327,328,331,325,331,324,324,320,323,322,327,324,323&chxl=0:|365|330|300|270|240|210|180|150|120|90|60|30|1|2:|%CE%BCSv%2Fh&chxp=0|2,0&chxr=0,0,100|1,0,50&chxs=0,676767,9,0,lt,676767|1,676767,9,0,l,676767&chxtc=0,3&chxt=x,y,t&chs=300x150&cht=lc&chco=0067B7&chds=0,5000&chls=1&chma=1,1,1,1&chg=8.333,10,1,1"
#【相双地方】杉内集会所の測定値       
#628 
# Measured value of the bi-phase region] Sugiuchi meeting place
ys628 <- "http://chart.apis.google.com/chart?chd=t:280,281,281,269,278,279,277,279,282,278,281,273,275,275,274,276,267,267,269,270,271,272,261,267,269,269,270,272,270,266,269,270,270,268,270,269,268,270,268,183,183,197,201,203,208,216,223,234,243,246,253,250,251,254,255,255,256,262,265,266,258,263,201,242,245,250,255,256,255,258,259,252,253,256,258,260,259,258,259,256,257,258,261,261,262,258,263,262,263,264,264,268,269,265,265,265,265,268,265,265,267,268,267,266,265,265,266,265,263,261,262,247,253,254,253,247,253,244,235,246,249,232,242,247,250,239,245,240,246,251,250,251,254,246,244,234,234,240,243,233,242,230,239,242,244,243,234,241,244,247,248,248,246,247,248,248,227,231,236,238,240,236,239,243,241,242,242,246,245,245,245,245,247,247,239,221,231,230,227,233,235,237,238,238,227,236,237,238,226,224,232,218,223,226,225,220,224,216,217,224,227,230,217,218,220,220,224,226,230,218,224,227,227,228,222,230,231,235,236,222,223,222,224,225,217,218,222,223,214,218,211,216,215,210,214,209,214,214,207,209,214,215,210,211,217,221,223,228,229,231,232,232,235,235,237,238,238,226,227,227,225,228,230,233,221,225,226,231,222,223,223,228,215,214,215,216,215,216,217,219,221,213,218,210,204,209,211,212,215,217,215,214,216,203,205,207,209,211,208,206,201,205,202,197,197,201,204,207,206,205,206,206,206,190,194,196,197,198,187,190,185,188,183,185,186,190,191,193,194,195,195,195,196,192,193,195,194,193,194,194,188,190,191,192,192,191,193,193,194,194,194,194,194,194,195,189,190,191,190,190,191,190,191,188,188&chxl=0:|365|330|300|270|240|210|180|150|120|90|60|30|1|2:|%CE%BCSv%2Fh&chxp=0|2,0&chxr=0,0,100|1,0,50&chxs=0,676767,9,0,lt,676767|1,676767,9,0,l,676767&chxtc=0,3&chxt=x,y,t&chs=300x150&cht=lc&chco=0067B7&chds=0,5000&chls=1&chma=1,1,1,1&chg=8.333,10,1,1"
# 【相双地方】前乗集会所の測定値   
#741 
# Measurement of the meeting place to the power [bi-phase region] before
ys741 <- 'http://chart.apis.google.com/chart?chd=t:114,113,112,110,98,101,102,105,109,112,115,114,114,114,105,105,104,105,102,92,91,93,92,95,99,100,100,103,99,98,97,100,99,98,99,98,98,99,100,77,56,56,56,51,48,49,50,48,48,49,49,46,43,45,45,46,46,49,66,67,68,70,59,68,64,65,68,66,67,64,66,63,61,62,65,64,65,63,66,64,57,58,59,64,68,75,77,78,81,84,90,96,99,98,98,97,100,101,99,100,101,103,104,103,105,101,103,102,102,102,100,100,103,102,101,99,102,100,94,97,99,96,96,98,99,98,98,95,97,101,100,101,102,97,98,92,64,90,96,96,97,94,97,98,98,99,95,96,97,98,99,100,98,100,101,102,95,98,99,100,100,97,98,100,98,99,99,100,100,100,100,101,101,101,96,93,96,96,95,97,98,99,100,99,93,96,98,97,97,96,97,90,91,93,93,90,92,92,93,93,94,94,87,90,90,90,91,92,92,88,90,90,91,91,88,91,92,92,92,87,89,87,90,85,86,88,89,90,87,90,87,88,88,85,86,85,85,86,84,84,87,88,85,84,88,89,89,91,91,92,92,93,94,94,94,94,95,88,90,91,89,90,87,87,87,89,90,91,87,86,87,87,87,84,86,86,83,85,86,86,87,87,88,84,83,85,85,86,88,88,87,86,87,86,83,84,85,86,85,85,83,83,83,81,81,82,84,85,85,86,84,85,85,79,80,80,81,82,77,79,77,78,78,78,78,79,79,80,79,79,80,80,81,79,79,80,79,79,79,78,76,76,76,78,78,77,78,79,78,79,78,78,78,78,75,77,77,76,76,76,76,76,77,76,76&chxl=0:|365|330|300|270|240|210|180|150|120|90|60|30|1|2:|%CE%BCSv%2Fh&chxp=0|2,0&chxr=0,0,100|1,0,50&chxs=0,676767,9,0,lt,676767|1,676767,9,0,l,676767&chxtc=0,3&chxt=x,y,t&chs=300x150&cht=lc&chco=0067B7&chds=0,5000&chls=1&chma=1,1,1,1&chg=8.333,10,1,1'
# 【相双地方】町区集落センターの測定値  
#1376 
# Measured value of the bi-phase region] township village center
ys1376 <- 'http://chart.apis.google.com/chart?chd=t:1023,1010,991,987,953,1027,1007,1008,1006,957,963,980,1000,899,1003,994,920,995,886,1008,884,1003,910,915,918,904,977,1006,881,993,888,984,908,976,986,861,869,995,986,612,651,683,711,709,728,733,770,813,834,853,849,851,981,969,994,968,963,986,978,988,967,965,742,967,833,855,946,973,978,853,970,836,821,840,975,835,841,977,854,971,970,824,982,849,984,885,971,861,976,985,883,902,985,878,865,983,986,992,846,867,987,901,909,987,887,984,966,875,983,957,849,841,953,855,849,814,875,956,827,844,849,825,830,849,865,826,838,825,834,869,853,872,949,854,834,1003,1004,1041,1042,1012,1024,865,989,1009,1015,1003,981,997,1024,902,1029,1027,1005,899,894,908,975,986,1000,1009,1014,996,985,998,1003,999,991,1008,1020,988,1003,1016,919,921,984,847,872,981,980,966,996,991,898,890,977,988,984,979,875,874,989,861,862,873,965,858,867,846,844,860,960,972,845,847,854,855,869,869,879,850,865,867,872,877,859,881,873,884,885,849,853,850,849,828,837,829,841,849,826,840,823,833,839,818,834,814,829,830,816,813,827,831,825,822,837,845,856,866,930,860,861,862,869,873,872,878,881,847,851,855,847,850,854,861,846,847,855,870,862,858,860,870,849,840,837,841,818,833,839,844,851,834,848,828,803,814,820,821,854,834,836,826,837,811,796,799,809,825,817,813,801,803,797,784,790,808,813,820,787,795,780,772,774,743,739,743,750,754,729,817,791,793,795,797,769,782,784,808,793,805,790,799,812,783,779,789,805,779,793,800,747,752,748,771,768,769,772,785,766,777,782,779,774,767,777,763,761,758,759,747,764,758,763,746,745&chxl=0:|365|330|300|270|240|210|180|150|120|90|60|30|1|2:|%CE%BCSv%2Fh&chxp=0|2,0&chxr=0,0,100|1,0,50&chxs=0,676767,9,0,lt,676767|1,676767,9,0,l,676767&chxtc=0,3&chxt=x,y,t&chs=300x150&cht=lc&chco=0067B7&chds=0,5000&chls=1&chma=1,1,1,1&chg=8.333,10,1,1'

linktodata <- function(ys,i) {
  ys <- gsub('.*t:','',ys)
  ys <- gsub('&chxl.*','',ys)
  dose <- scan(textConnection(ys),sep=',')
  data.frame(dose=dose,station=i,relday = -(length(dose):1))
}

locations <-data.frame(station=c(704,628,741,1376),
    loc=factor(c(‘community center Murohara’,’Sugiuchi meeting place’,’meeting place to the power before’, ‘township village center’)))

r1 <- rbind(linktodata(ys704,704),
    linktodata(ys628,628),
    linktodata(ys1376,1376),
    linktodata(ys741,741))
r1$doseus <- r1$dose/100
r1$day <- as.Date("05-12-2013",format='%d-%m-%Y')+r1$relday
Sys.setlocale(“LC_TIME”,’C’)

r1 <- merge(r1,locations)
r1$status <- 'Observed'
library(robust)

models <- lapply(levels(r1$loc),function(st) {
      r2 <- r1[r1$loc==st,]
      lmrob(log10(doseus) ~ relday,data=r2)
    })

topred <- r1[r1$station==r1$station[1],c('relday','day')]

preds <- lapply(1:nlevels(r1$loc),function(i) {
      preds <- topred
      preds$doseus <- 10^predict(models[[i]],preds)
      preds$status <- 'Predicted'
      preds$loc <- levels(r1$loc)[i]
      preds
    }
)
preds <- do.call(rbind,preds)
preds <- merge(preds,locations)
all <- rbind(subset(r1,select=-dose),preds)

library(lattice)    
png(‘fuk1.png’)
xyplot(doseus ~ day | loc,groups= status,data=all,type=’l’,
    scales=list(y=list(log=10) ),
    ylab=expression(mu * “Sv per hour” ))
dev.off()    

halflives <-     t(sa <- sapply(la,function(x) {
          point <- coef(x)[2]
          sd <- sqrt(vcov(x)[2,2])
          yy <- c(log10(.5)/c(`0.025`=point-1.96*sd,mean=point,`0.975`=point+1.96*sd))/365
          names(yy) <- c('Lower CI','mean','Upper CI')
          yy
        }))
rownames(halflives) <- levels(r1$loc)
halflives

To leave a comment for the author, please follow the link and comment on their blog: Wiekvoet.

R-bloggers.com offers daily e-mail updates about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

Never miss an update!
Subscribe to R-bloggers to receive
e-mails with the latest R posts.
(You will not see this message again.)

Click here to close (This popup will not appear again)