Sunday, December 8, 2013

Fukushima region radiation decrease

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.
Edit 9 December 2012: I got sent better translations from Ishigami, a Japanese reader, hence this post has been edited.

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
Maenori meeting place in Soso district     1.922496 2.074042 2.251524
Murohara community center in Soso district 3.171938 4.463582 7.529780
Sugiuchi meeting place in Soso district    1.764994 1.817172 1.872530
township village center                    2.602491 2.854096 3.159558

code

#【相双地方】室原公民館の測定値       
#704 
# Measurements [phase] of bi-district community center Murohara
# Measurements of Murohara community center in Soso district.
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
# Measured value of Sugiuchi meeting place in Soso district.
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
# Measurement of Maenori meeting place in Soso district.
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('Murohara community center in Soso district','Sugiuchi meeting place in Soso district',
            'Maenori meeting place in Soso district', '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('fuk2.png')
xyplot(doseus ~ day | loc,groups= status,data=all,type='l',
    scales=list(y=list(log=10)),
        par.strip.text=list(cex=.7),
    ylab=expression(mu * "Sv per hour" ))
dev.off()    

halflives <-     t(sa <- sapply(models,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

No comments:

Post a Comment