京橋のバイオインフォマティシャンの日常

南国のビーチパラソルの下で、Rプログラムを打ってる日常を求めて、、Daily Life of Bioinformatician in Kyobashi of Osaka

【R言語と米国株】S&P500指数が年初6ヶ月でマイナス20%以上だった年の後半戦、そして翌年はどうなるのか?

2007年9月-と2021年12月-の2つの GSPCプロットの重ね合わせ(2022年10月8日アップデート)

2007-2009年と2021-2022年のS&P500の株価を縮尺を合わせて、プロットを重ねています。 株価がピークであった、2007年9月と2021年12月とをスタートにしています。 2022年10月8日に結果をアップデートしました。

2007年9月-と2021年12月-の2つの GSPCプロットの重ね合わせ(2022年9月19日アップデート)

2007-2009年と2021-2022年のS&P500の株価を縮尺を合わせて、プロットを重ねています。 株価がピークであった、2007年9月と2021年12月とをスタートにしています。 2022年9月19日に結果をアップデートしました。

 

S&P500指数が、徐々に下がってきています。やはり9月22日のFOMCは鬼門ですね。

まずは、3650あたりを割れるかどうかを見守りましょうか。

Rコード: 2つのGSPCプロットの重ね合わせ

このRコードでは、plot、axis、axis.Date、mtext関数を使って、 「スケールが異なる、2つの折れ線グラフを重ねる」ということをやっています。

#ライブラリ準備
library(quantmod)

##データ準備
#2007-2009年 S&P 500【^GSPC】の株価取得
NN01 <- 2007
Date <- c(paste0(NN01, "-09-01"), paste0(NN01+2, "-12-31"))  
GSPC <- quantmod::getSymbols("^GSPC", src = "yahoo", verbose = T, 
                               auto.assign=FALSE, from = Date[1], to=Date[2])

#データフレームの作成
GSPC01 <- data.frame(Date=index(GSPC),
                    Close=as.numeric(GSPC[,c(4)]), 
                    Month=format(index(GSPC), "%m"), 
                    row.names=1:nrow(GSPC))

#2022年 S&P 500【^GSPC】の株価取得
NN01 <- 2021
Date <- c(paste0(NN01, "-12-01"), paste0(NN01+2, "-12-31"))  
GSPC <- quantmod::getSymbols("^GSPC", src = "yahoo", verbose = T, 
                               auto.assign=FALSE, from = Date[1], to=Date[2])

#データフレームの作成
GSPC02 <- data.frame(Date=index(GSPC),
                    Close=as.numeric(GSPC[,c(4)]), 
                    Month=format(index(GSPC), "%m"), 
                    row.names=1:nrow(GSPC))

##プロット作成
#2007-2009年【^GSPC】
Dat <- GSPC01
M <- max(Dat$Close) + abs(max(Dat$Close) - min(Dat$Close))*0.1
m <- min(Dat$Close) - abs(max(Dat$Close) - min(Dat$Close))*0.1

#プロット
par(mai=c(0.75, 1, 0.75, 1), family = "HiraKakuPro-W3")
plot(Dat$Date, Dat$Close, 
     ylim = c(m, M), 
     xlab = "", ylab = "", 
     type = "l", col = "#E98683", 
     lwd = 2,
     xaxt = "n",
     yaxt = "n")

#軸作成
axis(2, at = seq(600, 1600, by=100))
axis.Date(1, at=seq.Date(min(Dat$Date), max(Dat$Date)+30, by="month")[c(T,F)], format="%B")
mtext("2007-2009年 ^GSPC", side = 2, line=2.5)
mtext("2007年9月-", side = 1, line=2.5)

#2021-2022年【^GSPC】
par(new = TRUE)
Dat0 <- rbind(GSPC02, GSPC01[c((nrow(GSPC02)+1):nrow(GSPC01)),])
Dat0$Date <- GSPC01$Date
Dat0$Close[c((nrow(GSPC02)+1):nrow(GSPC01))] <- NA

Dat <- GSPC02
M <- max(Dat$Close) + abs(max(Dat$Close) - min(Dat$Close))*0.1
m <- min(Dat$Close)/2

#Dateのズレ
Dat0$Date <- Dat0$Date+90

#プロット
plot(Dat0$Date, Dat0$Close, 
     ylim = c(m, M), 
     xlab = "", ylab = "", 
     type = "l", col = "#73CACE", 
     lwd = 2,
     axes = FALSE)
axis(4, at = seq(1800, 5000, by=400))
axis.Date(3, at=(seq.Date(min(Dat0$Date), max(Dat0$Date)+30, by="month"))[c(T,F)], format="%B")
mtext("2021-2022年 ^GSPC", side = 4, line=2.5)  
mtext("2021年12月-", side = 3, line=2.5)  

#サポートライン
abline(h=3650, col="#5ac2c7", lwd=2)
text(Dat0$Date[length(Dat0$Date)]-20, 3650+(max(Dat0$Close, na.rm = T) - min(Dat0$Close, na.rm = T))*0.1, "3650", col="#5ac2c7")

#いったんの下値目処
abline(h=2800, col="#5ac2c7", lwd=2)
text(Dat0$Date[1]+20, 2800+(max(Dat0$Close, na.rm = T) - min(Dat0$Close, na.rm = T))*0.1, "2800", col="#5ac2c7")

#凡例
legend("topright", legend = c("2007-2009年", "2021-2022年"), 
       col = c("#E98683", "#73CACE"), lwd=2, lty = 1, cex=0.8)

年初より続く、S&P500指数の下落はかなり稀だった件

もうすぐ2022年8月も終わろうとしています。 S&P500指数のリターンは、2022年初来6カ月で、-20%以上の大きな下落を記録しました。

6月16日に記録した、3666ポイントが上半期の最安値でした。

https://kumes.github.io/Blog/SP500/GSPC/GSPC_2022.html

年初来、上半期のみで同じくらい下落したことは、過去に5度あったようです。。

1932年(-41.402%)、1940年(-20.982%)、1962年(-22.844%)、1970年(-21.806%)、 そして、今年2022年(-21.081%)です。つまり、52年ぶりの出来事ということです。

消費者物価指数(CPI)、要するに、インフレ率が高止まりしていて、 これほどまでにお金の価値がどんどん下がっているのは、1970年代以来のことらしいです。

消費者物価指数(CPI)についての注釈

消費者が購入する各種の消費やサービスの小売価格の変動を調査・算出した経済指標です。 ある時点を基準に、同等のものを購入した場合に費用がどのように変動したかを指数値で表したもので、物価そのものの変動を測定することを目的としています。 季節的な影響で価格が変動しやすい生鮮食品を除いた「コア指数」が特に注目されています。 世界各国で発表され、各国(地域)のインフレ動向を測る重要な経済指標として、生産者物価指数(PPI)とともにマーケットでも注目されています。

www.daiwa.jp

1970-1971年のS&P500指数のチャート

https://kumes.github.io/Blog/SP500/GSPC/GSPC_1970.html

1970年から1971年のチャートを眺めてみると、どうも1970年の年末ごろには、年初来±0%近くまでで戻っています。 このストーリーは、まだありえるのでしょうか?

それとも、まだまだ金利が上がっている最中で、これは無いのでしょうかね。

2008年のリーマン・ショック時に似ているチャートという噂もチラホラと。。。

https://kumes.github.io/Blog/SP500/GSPC/GSPC_2008.html

2008年のリーマン・ショックで暴落が起こった時のチャートに似てきた、、、そういった巷の噂がチラホラと聴こえてきています。

2008年からのチャートと、2022年のチャートを重ねてみると、どことなく、2022年のチャートはトレースしているように見えだしてきました。

2008年からのチャートと、2022年のチャートを重ねてみると

2008年8月以降の下げで、初めのボトムは11月あたりでした。 そして、2回目のボトムは翌年の5月あたりでした。 人気ユーチューバーが言ってることと同じ?、、、そんなことを勝手に思い巡らせている今日この頃です。

補足

GSPC

^GSPCが、S&P 500のテッィカーコードです。

finance.yahoo.co.jp

Digital Color Meter.app をターミナルから開く

Digital Color Meter.appは、色のRGB値を抽出してくれるツールです。

ターミナルを起動します。以下をコピペして実行します。

open /System/Applications/Utilities/Digital\ Color\ Meter.app 

ggplotっぽい配色

#2色: シャープな赤色, 大人しい青緑色
col2 <- c( "#E98683", "#73CACE")

#3色
col3 <- c( "#e87d72", "#709bf8", "#53b74c")

#5色
col5 <- c( "#e87d72", "#a3a533", "#56bc82", "#50adf0", "#d972ec")

Rコード

パッケージの準備

quantmodとdygraphsのパッケージをセットアップします。あとは、好みで。

library(quantmod)
library(dygraphs)

S&P500 index: 1928から2022年までのデータを取得する

#S&P500 index 【^GSPC】の株価取得
GSPC.all <- quantmod::getSymbols("^GSPC", src = "yahoo", verbose = T, 
                                 from = "1900-01-01", to="2023-01-01",
                                 auto.assign=FALSE)

#表示
head(GSPC.all)
#           GSPC.Open GSPC.High GSPC.Low GSPC.Close GSPC.Volume GSPC.Adjusted
#1927-12-30     17.66     17.66    17.66      17.66           0         17.66
#1928-01-03     17.76     17.76    17.76      17.76           0         17.76
#1928-01-04     17.72     17.72    17.72      17.72           0         17.72
#1928-01-05     17.55     17.55    17.55      17.55           0         17.55
#1928-01-06     17.66     17.66    17.66      17.66           0         17.66
#1928-01-09     17.50     17.50    17.50      17.50           0         17.50

#対象年
ALLYear <- 1928:2022
#結果の格納変数
Performance <- data.frame(Year=ALLYear, Performance=NA)

#実行
for(n in seq_len(length(ALLYear))){
#検索する年月日の定義
NN00 <- ALLYear[n]
Date <- c(paste0(NN00, "-01-01"), paste0(NN00, "-07-01"))  

#S&P 500【^GSPC】の株価取得
GSPC <- quantmod::getSymbols("^GSPC", src = "yahoo", verbose = T, 
                            auto.assign=FALSE, from = Date[1], to=Date[2])

#データフレームの作成
Dat00 <- data.frame(Date=index(GSPC), 
                    Close=as.numeric(GSPC[,c(4)]), 
                    Month=format(index(GSPC), "%m"), 
                    row.names=1:nrow(GSPC))

#7月データを外す
Dat01 <- Dat00[Dat00$Month != "07",]

#6ヶ月の変動を計算
Dat01$ClosePercent <- Dat01$Close/Dat01$Close[1]*100
Performance_6month <- round(Dat01$ClosePercent[nrow(Dat01)] - 100, 3)
Performance$Performance[n] <- Performance_6month
}

年初6ヶ月間で-20%変動の年を見つける

head(Performance)
#  Year Performance
#1 1928       7.770
#2 1929      10.480
#3 1930      -3.399
#4 1931      -6.435
#5 1932     -41.402
#6 1933      59.736

#-20%未満を抽出する
Performance00 <- Performance[Performance$Performance < -20,]

Performance00
#   Year Performance
#5  1932     -41.402
#13 1940     -20.982
#35 1962     -22.844
#43 1970     -21.806
#95 2022     -21.081

1932-1933年のS&P 500 index

https://kumes.github.io/Blog/SP500/GSPC/GSPC_1932.html

#1932-1933年
NN01 <- Performance00$Year[1]
Date <- c(paste0(NN01, "-01-01"), paste0(NN01+1, "-12-31"))  

#S&P 500【^GSPC】の株価取得
GSPC <- quantmod::getSymbols("^GSPC", src = "yahoo", verbose = T, 
                            auto.assign=FALSE, from = Date[1], to=Date[2])

#全体プロット
GSPC[,4] %>%
  dygraphs::dygraph(main = paste0("S&P 500 index: ", NN01, "-", NN01+1)) %>% 
  dygraphs::dySeries("GSPC.Close", label = "GSPC") %>%
  dygraphs::dyRangeSelector(height = 50)  %>%
  dygraphs::dyAxis(name="x", axisLabelFontSize = 12, axisLabelWidth = 70)  %>%
  htmltools::save_html(file=paste0("GSPC_", NN01, ".html"))

1940-1941年のS&P 500 index

https://kumes.github.io/Blog/SP500/GSPC/GSPC_1940.html

#1940-1941年
NN01 <- Performance00$Year[2]
Date <- c(paste0(NN01, "-01-01"), paste0(NN01+1, "-12-31"))  

#S&P 500【^GSPC】の株価取得
GSPC <- quantmod::getSymbols("^GSPC", src = "yahoo", verbose = T, 
                            auto.assign=FALSE, from = Date[1], to=Date[2])

#全体プロット
GSPC[,4] %>%
  dygraphs::dygraph(main = paste0("S&P 500 index: ", NN01, "-", NN01+1)) %>% 
  dygraphs::dySeries("GSPC.Close", label = "GSPC") %>%
  dygraphs::dyRangeSelector(height = 50)  %>%
  dygraphs::dyAxis(name="x", axisLabelFontSize = 12, axisLabelWidth = 70)  %>%
  htmltools::save_html(file=paste0("GSPC_", NN01, ".html"))

1962-1963年のS&P 500 index

https://kumes.github.io/Blog/SP500/GSPC/GSPC_1962.html

#1962-1963年
NN01 <- Performance00$Year[3]
Date <- c(paste0(NN01, "-01-01"), paste0(NN01+1, "-12-31"))  

#S&P 500【^GSPC】の株価取得
GSPC <- quantmod::getSymbols("^GSPC", src = "yahoo", verbose = T, 
                            auto.assign=FALSE, from = Date[1], to=Date[2])

#全体プロット
GSPC[,4] %>%
  dygraphs::dygraph(main = paste0("S&P 500 index: ", NN01, "-", NN01+1)) %>% 
  dygraphs::dySeries("GSPC.Close", label = "GSPC") %>%
  dygraphs::dyRangeSelector(height = 50)  %>%
  dygraphs::dyAxis(name="x", axisLabelFontSize = 12, axisLabelWidth = 70)  %>%
  htmltools::save_html(file=paste0("GSPC_", NN01, ".html"))

1970-1971年のS&P 500 index

https://kumes.github.io/Blog/SP500/GSPC/GSPC_1970.html

#1970-1971年
NN01 <- Performance00$Year[4]
Date <- c(paste0(NN01, "-01-01"), paste0(NN01+1, "-12-31"))  

#S&P 500【^GSPC】の株価取得
GSPC <- quantmod::getSymbols("^GSPC", src = "yahoo", verbose = T, 
                            auto.assign=FALSE, from = Date[1], to=Date[2])

#全体プロット
GSPC[,4] %>%
  dygraphs::dygraph(main = paste0("S&P 500 index: ", NN01, "-", NN01+1)) %>% 
  dygraphs::dySeries("GSPC.Close", label = "GSPC") %>%
  dygraphs::dyRangeSelector(height = 50)  %>%
  dygraphs::dyAxis(name="x", axisLabelFontSize = 12, axisLabelWidth = 70)  %>%
  htmltools::save_html(file=paste0("GSPC_", NN01, ".html"))

2022年のS&P 500 index

https://kumes.github.io/Blog/SP500/GSPC/GSPC_2022.html

#2022年
NN01 <- Performance00$Year[5]
Date <- c(paste0(NN01, "-01-01"), paste0(NN01+1, "-12-31"))  

#S&P 500【^GSPC】の株価取得
GSPC <- quantmod::getSymbols("^GSPC", src = "yahoo", verbose = T, 
                            auto.assign=FALSE, from = Date[1], to=Date[2])

#全体プロット
GSPC[,4] %>%
  dygraphs::dygraph(main = paste0("S&P 500 index: ", NN01)) %>% 
  dygraphs::dySeries("GSPC.Close", label = "GSPC") %>%
  dygraphs::dyRangeSelector(height = 50)  %>%
  dygraphs::dyAxis(name="x", axisLabelFontSize = 12, axisLabelWidth = 70)  %>%
  htmltools::save_html(file=paste0("GSPC_", NN01, ".html"))

2008-2009年のS&P 500 index

https://kumes.github.io/Blog/SP500/GSPC/GSPC_2008.html

#2008-2009年
NN01 <- 2008
Date <- c(paste0(NN01, "-01-01"), paste0(NN01+1, "-12-31"))  

#S&P 500【^GSPC】の株価取得
GSPC <- quantmod::getSymbols("^GSPC", src = "yahoo", verbose = T, 
                            auto.assign=FALSE, from = Date[1], to=Date[2])

#全体プロット
GSPC[,4] %>%
  dygraphs::dygraph(main = paste0("S&P 500 index: ", NN01, "-", NN01+1)) %>% 
  dygraphs::dySeries("GSPC.Close", label = "GSPC") %>%
  dygraphs::dyRangeSelector(height = 50)  %>%
  dygraphs::dyAxis(name="x", axisLabelFontSize = 12, axisLabelWidth = 70)  %>%
  htmltools::save_html(file=paste0("GSPC_", NN01, ".html"))

2008-2009年と2022年の2つのプロットの重ね合わせ

このRコードでは、plot、axis、axis.Date、mtext関数を使って、 「スケールが異なる、2つの折れ線グラフを重ねる」ということをやっています。

##データ準備
#2008-2009年 S&P 500【^GSPC】の株価取得
NN01 <- 2008
Date <- c(paste0(NN01, "-01-01"), paste0(NN01+1, "-12-31"))  
GSPC <- quantmod::getSymbols("^GSPC", src = "yahoo", verbose = T, 
                               auto.assign=FALSE, from = Date[1], to=Date[2])
#データフレームの作成
GSPC01 <- data.frame(Date=index(GSPC),
                    Close=as.numeric(GSPC[,c(4)]), 
                    Month=format(index(GSPC), "%m"), 
                    row.names=1:nrow(GSPC))

#2022年 S&P 500【^GSPC】の株価取得
NN01 <- 2022
Date <- c(paste0(NN01, "-01-01"), paste0(NN01+2, "-12-31"))  
GSPC <- quantmod::getSymbols("^GSPC", src = "yahoo", verbose = T, 
                               auto.assign=FALSE, from = Date[1], to=Date[2])

#データフレームの作成
GSPC02 <- data.frame(Date=index(GSPC),
                    Close=as.numeric(GSPC[,c(4)]), 
                    Month=format(index(GSPC), "%m"), 
                    row.names=1:nrow(GSPC))

##プロット作成
#2008-2009年【^GSPC】
Dat <- GSPC01
M <- max(Dat$Close) + abs(max(Dat$Close) - min(Dat$Close))*0.1
m <- min(Dat$Close) - abs(max(Dat$Close) - min(Dat$Close))*0.1

#プロット
par(mai=c(0.75, 1, 0.5, 1), family = "HiraKakuPro-W3")
plot(Dat$Date, Dat$Close, 
     ylim = c(m, M), 
     xlab = "", ylab = "", 
     type = "l", col = "#E98683", 
     lwd = 2,
     axes = FALSE)
axis(2, at = seq(600, 1600, by=100))
mtext("2008-2009年 ^GSPC", side = 2, line=2.5)  
axis.Date(1, at=seq.Date(min(Dat$Date), max(Dat$Date)+30, by="month")[c(T,F)], format="%B")

#2022年【^GSPC】
par(new = TRUE)
Dat0 <- rbind(GSPC02, GSPC01[c((nrow(GSPC02)+1):nrow(GSPC01)),])
Dat0$Date <- GSPC01$Date
Dat0$Close[c((nrow(GSPC02)+1):nrow(GSPC01))] <- NA

Dat <- GSPC02
M <- max(Dat$Close) + abs(max(Dat$Close) - min(Dat$Close))*0.1
m <- min(Dat$Close)/2

#プロット
plot(Dat0$Date, Dat0$Close, 
     ylim = c(m, M), 
     xlab = "", ylab = "", 
     type = "l", col = "#73CACE", 
     lwd = 2,
     axes = FALSE)
axis(4, at = seq(1800, 5000, by=400))
mtext("2022年 ^GSPC", side = 4, line=2.5)  

#凡例
legend("topright", legend = c("2008-2009年", "2022年"), 
       col = c("#E98683", "#73CACE"), lwd=2, lty = 1, cex=0.8)

参考資料

skume.net

finance.yahoo.co.jp

multivariate-statistics.com

qiita.com

hira-labo.com

stats.biopapyrus.jp

【R言語と株価プロット】S&P500インデックス (SPY: SPDR S&P500 ETF) と仮想通貨(ビットコイン/BTC、イーサリアム/ETH) の季節性変動(Seasonality)を考えてみた件

(この記事は、2021年3月に執筆した記事を再校正しています。)

はじめに

2021年3月のアメリカ株式市場は、初旬からテック株、小型株がかなりの暴落でした。まさか、ここまで下がるとは、、、*1

こうなると、常識的な相場のSeasonality(季節性)をちゃんと理解して、相場の乗り降りのタイミングを考えないといけないな、身を持って強く思えてきました(泣泣泣)。

この記事では、S&P500インデックスに連動する S&P500 ETF(ティッカーシンボル: SPY)に加えて、仮想通貨のビットコイン(BTC)やイーサリアム(ETH)の価格情報を取得して、Seasonality(季節変動)をグラフ化してみることにしました。 これらの実行には、無料の解析環境である、R言語を用いています。

また、このアイディアをもとに、 月ごとの上昇・下落を色分けしてプロット作成できる、seasonalityPlotというRパッケージも作成しました。

seasonalityPlot package (version 0.99.4) | seasonalityPlot

CRAN - Package seasonalityPlot

また、記事の後半部分に、Rでの実行コードと補足事項を記載しています。

seasonalityPlotパッケージを用いた季節変動性のプロット

関連パッケージのインストール

株価などの季節変動を可視化する、seasonalityPlotパッケージを用います。

インストール方法は選べる2タイプです。

#CRANパッケージのインストール(旧バージョン)
install.packages("seasonalityPlot", repos="http://cran.r-project.org")
library(seasonalityPlot)

##GitHubパッケージのインストール(新バージョン、開発版)
install.packages("devtools", repos="http://cran.r-project.org")
devtools::install_github("kumeS/seasonalityPlot")
library(seasonalityPlot)

seasonalityPlotパッケージを使うのに、その他のパッケージは、特段必要ないです。

SPYチャートの季節変動性プロット

まずは、SPYの季節変動(2000-2021年)をプロットします。

seasonPlot(Symbols="SPY",
           StartYear = 2000,
           EndYear = 2021)

続いて、SPYの季節変動(2010-2021年)をプロットします。

seasonPlot(Symbols="SPY",
           StartYear = 2010,
           EndYear = 2021)

ビットコイン・チャートの季節変動性プロット

ビットコイン(BTC-USD)の季節変動(2014-2021年)をプロットします。

seasonPlot(Symbols="BTC-USD",
           StartYear = 2014,
           EndYear = 2021)

イーサリアム・チャートの季節変動性プロット

イーサリアム(ETH-USD)の季節変動(2015-2021年)をプロットします。

seasonPlot(Symbols="ETH-USD",
           StartYear = 2015,
           EndYear = 2021)

R言語/quantmodを用いた、SPYチャートのプロット(過去バージョン)

1994年から2021年までのSPYチャート

まずは、1994年から2021年までのSPYの株価チャートを示します。圧巻の右上がり、、2021年3月の下がりもほぼノイズの域です、、

2021年1-3月のSPYチャート

こっちが、2021年1-3月のSPYの株価チャートです。

月中旬でピークで、月終わりに下がってくるパターンではあるが、あれ、上がってるけど、、、大型セクターは結構安定なことをよく物語っています。

SPYチャートのいろんなバリエーション

次に、3パターンの図を作成してみました*2

2000-2020年平均のパターンだと、2月中旬がピークで、3月3週目あたりが底となるみたいです。

2010-2020年平均のパターンも上記と形が似てるけど、3月初めからの下りが激しく、3月最終週あたりから上がりだすようです。

最後に、大統領選翌年平均のパターンだけど、最悪想定ではあるが、4月中旬まで底が続くパターンがあるようです。

まとめ

これくらいのパターンは想定して、春相場を見ていかないとダメと言うことなのね。

ネットを見てると、3月は ・ 債券のショート ・ 株 => 債券へのリバランス の傾向が強いらしく、2月までに株で儲かり過ぎたのだろう。

次は、8月あたりで売りつけれるかが山場になりそうだけど。。。

補足

strftime関数などにおけるformatについて

%A: 曜日の名称(Sunday, Monday ... )
%a: 曜日の省略名(Sun, Mon ... )
%B: 月の名称(January, February ... )
%b: 月の省略名(Jan, Feb ... )
%j: 年中の通算日(001-366)
%U: 週を表す数。最初の日曜日が第1週の始まり(00-53)
%u: 月曜日を1とした、曜日の数値表現 (1..7)
%V: ISO 8601形式の暦週 (01..53)
%v: VMS形式の日付 (%e-%^b-%4Y)
%W: 週を表す数。最初の月曜日が第1週の始まり(00-53)
%w: 曜日を表す数。日曜日が0(0-6)

plot図XY軸の遊びをなくする引数

このオプションを用いることで、プロット領域の両端をxlim、およびylimで指定した範囲のプロットを生成します。

xaxs="i"
yaxs="i"

text()のちょっと高級版: 背景色とかが設定できる。

plotrix::boxed.labels(x,y=NULL,labels,
  bg=ifelse(match(par("bg"),"transparent",0),"white",par("bg")),
  border=TRUE,xpad=1.2,ypad=1.2,srt=0,cex=1,adj=0.5,xlog=FALSE,ylog=FALSE)

月の英語略語

month.abb
# [1] "Jan" "Feb" "Mar" "Apr" "May" "Jun" "Jul" "Aug" "Sep"
#[10] "Oct" "Nov" "Dec"

Rでの実行コード

#パッケージ・ロード
library(quantmod)
library(magrittr)
library(dygraphs)
library(plotrix)
library(htmltools)

#SPYの株価取得
NN01 <- 1994; NN02 <- 2021
Date <- c(paste0(NN01, "-01-01"), paste0(NN02, "-12-31"))
SPY <- quantmod::getSymbols("SPY", src = "yahoo", verbose = T, auto.assign=FALSE, from = Date[1], to=Date[2])

head(SPY)
#           SPY.Open SPY.High  SPY.Low SPY.Close SPY.Volume SPY.Adjusted
#1994-01-03 46.59375 46.65625 46.40625  46.46875     960900     28.06298
#1994-01-04 46.53125 46.65625 46.46875  46.65625     164300     28.17622
#1994-01-05 46.71875 46.78125 46.53125  46.75000     710900     28.23285
#1994-01-06 46.81250 46.84375 46.68750  46.75000     201000     28.23285
#1994-01-07 46.84375 47.06250 46.71875  47.03125     775500     28.40267
#1994-01-10 47.09375 47.59375 46.96875  47.59375     593700     28.74237

str(SPY)
#An ‘xts’ object on 1994-01-03/2021-03-25 containing:
#  Data: num [1:6856, 1:6] 46.6 46.5 46.7 46.8 46.8 ...
# - attr(*, "dimnames")=List of 2
#  ..$ : NULL
#  ..$ : chr [1:6] "SPY.Open" "SPY.High" "SPY.Low" "SPY.Close" ...
#  Indexed by objects of class: [Date] TZ: UTC
#  xts Attributes:  
#List of 2
# $ src    : chr "yahoo"
# $ updated: POSIXct[1:1], format: "2021-03-26 00:43:34"

#全体プロット
SPY[,4] %>%
  dygraphs::dygraph(main = "SPDR S&P 500 (SPY) Stock Price 26-years") %>% 
  dygraphs::dySeries("SPY.Close", label = "SPY") %>%
  dygraphs::dyRangeSelector(height = 40)  %>%
  htmltools::save_html(file="SPY_00.html")

#xtsオブジェクトをデータフレームに変換する
#データは終値だけにする
SPY00 <- data.frame(Date=index(SPY), SPY[,c(4)], row.names=1:nrow(SPY))
head(SPY00)
#        Date SPY.Close
#1 1994-01-03  46.46875
#2 1994-01-04  46.65625
#3 1994-01-05  46.75000
#4 1994-01-06  46.75000
#5 1994-01-07  47.03125
#6 1994-01-10  47.59375

#年、週目などの情報取得
SPY00$Year <- substring(SPY00$Date, 1, 4)
SPY00$Week <- strftime(SPY00$Date, format = "%V")
SPY00$Year.Week <- paste0(SPY00$Year, SPY00$Week)
SPY00$Day <- sub("-", "", substring(SPY00$Date, 6, 10))
head(SPY00)
#        Date SPY.Close Year Week Year.Week  Day
#1 1994-01-03  46.46875 1994   01    199401 0103
#2 1994-01-04  46.65625 1994   01    199401 0104
#3 1994-01-05  46.75000 1994   01    199401 0105
#4 1994-01-06  46.75000 1994   01    199401 0106
#5 1994-01-07  47.03125 1994   01    199401 0107
#6 1994-01-10  47.59375 1994   02    199402 0110
tail(SPY00)
#           Date SPY.Close Year Week Year.Week  Day
#6851 2021-03-18    391.48 2021   11    202111 0318
#6852 2021-03-19    389.48 2021   11    202111 0319
#6853 2021-03-22    392.59 2021   12    202112 0322
#6854 2021-03-23    389.50 2021   12    202112 0323
#6855 2021-03-24    387.52 2021   12    202112 0324
#6856 2021-03-25    387.90 2021   12    202112 0325

#いったん、1年372日分作成
Days <- c()
for(n in 1:12){
Days <- c(Days, paste0(formatC(n, width=2, flag = "0"), formatC(seq(1, 31, 1), width=2, flag = "0")))
}
SPY01 <- data.frame(Date=Days)
Years <- paste0("Y", seq(NN01, NN02, 1))
SPY01[,Years] <- NA
head(SPY01)
#  Date Y1994 Y1995 Y1996 Y1997 Y1998 Y1999 Y2000 Y2001 Y2002 Y2003 Y2004 Y2005 Y2006 Y2007 Y2008 Y2009 Y2010
#1 0101    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA
#2 0102    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA
#3 0103    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA
#  Y2011 Y2012 Y2013 Y2014 Y2015 Y2016 Y2017 Y2018 Y2019 Y2020 Y2021
#1    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA
#2    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA
#3    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA    NA

#日付、年ごとにデータ整理
for(m in 1:length(Years)){
a <- SPY00[SPY00$Year == sub("Y", "", Years[m]),]
SPY01[c(SPY01$Date %in% a$Day),m+1] <- a$SPY.Close
}

#ここは、特に実行しなくて良い
#欠損数のカウント
#SPY01$MissingNum <- sapply(data.frame(t(SPY01[,-1])), function(x) sum(is.na(x)))
#株価データが無い日を消す。
#table(SPY01$MissingNum == length(Years))
#SPY01 <- SPY01[SPY01$MissingNum != length(Years),]
#head(SPY01)

#それぞれの年初日を「0%」とする
SPY02 <- SPY01
for(m in 1:length(Years)){
b <- SPY02[,m+1]
SPY02[,m+1] <- SPY02[,m+1]/b[!is.na(b)][1]*100 - 100
}
head(SPY02)
#  Date     Y1994     Y1995      Y1996     Y1997      Y1998    Y1999
#1 0101        NA        NA         NA        NA         NA       NA
#2 0102        NA        NA  0.0000000 0.0000000  0.0000000       NA
#3 0103 0.0000000 0.0000000  0.2765904 1.4352047         NA       NA
#4 0104 0.4034970 0.4778157 -0.6789037        NA         NA 0.000000
#5 0105 0.6052455 0.4778157 -0.8800603        NA  0.2242152 1.143002
#6 0106 0.6052455 0.5802048         NA 0.5487547 -1.3773222 3.581407

#株価データが無い日を削除して、株価データを繋げる
SPY03 <- data.frame(Date=1:nrow(SPY01))
SPY03[,Years] <- NA
for(m in 1:length(Years)){
d <- SPY02[,m+1][!is.na(SPY02[,m+1])]
SPY03[,m+1] <- c(d, rep(NA, nrow(SPY01)-length(d)))
}

#株価データが無い(欠損数2以上)行を削除する
SPY03$MissingNum <- sapply(data.frame(t(SPY03[,-1])), function(x) sum(is.na(x)))
SPY03 <- SPY03[c(SPY03$MissingNum < 2),]

SPY04 <- SPY03
colnames(SPY04)
# [1] "Date"       "Y1994"      "Y1995"      "Y1996"      "Y1997"      "Y1998"      "Y1999"      "Y2000"     
# [9] "Y2001"      "Y2002"      "Y2003"      "Y2004"      "Y2005"      "Y2006"      "Y2007"      "Y2008"     
#[17] "Y2009"      "Y2010"      "Y2011"      "Y2012"      "Y2013"      "Y2014"      "Y2015"      "Y2016"     
#[25] "Y2017"      "Y2018"      "Y2019"      "Y2020"      "Y2021"      "MissingNum"
dim(SPY04)
#[1] 248  30

#図 1. 2021年 SPY のSeasonality
Main <- "図 1. 2021年 SPY のSeasonality"
Y <- c(29)
SPY04$Mean <- NA
SPY04$Mean <- apply(data.frame(SPY04[,c(Y)]), 1, function(x) mean(x[!is.na(x)]))

COL <- "blue1"
St <- SPY04$Mean; St00 <- St[!is.na(St)]
par(family="HiraKakuProN-W3", lwd=1, xpd=F, cex=1, mgp=c(0.5, 1, 0), mai=c(0.5, 0.75, 0.5, 0.5))
plot(St, type="n", axes = F, xlab="© 2021 京橋のバイオインフォマティシャンの日常 by skume", 
     ylab="", cex.lab=0.75, xlim=c(0,nrow(SPY04)), 
     ylim=c(min(St00) - (max(St00)- min(St00))*0.1, max(St00) + (max(St00)- min(St00))*0.15), 
     xaxs="i", yaxs="i", main=Main)
axis(side=2, labels=paste0(seq(-30, 30, by=1), "%"), at=seq(-30, 30, by=1), las=2)
abline(v=seq(0, nrow(SPY04), length.out = 13), col="grey", lty=3, lwd=0.5)
abline(h=seq(-30, 30, by=1), col="black", lty=1, lwd=0.3)
lines(St, col=COL, lwd=1.2)
plotrix::boxed.labels(seq(0, nrow(SPY04), length.out = 13)[1:12]+10,  min(St00) - (max(St00)- min(St00))*0.05, 
                      month.abb, cex=0.7, bg = "grey20", xpad = 1.2, ypad = 1.2)
legend("topleft", legend="SPDR S&P500 ETF", col=COL, lwd=1, cex=0.7)
quartz.save(file = paste0("./SPY_01.png"), type = "png", dpi = 300); dev.off()

#図 2. 2000-2020年 SPY のSeasonality
Main <- "図 2. 2000-2020年 SPY のSeasonality"
Y <- c(8:28)
#head(SPY04[,c(Y)])
SPY04$Mean <- NA
SPY04$Mean <- apply(data.frame(SPY04[,c(Y)]), 1, function(x) mean(x[!is.na(x)]))

COL <- "blue1"
St <- SPY04$Mean; St00 <- St[!is.na(St)]
par(family="HiraKakuProN-W3", lwd=1, xpd=F, cex=1, mgp=c(0.5, 1, 0), mai=c(0.5, 0.75, 0.5, 0.5))
plot(St, type="n", axes = F, xlab="© 2021 京橋のバイオインフォマティシャンの日常 by skume", 
     ylab="", cex.lab=0.75, xlim=c(0,nrow(SPY04)), 
     ylim=c(min(St00) - (max(St00)- min(St00))*0.1, max(St00) + (max(St00)- min(St00))*0.15), 
     xaxs="i", yaxs="i", main=Main)
axis(side=2, labels=paste0(seq(-30, 30, by=1), "%"), at=seq(-30, 30, by=1), las=2)
abline(v=seq(0, nrow(SPY04), length.out = 13), col="grey", lty=3, lwd=0.5)
abline(h=seq(-30, 30, by=1), col="black", lty=1, lwd=0.3)
lines(St, col=COL, lwd=1.2)
plotrix::boxed.labels(seq(0, nrow(SPY04), length.out = 13)[1:12]+10,  min(St00) - (max(St00)- min(St00))*0.05, 
                      month.abb, cex=0.7, bg = "grey20", xpad = 1.2, ypad = 1.2)
legend("topleft", legend="SPDR S&P500 ETF", col=COL, lwd=1, cex=0.7)
quartz.save(file = paste0("./SPY_02.png"), type = "png", dpi = 300); dev.off()

#図 3. 2010-2020年 SPY のSeasonality
Main <- "図 3. 2010-2020年 SPY のSeasonality"
Y <- c(18:28)
#head(SPY04[,c(Y)])
SPY04$Mean <- NA
SPY04$Mean <- apply(data.frame(SPY04[,c(Y)]), 1, function(x) mean(x[!is.na(x)]))

COL <- "red1"
St <- SPY04$Mean; St00 <- St[!is.na(St)]
par(family="HiraKakuProN-W3", lwd=1, xpd=F, cex=1, mgp=c(0.5, 1, 0), mai=c(0.5, 0.75, 0.5, 0.5))
plot(St, type="n", axes = F, xlab="© 2021 京橋のバイオインフォマティシャンの日常 by skume", 
     ylab="", cex.lab=0.75, xlim=c(0,nrow(SPY04)), 
     ylim=c(min(St00) - (max(St00)- min(St00))*0.1, max(St00) + (max(St00)- min(St00))*0.15), 
     xaxs="i", yaxs="i", main=Main)
axis(side=2, labels=paste0(seq(-30, 30, by=1), "%"), at=seq(-30, 30, by=1), las=2)
abline(v=seq(0, nrow(SPY04), length.out = 13), col="grey", lty=3, lwd=0.5)
abline(h=seq(-30, 30, by=1), col="black", lty=1, lwd=0.3)
lines(St, col=COL, lwd=1.2)
plotrix::boxed.labels(seq(0, nrow(SPY04), length.out = 13)[1:12]+10,  min(St00) - (max(St00)- min(St00))*0.05, 
                      month.abb, cex=0.7, bg = "grey20", xpad = 1.2, ypad = 1.2)
legend("topleft", legend="SPDR S&P500 ETF", col=COL, lwd=1, cex=0.7)
quartz.save(file = paste0("./SPY_03.png"), type = "png", dpi = 300); dev.off()

#図 4. 大統領選翌年 SPYのSeasonality
#'97 '01 '05 '09 '13 '17
Main <- "図 4. 大統領選翌年 SPYのSeasonality"
Y <- c(5, 9, 13, 17, 21, 25)
head(SPY04[,c(Y)])
#      Y1997     Y2001     Y2005      Y2009       Y2013     Y2017
#1 0.0000000 0.0000000  0.000000  0.0000000  0.00000000 0.0000000
#2 1.4352047 4.8034934 -1.221946 -0.1183315 -0.22593592 0.5949196
#3 0.5487547 3.6754003 -1.903575  0.5486252  0.21223949 0.5150013

SPY04$Mean <- NA
SPY04$Mean <- apply(data.frame(SPY04[,c(Y)]), 1, function(x) mean(x[!is.na(x)]))
colnames(SPY04)

COL <- "green1"
St <- SPY04$Mean; St00 <- St[!is.na(St)]
par(family="HiraKakuProN-W3", lwd=1, xpd=F, cex=1, mgp=c(0.5, 1, 0), mai=c(0.5, 0.75, 0.5, 0.5))
plot(St, type="n", axes = F, xlab="© 2021 京橋のバイオインフォマティシャンの日常 by skume", 
     ylab="", cex.lab=0.75, xlim=c(0,nrow(SPY04)), 
     ylim=c(min(St00) - (max(St00)- min(St00))*0.1, max(St00) + (max(St00)- min(St00))*0.15), 
     xaxs="i", yaxs="i", main=Main)
axis(side=2, labels=paste0(seq(-30, 30, by=1), "%"), at=seq(-30, 30, by=1), las=2)
abline(v=seq(0, nrow(SPY04), length.out = 13), col="grey", lty=3, lwd=0.5)
abline(h=seq(-30, 30, by=1), col="black", lty=1, lwd=0.3)
lines(St, col=COL, lwd=1.2)
plotrix::boxed.labels(seq(0, nrow(SPY04), length.out = 13)[1:12]+10,  min(St00) - (max(St00)- min(St00))*0.05, 
                      month.abb, cex=0.7, bg = "grey20", xpad = 1.2, ypad = 1.2)
legend("topleft", legend="SPDR S&P500 ETF", col=COL, lwd=1, cex=0.7)
quartz.save(file = paste0("./SPY_04.png"), type = "png", dpi = 300); dev.off()

参考資料

stackoverflow.com

rstudio.github.io

docs.ruby-lang.org

stackoverflow.com

*1:じっちゃま曰く、阿波踊りして、降りなかったヤツが悪いのですがね

*2:自由に転載可です

【R言語とオープンデータ】アメリカ地質調査所が提供する世界中の地震データを地図表示してみた件 〜初心者でもできるオープンデータを使った世界地図の簡単プロット〜

はじめに

R言語を使った、オープンデータと地図可視化というテーマで記事執筆してみました。

オープンデータとしては、 アメリカ地質調査所が提供する、世界中の地震データを活用させてもらいました。 このオープンデータは、CSVデータのURL指定で簡単なスクレイピングで取得できるので楽ちんで利用できます。

オープンデータの関連記事

skume.net

アメリカ地質調査所、および同研究所が提供する世界で発生した地震データのスプレッドシート

 

アメリカ地質調査所(USGS)は、 アメリカ合衆国政府の科学研究機関の一つで、 水文科学、生物学、地質学、地理学の4つの主要な科学分野についての ナチュラル・ハザードを対象とする調査研究を幅広く行なっているらしいです。

今回、USGSが提供している地震データのスプレッドシートを利用させてもらいます。

同研究所HPの「Real-time Notifications, Feeds, and Web Services」 の「Spreadsheet Format」項目から地震データが得られます。

earthquake.usgs.gov

過去何日かと、どの程度のマグニチュード(M)範囲の地震データを取得するかで、 対象データのURLが変わってきます。

以下に、各地震データのURLを示します。

過去1日地震データのスプレッドシートURL

Past Day URL
Significant Earthquakes https://earthquake.usgs.gov/earthquakes/feed/v1.0/summary/significant_day.csv
M4.5+ Earthquakes https://earthquake.usgs.gov/earthquakes/feed/v1.0/summary/4.5_day.csv
M2.5+ Earthquakes https://earthquake.usgs.gov/earthquakes/feed/v1.0/summary/2.5_day.csv
M1.0+ Earthquakes https://earthquake.usgs.gov/earthquakes/feed/v1.0/summary/1.0_day.csv
All Earthquakes https://earthquake.usgs.gov/earthquakes/feed/v1.0/summary/all_day.csv

過去7日地震データのスプレッドシート

Past 7 Days URL
Significant Earthquakes https://earthquake.usgs.gov/earthquakes/feed/v1.0/summary/significant_week.csv
M4.5+ Earthquakes https://earthquake.usgs.gov/earthquakes/feed/v1.0/summary/4.5_week.csv
M2.5+ Earthquakes https://earthquake.usgs.gov/earthquakes/feed/v1.0/summary/2.5_week.csv
M1.0+ Earthquakes https://earthquake.usgs.gov/earthquakes/feed/v1.0/summary/1.0_week.csv
All Earthquakes https://earthquake.usgs.gov/earthquakes/feed/v1.0/summary/all_week.csv

過去30日地震データのスプレッドシート

Past 30 Days URL
Significant Earthquakes https://earthquake.usgs.gov/earthquakes/feed/v1.0/summary/significant_month.csv
M4.5+ Earthquakes https://earthquake.usgs.gov/earthquakes/feed/v1.0/summary/4.5_month.csv
M2.5+ Earthquakes https://earthquake.usgs.gov/earthquakes/feed/v1.0/summary/2.5_month.csv
M1.0+ Earthquakes https://earthquake.usgs.gov/earthquakes/feed/v1.0/summary/1.0_month.csv
All Earthquakes https://earthquake.usgs.gov/earthquakes/feed/v1.0/summary/all_month.csv

 

今回、過去30日のM4.5+の地震データを使用することにします。

スプレッドシートには、 発生日時、緯度・軽度、マグニチュードに加えて、 関連するメタデータも付随しています。

leafletパッケージを使って、地震データを地図上に表示させる

実行環境

以下のMac/RStudioの実行環境で行なっています。

#RStudioターミナル上のR環境
R version 4.1.2 (2021-11-01) -- "Bird Hippie"
Copyright (C) 2021 The R Foundation for Statistical Computing
Platform: x86_64-apple-darwin17.0 (64-bit)

関連パッケージのインストール

地図表示には、leafletパッケージを用いています。

Leafletは、インタラクティブマップ表示のための人気のJavaScriptライブラリの1つで、Leaflet for RはR環境で簡単に使えるように設計されています。

The New York Times や The Washington Post から GitHub や Flickr まで、 またOpenStreetMap や Mapbox、CartoDB などの GIS 専門サイトでも使用されています。

インストール方法はCRANとGitHubからの選べる2タイプです。

#CRANパッケージのインストール(安定版)
install.packages("leaflet")
library(leaflet)

##GitHubパッケージのインストール(開発版)
install.packages("devtools")
devtools::install_github("rstudio/leaflet")
library(leaflet)

#その他パッケージ
install.packages("htmltools")
library(htmltools)

今回は、CRANパッケージの方を使用しました。

USGSの地震スプレッドシート(csvダイル)を読み込み

上記の通り、過去30日のM4.5+地震のCSVデータを使います。

地震データの緯度・経度の位置情報をもとに、 地震の発生場所を地図上にプロットして、各メタデータもポップアップできるようにします。

要するに、 世界のどこで、どの程度の地震が発生したのかを地図上に可視化してみます。

さて、R環境で実行していきます。

#地震データのURL定義
csv_path <- "https://earthquake.usgs.gov/earthquakes/feed/v1.0/summary/4.5_month.csv"

#CSV読み込み
Earthquakes <- read.table(file=csv_path, header = T, sep = ",")

#ヘッダー表示
head(Earthquakes)

#                      time latitude longitude   depth mag magType nst gap  dmin  rms net
#1 2022-08-16T08:34:12.247Z -22.4001  170.3261  34.758 5.0      mb  37 115 2.671 1.18  us
#2 2022-08-16T07:13:09.474Z  10.2714  125.5214  78.941 4.6      mb  49 116 3.182 1.14  us
#3 2022-08-16T07:09:48.807Z -54.8013  158.5708  10.000 5.2      mb  27  93 0.377 0.67  us
#4 2022-08-16T07:02:09.836Z -21.6775 -178.5008 438.673 4.5      mb  64  91 3.113 0.65  us
#5 2022-08-16T03:46:26.016Z -31.3193 -177.8050  10.000 4.7      mb  20 151 2.053 0.97  us
#6 2022-08-16T03:44:36.298Z  13.2068  -89.5779  59.953 4.6     mwr  38 193 0.462 0.57  us
#          id                  updated                            place       type
#1 us6000ib98 2022-08-16T08:53:12.040Z southeast of the Loyalty Islands earthquake
#2 us6000ib8x 2022-08-16T07:58:07.066Z 8 km SSW of Tubajon, Philippines earthquake
#3 us6000ib8t 2022-08-16T07:52:45.069Z          Macquarie Island region earthquake
#4 us6000ib8s 2022-08-16T07:41:25.040Z                      Fiji region earthquake
#5 us6000ib81 2022-08-16T04:17:21.040Z          Kermadec Islands region earthquake
#6 us6000ib7v 2022-08-16T05:58:35.878Z                      El Salvador earthquake
#  horizontalError depthError magError magNst   status locationSource magSource
#1           11.41      5.995    0.094     38 reviewed             us        us
#2           10.96      7.416    0.066     69 reviewed             us        us
#3            7.60      1.869    0.122     22 reviewed             us        us
#4           12.56      6.868    0.057     90 reviewed             us        us
#5           11.91      1.968    0.126     19 reviewed             us        us
#6            5.96      5.740    0.071     19 reviewed             us        us

地図プロットの作成には、実際に、 2列目のlatitude(緯度)、3列目のlongitude(経度)、5列目のmag(マグニチュード)、place(場所)などのメタデータを使っていきます。

地震データの地図表示 (1) clusterOptions無しで作図する

leafletパッケージで地図プロットを描く際に、 マーカーのclusterOptionsのON/OFFで結構見た目が変わってしまいます。

そこで、それぞれの条件で、マップ表示を作成することにしました。

まずは、clusterOptionsが「OFF」の場合で実行していきます。

#地図の下準備
#head(Earthquakes)
Earthquakes <- Earthquakes[rev(order(Earthquakes$mag)),]
map <- leaflet::leaflet(Earthquakes) %>% addTiles()

#取得時間: addControl用のタイトルhtml作成
tim <- base::substr(Earthquakes[rev(order(Earthquakes$time)),][1,1], 
                    start=1, stop=16)
title <- tags$p(tags$style("p {color: black; font-size:12px}"), tags$b(tim))

##ラベル作成
Lab <- paste0("日時: ", Earthquakes$time, "<br>",
              "場所: ", Earthquakes$place, "<br>",
              "マグニチュード: ", Earthquakes$mag, "<br>",
              "深さ: ", Earthquakes$depth)

#マーカーサイズ(1)
size1 <- (2^Earthquakes$mag)*600

#カラーパレット
pal <- leaflet::colorNumeric(palette="Reds", domain=Earthquakes$mag)

#地図作成: clusterOptions無し
map1 <- map %>% 
    addProviderTiles(providers$OpenMapSurfer) %>% 
    addCircles(lng=~longitude,lat=~latitude,
                     radius=size1, color=~pal(Earthquakes$mag), weight=1, 
                     stroke = TRUE, fillOpacity = 0.4, 
                     popup = Lab,
                     label =  ~as.character(mag),
                     labelOptions = labelOptions(noHide = F, direction = 'center', textOnly = T)) %>% 
     addMiniMap(position="bottomright", width = 75, height = 75) %>% 
     addLegend(position='topright', pal=pal, values=~Earthquakes$mag, title="mag", opacity = 0.6) %>% 
    addControl(title, position = "bottomleft" )

#地図表示
map1

#地図の保存
htmlwidgets::saveWidget(map1, "./map1.html", selfcontained = F)

地震データの地図表示 (1)の出力例

図下のURLをクリックするとhtmlマップが表示され、ズームイン・ズームアウトもできます。

https://kumes.github.io/Blog/R_map_leaflet/map1.html

バヌアツあたりの拡大マップです。

地震データの地図表示 (2) clusterOptions有りで作図する

続いて、clusterOptionsが「ON」の場合で実行していきます。

#マーカーサイズ(2)
size2 <- (2^Earthquakes$mag)/2

#地図作成: clusterOptions有り
map2 <- map %>% 
    addProviderTiles(providers$OpenMapSurfer) %>% 
    addCircleMarkers(lng=~longitude,lat=~latitude,
                     radius=size2, color="#09f", weight=1,
                     clusterOptions = markerClusterOptions(),
                     popup=Lab,
                     label =  ~as.character(mag),
                     labelOptions = labelOptions(noHide = T, direction = 'center', textOnly = T)
                     ) %>% 
    addMiniMap(position="bottomright", width = 75, height = 75) %>% 
    addControl(title, position = "bottomleft" )

#地図表示
map2

#地図の保存
htmlwidgets::saveWidget(map2, "./map2.html", selfcontained = F)

地震データの地図表示 (2)の出力例

図下のURLをクリックするとhtmlマップが表示され、ズームイン・ズームアウトもできます。

https://kumes.github.io/Blog/R_map_leaflet/map2.html

日本の関東地方を拡大してみました。

地震マップ可視化の自作関数化

以下の実行で、上記と同じマップが得られます。

#ブラウザでRコードを見る
browseURL("https://gist.github.com/kumeS/47785a19a29e0f943b0a3ea695ddf8cf")

#自作関数のソース
source("https://gist.githubusercontent.com/kumeS/47785a19a29e0f943b0a3ea695ddf8cf/raw/db58c921e62863b45566758c6afed61395512f4c/Earthquakes.map.R")

#Type=1: 地図作成でclusterOptions無し
Earthquakes.map(Type=1)

#Type=2: 地図作成でclusterOptions有り
Earthquakes.map(Type=2)

#Type=3: 地図作成 old-version
Earthquakes.map(Type=3)

まとめ

今回、R環境下で地図を描く時に便利なツールとして、 leafletパッケージを紹介しました。 地震データ以外のオープンデータと組み合わせることで、 様々なシーンで活用できそうです。

また、kazutanさんの公開ページや公式ドキュメントも見ると、色々と参考になります。

参考資料

rstudio.github.io

kazutan.github.io

stackoverflow.com