RgoogleMapsを使った放射線地図

RgoogleMaps はGoogle MapsをRから使うためのパッケージです。

options(repos="http://cran.ism.ac.jp") # 適当にCRANミラーを指定
update.packages()                      # 念のため
install.packages("RgoogleMaps")        # RgoogleMapsをインストール
install.packages("RColorBrewer")       # RColorBrewerをインストール

Googleのスタティックな地図を取得する関数 GetMap() で言語指定するための hl オプションが壊れているので修正した GetMap.R を置いておきます。作者には連絡済みで,次のバージョンでは直るはずです。

@nnistar さんのまとめられた福島県のデータを全部取得して放射線地図を描いてみましょう。CSVで取得したら一つにつなげて clean.rb のようなものできれいにしておきます。42万数千行ありますので,時間順にソートして,新しいデータが上に描かれるようにします。念のため ここに fukushima.csv という名前で置いておきました(データそのものは賞味期限切れです)。

下の「高速化」のところも,あとでお読みください。

library(RgoogleMaps)
library(RColorBrewer)
# source("GetMap.R") # 古いRgoogleMapsの場合。今は不要
fukushima = read.csv("fukushima.csv", as.is=TRUE)
t = as.POSIXct(fukushima$datetime)
o = order(t)
cols=c("white",brewer.pal(9, "YlOrRd"),rep("black",100))
FukushimaMap = GetMap(c(37.38,140.2), destfile="fukushima.png",
                      zoom=9, sensor="false", hl="ja", maptype="terrain")
# ここでURLが表示されるが無視。まだ地図は表示されない
tmp = PlotOnStaticMap(FukushimaMap, lat=fukushima$lat[o],
                      lon=fukushima$lon[o], pch=16,
                      col=cols[floor(fukushima$radiation[o]*2)+1])
# この時点で地図が表示されるはず
tmp = PlotOnStaticMap(FukushimaMap, lat=37.422778,
                      lon=141.032339, pch="×", cex=2,
                      add=TRUE) # 福島第一原発
legend(-300, -290, xjust=0, yjust=0,
       legend=paste("≧",format((10:0)/2,digits=2,nsmall=1),"μSv/h",sep=""),
       fill=cols[11:1], bg="white")
福島県の放射線地図

現在の GetMap()maptype="roadmap" がデフォルトです。デフォルトのほうが地図らしくなります。

別のサンプルも載せておきます。

福島市の放射線地図

描くものは点に限りません。文字列を描くには PlotOnStaticMap(..., FUN=text, labels=..., col=..., cex=...) のようにします。labels に文字列のベクトルを与え,cex に文字サイズ(標準は1)を与えます。

data.tableを使って書き直す

上のコードで

fukushima = read.csv("fukushima.csv", as.is=TRUE)
t = as.POSIXct(fukushima$datetime)

のところが非常に時間がかかっています。data.table のところにも書きましたが,次のようにするとずっと速くなります。

library(data.table)
library(fasttime)

fukushima = fread("fukushima.csv")
t = fastPOSIXct(fukushima$datetime) - 9*3600

ちなみに,元データの日付が YYYY/MM/DD だったり YYYY/M/D だったりするので,文字列のソートでは正しい順番に並びません。

どうせ data.table を使うなら,もうちょっとそれらしくしてみましょう。

fukushima = fread("fukushima.csv")
fukushima$datetime = fastPOSIXct(fukushima$datetime) - 9*3600
setkey(fukushima, datetime)

これで fukushima という data.table は datetime 欄をキーに設定しました。こうすると,並び順が datetime の昇順になりますので,o = order(t) を求める必要がなく,その後の部分でも [o] がすべて不要になります。

日本測地系と世界測地系

国の発表資料は日本測地系の場合がありますが,Google Mapsは世界測地系です。東京付近なら東経から12秒引き北緯に12秒足すと世界測地系になります。詳しくは国土地理院の日本測地系と世界測地系をご覧ください。Web上で変換するツールが,陸域は国土地理院のページ,海域は海上保安庁のページで提供されています。また,緯度 x' = x - 0.00010695 * x + 0.000017464 * y + 0.0046017,経度 y' = y - 0.000046038 * x - 0.000083043 * y + 0.01004 で近似計算できるそうです(Thanks: @nnistar さん)。

国土地理院が提供するTKY2JGDはWeb上・Windows上のツールですが,Macで使えるものとして BLconverter があります(Thanks: @Ishihara_Y さん)。GUIツールですがCUIでも使えます。

参考サイト(工事中)

RgoogleMapsより新しいものにggmapパッケージがあります。


Last modified: