Deprecated: The each() function is deprecated. This message will be suppressed on further calls in /home/zhenxiangba/zhenxiangba.com/public_html/phproxy-improved-master/index.php on line 456
# par(family="HiraKakuProN-W3") # Mac専用 par(mgp=c(2,0.8,0)) X = read.csv("http://oku.edu.mie-u.ac.jp/~okumura/stat/data/20-50km.csv", fileEncoding="SJIS", na.strings = c("-","--"), as.is=TRUE) t = as.POSIXct(X$日時) plot(NULL, type="l", xaxt="n", xlab="", ylab="線量率（μSv/h）", # log="y", xlim=range(t), ylim=range(X[,2:9], na.rm=TRUE)) for (i in 1:8) points(t, X[,i+1], type="o", pch=i, col=rainbow(8)[i]) r = as.POSIXct(round(range(t), "days")) # axis.POSIXct(1, at=seq(r[1],r[2],by="day"), format="%b月%d日") axis.POSIXct(1, at=seq(r[1],r[2],by="day"), format="%d日") legend(as.POSIXct("2011-04-03"), 45, names(X[2:9]), col=rainbow(8), lty=1, pch=1:8, ncol=2) # 半減期8日の直線を引く x0 = as.POSIXct("2011-03-15 18:20:00") y0 = log2(X[t==x0,4]) curve(2^(y0-(x-as.numeric(x0))/(8*60*60*24)), add=TRUE) #--------------------------------------------------------------------- # 飯舘村だけ # par(family="HiraKakuProN-W3") # Mac専用 par(mgp=c(2,0.8,0)) p = 14 # 飯舘村役場 X = read.csv("http://oku.edu.mie-u.ac.jp/~okumura/stat/data/20-50km.csv", header=TRUE, fileEncoding="SJIS", na.strings = c("-","--"), as.is=TRUE) t = as.POSIXct(X$日時) t0 = as.POSIXct("2011-03-11 00:00:00") # t0 = as.POSIXct("2011-05-01 00:00:00") t1 = t[t >= t0 & !is.na(X[,p])] x1 = X[t >= t0 & !is.na(X[,p]), p] plot(t1, x1, type="n", xaxt="n", xlab="", ylab="線量率（μSv/h）") r = as.POSIXct(round(range(t1), "days")) axis.POSIXct(1, at=seq(r[1],r[2],by="day"), format="%b月%d日") # axis.POSIXct(1, at=seq(r[1],r[2],by="day"), format="%d日") for (i in seq(r[1], r[2], by="day")) { abline(v=i, col=gray(0.8)) } for (i in seq(0, max(x1), by=10)) { abline(h=i, col=gray(0.8)) } points(t1, x1, type="l", col="#0068b7") text(t1[length(t1)], x1[length(x1)], x1[length(x1)], pos=4, offset=0.1, cex=0.8) #--------------------------------------------------------------------- # 飯舘村だけ # par(family="HiraKakuProN-W3") # Mac専用 par(mgp=c(2,0.8,0)) p = 13 # 飯舘村役場 # p = 4 # 飯舘村長泥 X = read.csv("http://oku.edu.mie-u.ac.jp/~okumura/stat/data/20-50km5.csv", header=TRUE, fileEncoding="SJIS", na.strings = c("-","--"), as.is=TRUE) t = as.POSIXct(X$日時) # t0 = as.POSIXct("2011-03-11 00:00:00") t0 = as.POSIXct("2011-05-01 00:00:00") t1 = t[t >= t0 & !is.na(X[,p])] x1 = X[t >= t0 & !is.na(X[,p]), p] plot(t1, x1, type="n", xaxt="n", xlab="", ylab="線量率（μSv/h）") r = as.POSIXct(round(range(t1), "days")) axis.POSIXct(1, at=seq(r[1],r[2],by="day"), format="%b月%d日") # axis.POSIXct(1, at=seq(r[1],r[2],by="day"), format="%d日") for (i in seq(r[1], r[2], by="day")) { abline(v=i, col=gray(0.8)) } for (i in seq(0, max(x1), by=10)) { abline(h=i, col=gray(0.8)) } points(t1, x1, type="l") text(t1[length(t1)], x1[length(x1)], x1[length(x1)], pos=4, offset=0.1, cex=0.8) #--------------------------------------------------------------------- # 飯舘村の積算線量を求める X = read.csv("http://oku.edu.mie-u.ac.jp/~okumura/stat/data/20-50km5.csv", fileEncoding="SJIS", na.strings = c("-","--",""), as.is=TRUE) t = as.POSIXct(X$日時) p = 3 a = 0 n = length(X[,p]) cum = rep(NA, n) i = 1 cum[1] = 0 while (i < n && is.na(X[i,p])) { i = i + 1 } while (i < n) { j = i + 1 while (j < n && is.na(X[j,p])) j = j + 1 dt = (as.numeric(t[j]) - as.numeric(t[i])) / (60 * 60) x = (X[i,p] + X[j,p]) / 2 a = a + x * dt i = j cum[i] = a } # par(family="HiraKakuProN-W3") # Mac専用 par(mgp=c(2,0.8,0)) plot(t, X[,p], type="o", xaxt="n", xlab="", # log="y", pch=16, ylab="線量率（μSv/h）", col="#0068b7") r = as.POSIXct(round(range(t), "days")) axis.POSIXct(1, at=seq(r[1],r[2],by="day"), format="%d日") # par(family="HiraKakuProN-W3") # Mac専用 par(mgp=c(2,0.8,0)) plot(t, cum/1000, type="o", xaxt="n", xlab="", # log="y", pch=16, ylab="線量（mSv）", col="#0068b7") r = as.POSIXct(round(range(t), "days")) axis.POSIXct(1, at=seq(r[1],r[2],by="day"), format="%d日")