気象庁の二酸化炭素濃度データ

Python版気象庁の二酸化炭素濃度データのR（tidyverse）版である。

library(tidyverse)

df = read_csv("co2_monthave_ryo.csv", locale=locale(encoding="CP932"))

これだけでは「年」の列の末尾に説明の文章が入ってしまう。そこでは「月」は NA になるので，drop_na(df, 月) で末尾が外せるが，数値であるべきところが文字列になるという問題が残る。

そこで，col_types="iidc" とすることにより，各列のデータ型をそれぞれinteger・integer・double・character（string）と指定すれば，これ以外の型の内容はすべて NA になる。

df = read_csv("co2_monthave_ryo.csv",
              locale=locale(encoding="CP932"), col_types="iidc")

NA は描画されないので特に取り除く必要はないが，drop_na() で取り除くほうがわかりやすい。

df = drop_na(df, 年)
t = as.Date(paste0(df[[1]], "-", df[[2]], "-15"))

as.Date() は as.POSIXct() でもよい。

# svg("co2_monthave_ryo.svg", width=7, height=5)
# par(mgp=c(2, 0.8, 0))
plot(t, df[[3]], type="l", xlab="year", ylab=expression(paste(CO[2], " (ppm)")))
# dev.off()
綾里の二酸化炭素濃度

あるいは，時系列を表す ts クラスに変換するほうが扱いが楽かもしれない：

co2 = ts(df[[3]], c(1987,1), frequency=12)
plot(co2)

これと，stl（Seasonal Decomposition of Time Series by Loess）で seasonal と trend に分解したときの trend とを重ね書きしたグラフを描く。欠測値は前後の平均とする。

# svg("co2_trend_ryo.svg", width=7, height=5)
# par(mgp=c(2, 0.8, 0))
plot(co2)
x = which(is.na(co2))
co2[x] = (co2[x-1] + co2[x+1]) / 2
points(stl(co2, "periodic")$time.series[,2], type="l")
# dev.off()
綾里の二酸化炭素濃度

Last modified: