そらまめ君をRでWebスクレイピング（１）
題名のまんまの話です。
環境省大気汚染物質広域監視システム
（通称：そらまめ君）では、日本全国の様々な地域の大気汚染物質について
情報の公開を行っていますが、過去のものが1週間くらいしか遡れない。
大抵こういうのって、しばらくしてから欲しくなったりするわけで
それならば、ってことでWebスクレイピングで情報をゲットして
保存しつつ、ついでに図示もできたらな、という話。
library(rvest)
library(dplyr)
library(ggplot2)
library(maptools)
使ったパッケージはこんなところ。
で、自動化して定期的に実行したかったので
実行時の時刻を、以下の手法で取得
nowtime <- Sys.time()
nowtime2 <- as.POSIXlt(nowtime)
year <- nowtime2$year +1900
month <- nowtime2$mon + 1 #月データはなぜか0が１月でのスタート
day <- nowtime2$mday
hour <- nowtime2$hour-1
そらまめ君データが、どうも毎時30分の段階で更新されるようなので
自動化の実施時間が30分より後なら、上の
hour <- nowtime2$hour-1
は、最後が "-1" でいいけど、実施が30分より前なら "-2" がよろしいかと。
地域としては中部（東海）を想定していたので、後でデータのURLを
決定するのに必要な地域のID（今回は 04)を以下のように指定
block <- "04"
ちゃんと "04" と認識させたいので、数字じゃなくて文字として指定。
次に、得たい情報の時刻にのっとった文字列を作成
titleday <- paste(formatC(year, width=2, format = "d", flag = "0"),
formatC(month, width=2, format = "d", flag = "0"),
formatC(day, width=2, format = "d", flag = "0"),
sep="")
titledaytime <- paste(formatC(year, width=2, format = "d", flag = "0"),
formatC(month, width=2, format = "d", flag = "0"),
formatC(day, width=2, format = "d", flag = "0"),
formatC(hour, width=2, format = "d", flag = "0"),
formatC(block, width=2, format = "d", flag = "0"),
sep="")
地域IDを含むのと含まないのを作成。
なんでこんなことしたかというと、データのURLを
soraurl <- paste("http://soramame.taiki.go.jp/Gazou/Hyou/AllMst/", titleday,"/hb", titledaytime,".html", sep="")
で指定したいから。
でもって、これを使ってスクレイピング。
dat <- read_html(soraurl, encoding = "euc-jp")
ここで結構時間を使った！！
風向のデータが日本語なのだけど、それが文字化けしてしまって
解ればこんなものだけど、しばらく検索したりしてました。
tbls <- dat %>% html_table() #任意の時間のデータ
colnames(tbls1) <- c("測定局コード","市区町村名","測定局名称","SO2(ppm)","NO(ppm)",
"NO2(ppm)","NOx(ppm)","CO(ppm)","OX(ppm)","NHMC(ppmC)","CH4(ppmC)",
"THC(ppmC)","SPM(mg/m3)","PM2.5(ug/m3)","SP(mg/m3)","WD","WS(m/s)",
"TEMP(℃)","HUM(%)","局種別")
データを変換です。
呼び出したURLは項目名が無いので、自分で付け加え。
このtblsを適当にcsvで保存すれば、とりあえずデータのスクレイピングは完成。
この後は、自動化とかグラフ化とかあるけど
それはまた次へ（なぜなら、いまから忘年会だから）
信頼区間９５%が真値を含む様子の図
データの９５%信頼区間が
場合によっては真値を含まないケースを図示するので
以前、どこかを参考に書いたはずの絵が見つからないので
ggplot2で描いてみた。
もっと良い方法があるのかもしれないけど
ggplot2の練習もかねて、これで一応、良しとする。
dataset <- NULL
setmean <- 5
for (i in 1:100){
sampledata <- rnorm(100, setmean, 1)
testresu <- t.test(sampledata)
check <- NULL
if (testresu$conf.int[2] < setmean | testresu$conf.int[1] > setmean){
check <- 1
} else{
check <- 2
}
dataset <- rbind(dataset, c(i, testresu$estimate, testresu$conf.int[1], testresu$conf.int[2], check))
}
colnames(dataset) <- c("No", "ave", "min", "max", "check")
dataset <- as.data.frame(dataset)
ggplot() +
geom_line(aes(x=c(0,101), y =c(setmean, setmean)), size=2, colour="#00FF00") +
geom_point(data = subset(dataset, dataset$check == 1), aes(x=No, y=ave),
colour = "#FF0000") +
geom_errorbar(data = subset(dataset, dataset$check == 1),
aes(x=No, ymin = min, ymax = max), width =.2,
colour = "#FF0000") +
geom_point(data = subset(dataset, dataset$check == 2), aes(x=No, y=ave)) +
geom_errorbar(data = subset(dataset, dataset$check == 2),
aes(x=No, ymin = min, ymax = max), width =.2,
colour = "#000000") +
xlab("") + ylab("") +
theme(plot.title = element_text(size=20),
legend.title = element_text(size=15),
legend.text = element_text(size=15),
axis.text.x=element_blank(),
axis.text.y=element_text(size=15),
axis.title.y = element_text(size=15),
panel.border= element_rect(colour="#000000", fill=NA, size=1.5),
panel.background = element_rect(fill="#FFFFFF"),
panel.grid.major = element_blank())
R:検索結果がヒットしないのに、NAとか０にはならないとき
Rでのちょっと戸惑ったこと。
データからsubset()で検索値をとりだして
その値を使ってさらに足し算したいとき。
検索結果がヒットしなかったら
０を返してほしかったのだけど
そうじゃなくてlogical(0)なんて出てきて
これが煮ても焼いても食えなかった。
is.na()をやってもだめで
この数字自体に1をかけたり、1を割ったりして
何かエラーにしようとしても
numeric(0)とか出て、うまく扱えない。
で、length()で大きさを見ると、見事に数字の０に落ち着いた。
d <- subset(data, lon == basedata$lon[ i ] & lat == basedata$lat[ i ]
d2 <- subset(data, lon == basedata$lon[ j ] & lat == basedata$lat[ j ]
で、このdに該当する数字が見つからなかったときに
（d2はヒットして5だったとすると）
d + d2 とかしたい場合
そのままだと計算できないけど
length(d) + d2 だと、5が返される。
もっとも、dに数字がヒットした時に困るから
その前にifとかでちゃんと仕分けしておかないといけないけど。
強引かもしれないけど、なんとか解決。
後方流跡線をRで描いてみる・終
後方流跡線を描くとき
つまづいたのが、geom_line と geom_pathの違い
ずっとlineで描いていたけど、ある線が、おかしくなった。
（例のデータがどれだったのか、解らないので
例図がありません。そのうち見つけたら、更新）
なんだか、ギザギザになって明らかにおかしい。
で、細かく見てみると
lonのｘ軸が小さいほうから順番に連なっていることが発覚。
後方流跡線にありがちな、くるっと回っていくのが
うまく表現できていませんでした。
その解決のための変更がgeom_path
線を引くときは、どちらがいいのかちゃんと考えないと。
知ってる人は知っているんだろうけど
一瞬、絶望的な気分になった。
後方流跡線をRで描いてみる・続続
後方流跡線のデータを取り出す方法も分かった。
地図も描いた。
なので、あとはこれを
くっつけるだけ。
いろいろと捗るといいなぁと思って
データのファイル名を引数に関数にしてみた。
ちなみに、後方流跡線はNOAAのものです。
気象データはGDASの０．５度グリッド。
backtraje <- function(data) {
i <- 11
lat <- na.omit(as.numeric(unlist(strsplit(data[i]," ")))/1)[10]
lon <- na.omit(as.numeric(unlist(strsplit(data[i]," ")))/1)[11]
height <- na.omit(as.numeric(unlist(strsplit(data[i]," ")))/1)[12]
data3 <- as.numeric(c(lat, lon , height))
for (i in 12:length(data)){
lat <- na.omit(as.numeric(unlist(strsplit(data[i]," ")))/1)[10]
lon <- na.omit(as.numeric(unlist(strsplit(data[i]," ")))/1)[11]
height <- na.omit(as.numeric(unlist(strsplit(data[i]," ")))/1)[12]
data3 <- rbind(data3, as.numeric(c(lat, lon, height)))
}
colnames(data3) <- c("lat", "lon", "height")
data3 <- as.data.frame(data3)
#流跡線かくときはlineじゃなくてpath
ggmap3 <-
ggmap2 +
geom_path(data = data3, aes(x=lon, y = lat), size=0.5) +
geom_point(data = data3, aes(x=lon, y = lat, size=height), shape = 1, colour="red") +
labs(size="Height(m)") +
scale_radius(range=c(0,10),limits=c(0,max(data3$height))) +
xlab("") +
theme(axis.text.x = element_blank(),
legend.position = "top")
}
で、データを読み込む（データ名：12110103.txt)
そして関数で描画
data <- readLines("12110103.txt")
backtraje(data)
なんでGMTじゃないかというと
こうやって、高度の情報を組み込みたかったから。
（GMTでも、できるのかもしれないけど・・・・）
高度が高いところは円が大きく、低いところは小さい。
上空から見た感覚でイメージしてみた。
色で分けたのも作ったけど、解りにくくなったので
大きさで高度を表現。
これでサクサク見ていけば
なにか感じることが出てくるかもね～
後方流跡線をRで描いてみる・続
後方流跡線のデータを読み込んだら
あとはこっちのもの（なにが？）
地図を描いて、後方流跡線を乗っけるだけ。
今回、地図はNature Earthのものを使用しました。
荒めの国境データですが、今回にはこれで十分かと。
shapeファイルを読み込んで、ggplot2で描けるように
変換しました。
library(ggplot2)
library(maptools)
worldCountry10_shp <- readShapeLines("ne_10m_admin_0_countries.shp")
worldCountry10_map <- fortify(worldCountry10_shp)
ggmap1 <-
ggplot() +
geom_path(data = worldCountry10_map, aes(x = long, y = lat, group = group), colour="#999999") +
xlim(c(80, 150)) + ylim(c(20, 60)) +
coord_fixed(ratio=1.35/1)
coord_fixedは目視で設定した、個人的な趣味の数字です。
正解は、何なんだろう？
あと、中国は土地勘が無いので適当な都市を地図で表示。
pointdata <- read.csv("pointdata.csv", header=T)
> pointdata
name lat lon
1 Beijing 39.90449 116.3915
2 Shanghai 31.24787 121.4727
3 Hong_Kong 22.27838 114.1743
4 Ulan_Bator 47.92138 106.9055
5 Shenyang 41.80030 123.4338
6 Chongqing 29.57000 106.5800
ggmap2 <-
ggmap1 +
geom_point(data = pointdata,aes(x = lon, y = lat), size=4, shape=16, colour="blue") +
geom_text(data = pointdata, aes(x= lon+1, y = lat-1, label = name), size=5)
この段階で地図を描くと、こんな感じ。
あとは、これに後方流跡線のデータを載せるだけ。
それはまた、次に。
