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

南国のビーチパラソルの下で、Rプログラムを打ってる日常を求めて、、

【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)

出力例

クリックするとhtmlマップが表示され、ズームインもできます。

[f:id:skume:20220816232353p:plain]
バヌアツあたりの拡大マップです。

地震データの地図表示: (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)

出力例

クリックするとhtmlマップが表示され、ズームインもできます。

[f:id:skume:20220816232515p:plain]
日本の関東地方を拡大してみました。

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

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

#ブラウザで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さんの公開ページや公式ドキュメントも見ると、色々と参考になります。

参考資料

[https://rstudio.github.io/leaflet/:embed:cite]

[https://kazutan.github.io/JapanR2015/leaflet_d.html:embed:cite]

[https://stackoverflow.com/questions/53465072/adding-blank-legend-in-r-leaflet-package:embed:cite]

R環境で小説のテキストマイニングをやってみたら、○○○な結末になった件【その4: テキストマイニングと形態素のワードクラウド】

はじめに: 『R環境で小説のテキストマイニング』の連載シリーズ

テキストマイニングは、テキストデータから、有益な情報を取り出すデータマイニング手法の1つです。 テキストデータに対する情報解析では、自然言語処理、形態素解析、キーワード抽出、共起分析、ネットワーク可視化、機械学習、感情分析など、様々な分析手法が用いられます。 近年では、TwitterのつぶやきやSNS、地図上のクチコミなどのテキストデータも活用されています。

この記事は、夏目漱石の「坊っちゃん」が対象小説で、 テキストマイニングと形態素のワードクラウドを扱った内容となっています。

連載シリーズの目次

1. 青空文庫、対象小説の紹介
2. 「坊っちゃん」のテキストの前処理
3. 形態素解析と辞書設定
4. 形態素解析と複合語抽出 (ルールベース抽出、pytermextract)
5. テキストマイニングと形態素のワードクラウド
6. テキストマイニングとN-gramのネットワーク表現(章ごとの関係性とか)
7. 文章の感情分析(感情ポジネガ分析とか、英語感情対応表とか、ML-askとか)
8. ミクロとマクロなテキスト解析〜応用篇として内容を考え中〜


この記事では、「5. テキストマイニングと形態素のワードクラウド」に関して、

の内容をやっていきます。

まずは、実行環境

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

#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)

今回、M1 Mac・RStudioでRMeCabパッケージがうまく動作せず、 形態素解析の部分はターミナルからR実行した結果を使っています。

形態素のワードクラウドに関するイントロダクション

ワードクラウド(Wordcloud)は、 テキスト中の語彙の出現頻度をもとに、 テキスト中の頻出語を視覚的に表現する手法の1つです。 各語彙を頻度に比例する大きさで表す、視覚表現をとっています。

さて、「ワードクラウドをどのように作るか?」ということですが、

  1. テキストの取得・前処理

  2. テキストの形態素解析

  3. ワードクラウドによる形態素頻度の可視化

の3ステップになります。

「1」のテキストの取得と前処理が終われば、 次に、「2」でテキストをそれぞれ語彙に分解する必要があります。 いわゆる、形態素解析のステップです。 具体的には、RMeCab + neologd辞書を使って、形態素解析していきます。

「1. テキストの取得・前処理」と「2. テキストの形態素解析」の具体的な実行方法については、 過去の記事で大変分かりやすくまとめていますので、そちらを参照してください。

skume.net

skume.net

さっそく、「3. ワードクラウドによる形態素頻度の可視化」のステップに進みます。

はじめに、「坊っちゃん」の第1章テキストでの形態素解析済みの結果(result1.Rds or result1.txt)を GitHubからダウンロードして、R環境に読み込むところをやります。

RMeCab形態素解析済みの結果の読み込み

.Rdsをロードする場合

.RdsファイルをR環境に読み込む場合には、一度ローカルにファイルをダウンロードする必要があります。

##形態素解析の結果(result1.Rds)をR環境にロードする
#URL先の定義
Rds_path <- "https://github.com/kumeS/Blog/raw/master/R_text_analysis/R_02/result1.Rds"

#ローカルディレクトリにダウンロード
download.file(url=Rds_path, destfile = basename(Rds_path))

#Rdsの読み込み
result1r <- readRDS(file = "result1.Rds")

.txtをロードする場合

.txtファイルの場合には、WebのURLから直接読み込むことができます。 そのため、テキストの方が手数としては少なくなります。

##形態素解析の結果(result1.txt)をR環境にロードする
#URL先の定義
Txt_path <- "https://github.com/kumeS/Blog/raw/master/R_text_analysis/R_02/result1.txt"

#txtの読み込み
result1t <- read.table(file=Txt_path, header = T, sep = "\t")

結局は、以下のように同じ結果がロードされます。 好みで、ファイル読み込みの方法を決めてください。

#ヘッド表示
head(result1r)
#  parts    mor
#1  名詞 親譲り
#2  助詞     の
#3  名詞 無鉄砲
#4  助詞     で
#5  名詞   小供
#6  助詞     の

#詳細表示
str(result1r)
#'data.frame': 4870 obs. of  2 variables:
# $ parts: chr  "名詞" "助詞" "名詞" "助詞" ...
# $ mor  : chr  "親譲り" "の" "無鉄砲" "で" ...

#ヘッド表示
head(result1t)
#  parts    mor
#1  名詞 親譲り
#2  助詞     の
#3  名詞 無鉄砲
#4  助詞     で
#5  名詞   小供
#6  助詞     の

#詳細表示
str(result1t)
#'data.frame': 4870 obs. of  2 variables:
# $ parts: chr  "名詞" "助詞" "名詞" "助詞" ...
# $ mor  : chr  "親譲り" "の" "無鉄砲" "で" ...

変数result1の中身は、2列 x 4870行のデータフレーム形式のデータで、 4870個の形態素(mor列)とそれに対応する品詞(parts列)の情報が含まれています。 これは、RMeCab + neologd辞書を使った形態素解析の結果です。

ワードクラウドによる形態素頻度の可視化

さて、ワードクラウドを作成するのに使う、関連パッケージをインストールしていきます。

#wordcloud2とhtmlwidgetsのインストール
install.packages(c("wordcloud2", "htmlwidgets"))
library(wordcloud2); library(htmlwidgets)

#カラーパレット
install.packages("colorspace")
library("colorspace")

次に、ワードクラウドによる可視化を行う、実際のRコードを紹介します。

品詞情報をもとに、名詞のみ抽出した場合、形容詞や動詞も含めた場合、複合語判定の結果を含めた場合など、 複数パターンのワードクラウドを作成していきます。 また、「平仮名1文字の形態素を除く」とか、「形態素頻度の補正」とか、「出現数上位20%のみ抽出」とかの前処理もやっています。

品詞情報から「名詞」のみ抽出した場合のワードクラウド

まずはじめに、坊ちゃんの第1章のテキストデータを使って、 「名詞」のみ抽出した形態素の結果 に対して、ワードクラウドを実行します。

#Rdsのロード結果を使用します。
result1 <- result1r

#品詞から「名詞」のみ抽出した場合
result2 <- result1[result1$parts %in% c("名詞"),]

#品詞の集計
table(result2$parts)
#名詞 
#1292

#平仮名1文字の形態素を除く
result3 <- result2[!base::grepl(pattern="^[あ-ん]$", result2$mor),]

#集計
Dat <- data.frame(mor=names(table(result3$mor)),
                  Freq=as.numeric(table(result3$mor)))

#頻度1回を除く
Dat <- Dat[Dat$Freq != 1,]

#形態素頻度を0~1の値に補正
Dat$FreqR <- round(Dat$Freq/max(Dat$Freq), 4)

#出現数上位20%のみ抽出
res <- quantile(Dat$FreqR, probs = seq(0, 1, 0.1))
Dat0 <- Dat[Dat$FreqR > as.numeric(res[9]),]

#抽出の結果
dim(Dat0)
#[1] 38  3

#表示
head(Dat0[order(Dat0$Freq, decreasing = T),])
#     mor Freq  FreqR
#12  おれ   43 1.0000
#356   清   37 0.8605
#183   兄   27 0.6279
#63  もの   26 0.6047
#115   何   25 0.5814
#255   事   22 0.5116

#wordcloud2.jsによるワードクラウド作成
Dat1 <- wordcloud2(Dat0[,c("mor", "FreqR")],
                   shape="circle",
                   size=0.8,
                   gridSize =  1,
                   fontFamily="YuGothic",
                   color = "random-light",
                   ellipticity = 0.75,
                   backgroundColor = "black")

#ワードクラウド結果の表示
Dat1

#ワードクラウドの保存
htmlwidgets::saveWidget(Dat1, "./WordCloud_Dat0.html", selfcontained = F)
20220815011642
クリックすると、結果がブラウザで表示されます。

簡単にRコードの解説をします。 まず、%in%を使って、名詞だけ抽出しています。 次に、平仮名1文字の形態素を除き、単語のリストとカウントを作成しています。 そして、各前処理を行い、各単語、それぞれに対応する頻度と補正頻度を持つデータフレームを作成しています。

最終的に、ワードクラウドを実行する、wordcloud2関数には、 形態素とその補正頻度「Dat0[,c("mor", "FreqR")]」をINPUTとして与えています。 また、オプションとしては、shape、fontFamily、ellipticity、backgroundColorなどを設定しています。

品詞情報から「名詞」「形容詞」「動詞」を抽出した場合のワードクラウド

坊ちゃんの第1章のテキストデータを使って、 「名詞」と「形容詞」と「動詞」を抽出した形態素の結果 に対して、ワードクラウドを実行します。

#Rdsのロード結果を使用します。
result1 <- result1r

#品詞から「名詞」と「形容詞」と「動詞」を抽出
result2 <- result1[result1$parts %in% c("名詞", "形容詞", "動詞"),]

#品詞の集計
table(result2$parts)
#形容詞   動詞   名詞 
#    87    775   1292

#平仮名1文字の形態素を除く
result3 <- result2[!base::grepl(pattern="^[あ-ん]$", result2$mor),]

#品詞の集計
table(result3$parts)
#形容詞   動詞   名詞 
#    87    774   1255

#集計
Dat <- data.frame(mor=names(table(result3$mor)),
                  Freq=as.numeric(table(result3$mor)))

#頻度1回を除く
Dat <- Dat[Dat$Freq != 1,]

#形態素頻度を0~1の値に補正
Dat$FreqR <- round(Dat$Freq/max(Dat$Freq), 4)

#出現数上位20%のみ抽出
res <- quantile(Dat$FreqR, probs = seq(0, 1, 0.1))
Dat0 <- Dat[Dat$FreqR > as.numeric(res[9]),]

#抽出の結果
dim(Dat0)
#[1] 52  3

#表示
head(Dat0[order(Dat0$Freq, decreasing = T),])
#     mor Freq  FreqR
#57  する   89 1.0000
#173 云う   50 0.5618
#16  いる   46 0.5169
#22  おれ   43 0.4831
#98  なる   37 0.4157
#542   清   37 0.4157

#wordcloud2.jsによるワードクラウド作成
Dat1 <- wordcloud2(Dat0[,c("mor", "FreqR")],
                   shape="circle",
                   size=0.8,
                   gridSize =  1,
                   fontFamily="YuGothic",
                   color = "random-light",
                   ellipticity = 0.75,
                   backgroundColor = "black")

#ワードクラウド結果の表示
Dat1

#ワードクラウドの保存
htmlwidgets::saveWidget(Dat1, "./WordCloud_Dat1.html", selfcontained = F)
20220815012628
クリックすると、結果がブラウザで表示されます。

品詞情報から「名詞」と「複合語」を抽出した場合のワードクラウド

坊ちゃんの第1章のテキストデータを使って、 「名詞」と「複合語」を抽出した形態素の結果 に対して、ワードクラウドを実行します。

#Rdsのロード結果を使用します。
result1 <- result1r

#複合語抽出: Compound_calc()関数を使います
source("https://raw.githubusercontent.com/kumeS/Blog/master/R_text_analysis/R_03/Compound_calc.R")

#実行
result1c <- Compound_calc(result = result1)

#品詞から「名詞」と「形容詞」と「動詞」を抽出
result2 <- result1c[result1c$parts %in% c("名詞", "複合語"),]

#品詞の集計
table(result2$parts)
#複合語   名詞 
#   137   1003

#平仮名1文字の形態素を除く
result3 <- result2[!base::grepl(pattern="^[あ-ん]$", result2$mor),]

#品詞の集計
table(result3$parts)
#複合語   名詞 
#   137    972

#集計
Dat <- data.frame(mor=names(table(result3$mor)),
                  Freq=as.numeric(table(result3$mor)))

#頻度1回を除く
Dat <- Dat[Dat$Freq != 1,]

#形態素頻度を0~1の値に補正
Dat$FreqR <- round(Dat$Freq/max(Dat$Freq), 4)

#出現数上位20%のみ抽出
res <- quantile(Dat$FreqR, probs = seq(0, 1, 0.1))
Dat0 <- Dat[Dat$FreqR > as.numeric(res[9]),]

#抽出の結果
dim(Dat0)
#[1] 30  3

#表示
head(Dat0[order(Dat0$Freq, decreasing = T),])
#     mor Freq  FreqR
#12  おれ   42 1.0000
#367   清   37 0.8810
#71  もの   26 0.6190
#189   兄   25 0.5952
#121   何   23 0.5476
#264   事   22 0.5238

#wordcloud2.jsによるワードクラウド作成
Dat1 <- wordcloud2(Dat0[,c("mor", "FreqR")],
                   shape="circle",
                   size=0.8,
                   gridSize =  1,
                   fontFamily="YuGothic",
                   color = colorspace::rainbow_hcl(nrow(Dat0), c = 60, l = 60),
                   ellipticity = 0.75,
                   backgroundColor = "black")

#ワードクラウド結果の表示
Dat1

#ワードクラウドの保存
htmlwidgets::saveWidget(Dat1, "./WordCloud_Dat2_1st.html", selfcontained = F)
20220815012817
クリックすると、結果がブラウザで表示されます。

坊ちゃん第2章以降のワードクラウドによる可視化

次に、第2章以降でワードクラウドを試してみた結果を紹介します。

坊ちゃんの第2章のテキストデータを使って、 「名詞」と「複合語」を抽出した形態素の結果 に対して、ワードクラウドを実行します。

上記の実行が冗長なので、 ワークフローを自作関数化して(wdCloud.R関数という名前)、 実行を簡略化しています。

坊ちゃん第2章のワードクラウド

#txtファイルを読み込んで使用します。
Txt_path <- "https://github.com/kumeS/Blog/raw/master/R_text_analysis/R_02/result2.txt"
result1 <- read.table(file=Txt_path, header = T, sep = "\t")

#wdCloud.Rのロード
source("https://raw.githubusercontent.com/kumeS/Blog/master/R_text_analysis/R_04/wdCloud_test.R")

#ワードクラウド実行
Dat1 <- wdCloud.R(result1)

#ワードクラウド結果の表示
Dat1

#ワードクラウドの保存
htmlwidgets::saveWidget(Dat1, "./WordCloud_Dat2_2nd.html", selfcontained = F)
20220815013043
クリックすると、結果がブラウザで表示されます。

坊ちゃん第3章のワードクラウド

続いて、坊ちゃん第3章のワードクラウドです。

#txtファイルを読み込んで使用します。
Txt_path <- "https://github.com/kumeS/Blog/raw/master/R_text_analysis/R_02/result3.txt"
result1 <- read.table(file=Txt_path, header = T, sep = "\t")

#wdCloud.Rのロード
source("https://raw.githubusercontent.com/kumeS/Blog/master/R_text_analysis/R_04/wdCloud_test.R")

#ワードクラウド実行
Dat1 <- wdCloud.R(result1)

#ワードクラウド結果の表示
Dat1

#ワードクラウドの保存
htmlwidgets::saveWidget(Dat1, "./WordCloud_Dat2_3rd.html", selfcontained = F)
20220815013128
クリックすると、結果がブラウザで表示されます。

同じなので、以下省略。。。

https://kumes.github.io/Blog/R_text_analysis/R_04/WordCloud_Dat2_4th.html

https://kumes.github.io/Blog/R_text_analysis/R_04/WordCloud_Dat2_5th.html

https://kumes.github.io/Blog/R_text_analysis/R_04/WordCloud_Dat2_6th.html

https://kumes.github.io/Blog/R_text_analysis/R_04/WordCloud_Dat2_7th.html

https://kumes.github.io/Blog/R_text_analysis/R_04/WordCloud_Dat2_8th.html

https://kumes.github.io/Blog/R_text_analysis/R_04/WordCloud_Dat2_9th.html

https://kumes.github.io/Blog/R_text_analysis/R_04/WordCloud_Dat2_10th.html

https://kumes.github.io/Blog/R_text_analysis/R_04/WordCloud_Dat2_11th.html

坊ちゃん第1〜11章のワードクラウドの結果まとめ

 

過去に坊ちゃんを読んだことなるような、無いような、内容はほとんど覚えてませんので、 ワードクラウド見た感じの印象をツラツラと書いてます。

全体を通して、「おれ」はよく出てきてますね。「おれ」中心の話なんでしょうね。

第1章は、「おれ」「清」「兄」とかが目立つ。清が兄なんかな、分からんけど。

第2章はなんか和茶和茶してますね。中学校が舞台なのか、教師とか校長とか、山嵐も出てくるようです。

第3章のクラウドはやや控えめですね。学校の隣に東京とあり、誰か転校生がやってくるのか? ワードクラウドでは、語彙の位置関係は関係ないですけども。

第4章は「バッタ」と「宿直」、第5章は「赤シャツ」、第6章は「山嵐」、第7章は「模試」が目立ちます。 模試とあるので試験勉強してのかな。山嵐ってなんだっけ?風のこと?柔道の技のこと?

第8章で、もう一度「赤シャツ」が出てくる、「婆さん」「月給」もでてくるようだ。 「月給」とか、なんだか現実的な話っぽいね。

さらに、第9章で「うらなり君」と「送別会」、第10章で「喧嘩」、 第11章ではまたも「山嵐」で、「新聞屋・新聞」が新しく出てくる。 転校生の「うらなり君」喧嘩したのか?最後に、新聞バイトしてるのか?謎は深まるばかり。。   ワードクラウドを見ただけだと、テキトーな考察になってしまいましたね。

まとめ

今回は、テキストマイニングの一例として、 坊ちゃんの各章のテキストデータをワードクラウドで可視化してみました。

結構簡単にワードクラウドを実行できるので、 色々な条件で見方を変えて、ぜひ色々と試してみてください。

参考資料

www.karada-good.net

Wordcloud2 introduction

福岡の旅 x 離島猫 x レトロ港

先日は、久々に福岡に行ってきました。まぁ、4ヶ月ぶりくらいなんですが、、、

これまでは学会とか研究会とかの参加で訪れていたためか、純粋な観光で行ったのはこれが初めてかもです。


初日は、小郡市(「おごおりし」だそうです)というところで、伝説のラーメンを食しました。本場のとんこつラーメンはうまし!でした。


そのあとは、ラーメン屋の隣に位置するメルフェン境、カエルケロケロ寺を探索してきました。今夏の最高気温を更新した日とあってか、鈴の音が涼しくもあり、暑さで脚元はフラフラでもあり・・・


夕暮れあたりには福岡市内に移動し、海近くのホテルで滞在。潮風が漂う海の景色を眺めていました。


次の日は、博多湾に浮かぶ離島を探索。島の高台から見える、海は立ち止まってしまうほど、非常に綺麗でした。


また、とある島では、猫、ねこ、ネコに大量遭遇。おそらく、1年分くらいの猫を見た気がします。


最後には、中洲でボラれたことで凹んで嫌気がさし、福岡市内から少し足を伸ばすことに。 そして、レトロ港に出会しました。。 色んなものが眩しい港町でした。

中洲の屋台街の物悲しさ、、身包みをやや少なくする真夏のイヤラシサを除けば、、真夏の『福岡』の旅は大満足でした。

Mac版wgetコマンドのプログレス・ログ詳細を表示させないTips

はじめに - 問題提起

wgetコマンドを使ってたら、そのデフォルトログが冗長で嫌になりませんか??

今回の記事では、その対策として、ログを非表示にする方法を紹介します。

wgetの使い方の基本は、過去の記事を参考のこと。

skume.net

skume.net

Googleドライブからのファイルダウンロード

Googleドライブに置いてあるZipファイル(10MBくらい)のwgetダウンロードを試しにやってみます。

まずは、Macのターミナルを起動します。

#ターミナル上でのwgetダウンロード実行
wget -r "https://drive.google.com/uc?export=view&id=1MRdPXshHQshYszUpcLvNf-FjCkPCGdLo"
#OR
wget -r "https://drive.google.com/uc?export=view&id=1MRdPXshHQshYszUpcLvNf-FjCkPCGdLo" -O test.zip

#-r, --recursive: ダウンロードを再帰的に行う
#-O: 複数ファイルを一つに連結して別名保存

##(以下、出力ログ)
#--2021-08-09 17:28:54--  https://drive.google.com/uc?export=view&id=
#1MRdPXshHQshYszUpcLvNf-FjCkPCGdLo
#drive.google.com (drive.google.com) をDNSに問いあわせています... 
#2404:6800:4004:80a::200e, 172.217.161.78
#drive.google.com (drive.google.com)|2404:6800:4004:80a::200e|:443 
#に接続しています... 接続しました。
#HTTP による接続要求を送信しました、応答を待っています... 302 Moved Temporarily
#場所: https://doc-14-4o-docs.googleusercontent.com/docs/securesc/
#ha0ro937gcuc7l7deffksulhg5h7mbp1/48tuvarfgjr41bb25vc44vm2nots5rtq/
#1628497725000/15705379774912454503/*/1MRdPXshHQshYszUpcLvNf-FjCkPCGdLo?e=view [続く]
#警告: HTTPはワイルドカードに対応していません。
#--2021-08-09 17:28:55--  https://doc-14-4o-docs.googleusercontent.com/
#docs/securesc/ha0ro937gcuc7l7deffksulhg5h7mbp1/48tuvarfgjr41bb25vc44vm
#2nots5rtq/1628497725000/15705379774912454503/*/1MRdPXshHQshYszUpcLvNf-
#FjCkPCGdLo?e=view
#doc-14-4o-docs.googleusercontent.com (doc-14-4o-docs.googleusercontent.com) 
#をDNSに問いあわせています... 2404:6800:4004:821::2001, 172.217.25.97
#doc-14-4o-docs.googleusercontent.com (doc-14-4o-docs.googleusercontent.com)|
#2404:6800:4004:821::2001|:443 に接続しています... 接続しました。
#HTTP による接続要求を送信しました、応答を待っています... 200 OK
#長さ: 特定できません [application/zip]
#`drive.google.com/uc?export=view&id=1MRdPXshHQshYszUpcLvNf-FjCkPCGdLo' に保存中
#
#drive.google.com/u     [       <=>    ]  12.60M  8.62MB/s 時間 1.5s     
#
#2021-08-09 17:28:57 (8.62 MB/s) - `drive.google.com/uc?
#export=view&id=1MRdPXshHQshYszUpcLvNf-FjCkPCGdLo' へ保存終了 [13208722]
#
#終了しました --2021-08-09 17:28:57--
#経過時間: 3.7s
#ダウンロード完了: 1 ファイル、13M バイトを 1.5s で取得 (8.62 MB/s)

うーん、ログが多い。これは冗長だ。。。問題だ!!

デフォルトログ表示への対策

結局のところ、wgetのオプションで、-q-dあるいは-oオプションを付与すると解決します。

以下に、それらのオプションについて説明します。

  • -q, --quiet: 処理結果の非表示

  • -b, --background: 処理をバックグラウンド実行。また、-oオプションが非指定の場合、処理結果をwget-logに書き込む。

  • -o, --output-file=: 処理結果の別名保存

次に、それらのオプションを付けて、wgetの実行をしてみます。

#q実行
wget -q "https://drive.google.com/uc?export=view&id=1MRdPXshHQshYszUpcLvNf-FjCkPCGdLo" -O test.zip
#(何も出ない)

#b実行
wget -b "https://drive.google.com/uc?export=view&id=1MRdPXshHQshYszUpcLvNf-FjCkPCGdLo" -O test.zip
#バックグラウンドで継続します、pidは 2597。
#出力を `wget-log' に書き込みます。

#o実行
wget "https://drive.google.com/uc?export=view&id=1MRdPXshHQshYszUpcLvNf-FjCkPCGdLo" -o test.zip
#(何も出ない、、一時ファイル?)

-bだと出力ログのメッセージが少しでるものの、 これらオプションをつけることで、ログ表示はほぼ無くなりました!勝利!!

あと、-oだと、一時ファイルが現在のディレクトリに残るのが気になりますが。

以上をまとめると、ターミナル実行の場合には、-q(あるいは-b)をオプションとして付与しておくのが無難と思います。

R言語でのwget実行

次に、R環境での実行例を示します。

これで、Rコンソールでの実行エラーも減るので助かります。

#Rでの実行例
system('wget -q "https://drive.google.com/uc?export=view&id=1MRdPXshHQshYszUpcLvNf-FjCkPCGdLo" -O test.zip')

#OR

system('wget -b "https://drive.google.com/uc?export=view&id=1MRdPXshHQshYszUpcLvNf-FjCkPCGdLo" -O test.zip')

R環境だと、wget実行中に、Rコンソールが占有されるので、 -bでバックグラウンド実行しておくのがオススメです。

参考資料

www.itsenka.com

【R言語/速度論/数値計算】10倍量のヨーグルトを作りながら、乳酸菌の増殖曲線を考えてみた件

この記事は、2020年6月13日記事のアップデート版です。

はじめに

最近、ヨーグルト作りにハマっています。

ヨーグルト菌、いわゆる、0.1Lの飲むヨーグルトから、10倍量(1L)のヨーグルトを作っています。 飲むヨーグルト1本は、だいたい100円くらいなので、10倍量で1000円となって、必要経費を引いても、だいたい700円弱お得になります。。また、ヨーグルトの酸味や固さを自分の好みの具合に調整できます。

今回、「ガセリ菌を使った10倍量(1L)ヨーグルトの作製法」と「菌の増殖曲線」について検討したことを、それを実験ぽく紹介することにします*1

方法 (1): 10倍量ヨーグルト作製のマテリアル

  • ガセリ菌の種菌(= 某メグミルク製のガセリ菌SP株の飲むヨーグルト)
  • 培地(= 某メグミルク製の牛乳)
  • グルコース(= 三温糖)
  • 1L プラビーカー(= ヨーグルトメーカーの容器)
  • スパチュラ(いわゆる、スプーンてやつ)
  • 温度・時間可変の恒温器(いわゆる、ヨーグルトメーカー)

方法 (2): 10倍量ヨーグルトの作製プロトコール

1. 1L プラビーカーに対して、スパチュラ(いわゆる、スプーンですね)3-4杯分のグルコース(要するに、砂糖です)を添加します。おおよそ、15-20gくらい加えます。

2. ガセリ菌の種菌を全量加えて、振盪混合します(いわゆる、振って混ぜます)。

3. 培地(いわゆる、牛乳です)で、1L目盛までメスアップします(つまりは、牛乳を900mL入れる)。

4. 蓋を軽めに閉めて、35℃の恒温器(いわゆる、ヨーグルトメーカー)内で静置します。オーバーナイト・インキュベート(一晩培養)します。

ガセリ菌の場合、だいたい、12-15時間ほどインキュベートすると、良い感じの固形ヨーグルトが出来あがります。

5. 4℃で保管して(冷蔵庫で保管)して置いておきます。そして、官能評価を行います。いわゆる、「食う」というやつです。

結果: 種菌の増殖曲線の検討

微生物の増殖曲線について調べてみると、 ゴンペルツの増殖曲線 (Gompertz growth curve)というのがあるようです。

もともと、1832年に、Benjamin Gompertz によって、人間の死亡率を推定するためのモデル式として、「ゴンペルツ関数(Gompertz function)」が提案されたことを起源としているようです。

Gompertz, B. "On the Nature of the Function Expressive of the Law of Human Mortality, and on a New Mode of Determining the Value of Life Contingencies." Phil. Trans. Roy. Soc. London 123, 513-585, 1832.

それから、さらに100年後、1932年に、Charles Winsorによって、微生物の成長過程を説明するための方程式として使用されたようです。

www.pisces-conservation.com

ゴンペルツの増殖曲線の一般式は、以下のように定義されます。



A \cdot exp(-B \cdot exp(-k (t - I)))


このとき、

A は、飽和(上限)時の値, A > 0

B, k は、成長速度係数 (Growth velocity factor), k > 0, B > 0

I は、変曲点

t は、相対時間

をそれぞれ意味します。

今回、この増殖曲線について、 成長速度係数 k を、0.75、0.5、0.3、0.2、0.1 の5パターンに設定して、 R上で数値計算のシミュレーションしてみました。

## 固定パラメータの設定
# A: 飽和(上限)時の値, A > 0 
A <- 10
# B: 成長速度係数, B > 0
B <- 7
# I: 変曲点
I <- 5
# t: 相対時間
t <- seq(0, 20, 0.5)

## パラメータ k の設定
# 成長速度係数 k = 0.75
k1 <- 0.75
Lt1 <- round(A*exp(-B*exp(-k1*(t-I))), 3)

# 成長速度係数 k = 0.5
k2 <- 0.5
Lt2 <- round(A*exp(-B*exp(-k2*(t-I))), 3)

# 成長速度係数 k = 0.3
k3 <- 0.3
Lt3 <- round(A*exp(-B*exp(-k3*(t-I))), 3)

# 成長速度係数 k = 0.2
k4 <- 0.2
Lt4 <- round(A*exp(-B*exp(-k4*(t-I))), 3)

# 成長速度係数 k = 0.1
k5 <- 0.1
Lt5 <- round(A*exp(-B*exp(-k5*(t-I))), 3)

## データフレームに変換
Data <- data.frame(t, Lt1, Lt2, Lt3, Lt4, Lt5)

par(family="HiraKakuProN-W6", lwd=1.5, 
    cex=1, mgp=c(2.5, 1, 0), mai=c(0.75, 0.75, 0.8, 0.25))
plot(Data$t, Data$Lt1, type="l", col="#0078B9",
     main=expression(A %.% exp(-B %.% exp(-k*(t - I)))),
     xlab="time (相対時間)",ylab="Values of Gompertz growth model")
points(Data$t, Data$Lt2, type="l", col="#FF7700")
points(Data$t, Data$Lt3, type="l", col="#00A400")
points(Data$t, Data$Lt4, type="l", col="#E90017")
points(Data$t, Data$Lt5, type="l", col="#9D62C3")
20200528071149
ゴンペルツ増殖曲線の作図

さらに、plotlyパッケージを用いて、増殖曲線を可視化しました*2

if(!require("plotly")){install.packages("plotly")}; library(plotly)
fig <- plot_ly(Data, x = ~t)
fig <- fig %>% add_trace(y = ~Lt1, name = 'Lt1', mode = 'lines')
fig <- fig %>% add_trace(y = ~Lt2, name = 'Lt2', mode = 'lines')
fig <- fig %>% add_trace(y = ~Lt3, name = 'Lt3', mode = 'lines')
fig <- fig %>% add_trace(y = ~Lt4, name = 'Lt4', mode = 'lines')
fig <- fig %>% add_trace(y = ~Lt5, name = 'Lt5', mode = 'lines')
fig <- fig %>% layout(xaxis = list(title = 'time'),
                      yaxis = list(title = 'Values of Gompertz growth model'))
fig
20200528071149
Poltlyによるゴンペルツ増殖曲線の作図

まとめ

カゼリ菌の10倍量ヨーグルトは、だいたい12時間ほどで固まります。 シミュレーションの結果を考慮すると、種菌はLt1Lt2くらいスピードで増殖して、 牛乳に含まれるカゼインの酸凝固が進んでいくのではと推測されました。

参考文献

Gompertz sigmoidal model の参考文献

Rui Yuan et al, Modified Gompertz sigmoidal model removing fine-ending of grain-size distribution, 2019.

Modified Gompertz sigmoidal model removing fine-ending of grain-size distribution

journals.plos.org

www.peko-step.com

stat.ethz.ch

www27.cs.kobe-u.ac.jp

*1:某ガセリ菌株の方が、某R1菌株よりも、長期培養した際の酸味が少ないので、10倍量ヨーグルト作製時には扱いやすいと思います

*2:どちらも、plotly図にリンクしている