太陽黒点数の推移プロット

国立天文台のサイト：

　　http://solarwww.mtk.nao.ac.jp/mitaka_solar1/data03/sunspots/number

から、太陽黒点の月平均個数を取得して、折れ線グラフを描く。

もうだいぶ前に作ったものなので、稚拙なスクリプトではあるが、いちおう動作する。

library(ggplot2)
library(ggthemes)

last.year <- as.numeric(substr(Sys.time(), 1, 4))

url <- "http://solarwww.mtk.nao.ac.jp/mitaka_solar1/data03/sunspots/number/rz"
df <- data.frame(NULL, NULL, NULL)

for(year in 1929:last.year) {
   url2 <- paste(url, as.character(year), ".txt", sep = "")
   rz <- file(url2, "r")
   while(length(str <- readLines(con = rz, 1)) != 0) {
      if(length(grep("Mon", str)) == 1) {
        vec <- unlist(strsplit(str, " "))
        vec <- vec[vec != ""]
        df <- rbind(df, as.numeric(tail(vec, n = 3)))
     }
   }
}

seq <- 0:(nrow(df) - 1)
yy <- as.character(seq %/% 12 + 1929)
mm <- as.character(seq %% 12 + 1)
date <- as.Date(paste(yy, "-", mm, "-01", sep = ""))
sunspots.df <- cbind(date, df)
names(sunspots.df) <- c("Date", "North", "South", "Total")

str <- substr(tail(sunspots.df$Date, 1), 1, 7)
substring(str, 5, 5) <- "/"
title <- paste("太陽黒点相対数の推移（月平均; 1929/01 ～ ", str, "）", sep ="")

g <- ggplot(sunspots.df, aes(x = Date, y = Total))
g <- g + geom_line(color = "darkblue", size = 0.1)
g <- g + geom_point(size = 0.5)
g <- g + xlab("year") + ylab("Number of Sun Spots")
g <- g + ggtitle(title)
g <- g + theme_economist()
g <- g + theme(axis.title=element_text(size=14))
print(g)