R - 地震活動状況の地図表示！
Updated:
少し前から国土地理院や国土交通省提供の地図情報（Shape ファイル）から R(maptools) を用いて地図を表示してきました。
今回は、国土地理院・地球地図日本上に気象庁発表の地震情報（震源・地震規模）の情報を重ねて表示してみした。
（乱文ご容赦下さい）
0. 前提条件
- R 3.1.1 での作業を想定。
- SHAPE ファイルは国土地理院・地球地図日本を使用。（別途都道府県単位で地物融合済み）
1. 地震データ（CSV ファイル）の準備
使用する地震データを何かしらの形で準備する。
当方は、以前から取り貯めていた「気象庁防災情報 XML」の「震源・震度に関する情報」のデータを整形していかのような CSV ファイルを用意した。（2013/01/01から2014/06/30までの3,357件）
File: eq_data.csv
1
2
3
4
5
6
7
8
9
10
11
12
13
14
lat,lon,depth,magnitude
-23.1,-176.6,210,7.4
-03.8,154.0,380,7.3
39.1,142.1,50,3.6
40.5,141.7,100,3.9
27.1,127.5,90,4.0
38.9,142.0,50,3.5
40.1,141.8,10,3.1
35.7,141.1,10,4.2
39.0,142.1,50,3.7
:
===< 以降省略 >===
:
2. R スクリプト作成
以下のような R スクリプト（プログラム）を作成した。
処理の詳細は説明しない。スクリプト中のコメントで概要を理解してください。
File: generate_map_prefs.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
# 国土地理院・地球地図日本に地震情報をプロット
# ライブラリの読み込み
library(gpclib)
library(ggplot2)
library(maptools)
# gpclib ライセンス警告表示の抑止
gpclibPermit()
# Shapefile, CSV ファイルのフルパス
shp_file <- "/path/to/polbnda_jpn_pref.shp"
csv_file <- "/path/to/eq_data.csv"
# Shapefile の読み込み
shp <- readShapePoly(shp_file)
# Shapefile をデータフレーム形式に変換
df <- fortify(shp)
# CSV データ読み込み
csv <- read.csv(csv_file)
# タイトル・コピーライト文字列
title <- "地震の活動状況 [ 2013/01/01 - 2014/06/30 ]"
copy <- paste(
"（出典：国土地理院・地球地図日本 / 気象庁・防災情報XML）",
"© 2014 mk-mode.com", sep = "\n"
)
# 地図描画
g <- ggplot(df) # オブジェクト生成
g <- g + ggtitle(title) # グラフタイトル設定
g <- g + geom_polygon(
aes(long, lat, group = group),
colour = "gray20", fill = "darkolivegreen4", size = 0.1
) # 日本地図描画
g <- g + geom_point(
data = csv, shape = 1, alpha = 0.5,
aes(x = lon, y = lat, size = magnitude, colour = depth)
) # 地震情報 Plot
g <- g + scale_colour_gradient(low = "red", high="midnightblue") # 深さ色設定
g <- g + xlim(c(122, 150)) + ylim(c(23, 47)) # プロット範囲指定
g <- g + coord_equal() # メモリ刻みを等間隔設定
g <- g + theme(
plot.background = element_rect(
colour = "black", size = 0.5
)
) # プロット領域背景
g <- g + labs(
x = copy, y = "", size = "規模(M)", colour = "深さ(km)"
) # ラベル設定
g <- g + theme(legend.position = "right") # 凡例位置
g <- g + theme(
title = element_text(size = 10, colour = "black"), # タイトルのサイズ・色
axis.title.y = element_blank(), # y 軸ラベル領域非表示
axis.title.x = element_text(size = 10, colour = "gray20") # x 軸ラベルのサイズ・色
)
g <- g + guides(
colour = guide_legend(order = 1),
size = guide_legend(order = 2)
)
g <- g + theme(
panel.background = element_rect(
fill = "lightsteelblue", colour = "black",
size= 0.2 , linetype = 1
)
) # グラフ枠・背景
# 画像保存
ggsave(
file = paste("earthquakes.png"),
dpi = 100, width = 6.4, height = 6.4,
g
)
3. R スクリプト実行
R を起動し以下のように R スクリプトを指定して実行する。（表示範囲外の震源データは警告となる）
> source("earthquake.R")
コマンドラインから実行する場合は以下のように実行する。
$ R --vanilla --slave < earthquake.R
--vanilla は、R がこれまでに保存していたオブジェクトを「読み込まない」（プレーンな）状態で R を実行するオプション。
--slave は、R の出力（処理内容）を標準出力（ディスプレイ）に表示しないオプション。
4. 地図確認
成功すると “earthquakes.png” という画像ファイルが作成されるので表示して確認してみる。
「数値データを視覚化することで状況をより把握しやすくなる」ということがよく分かります。
以上。
Comments