温度の鉛直分布の表示を試みる - その 2 - ライブラリ忘備録編 -

　今回は、前回試した気象庁の MSM から一地点のみをわざわざ取り出して鉛直方向に表示するコードの忘備録編です。

　MSM の予報は 5 km メッシュで作成されています。よって、見たい地点の予報データを、その地点の緯度経度を含む 5 km メッシュのデータと見なして、それを取り出すこととしました。

Ⅰ. GRIB 2 から必要な範囲のデータを取り出す

　これまで予報データは、水平分布を把握するべく地図上にしか図示しませんでした。そのため、描画範囲は Cartopy の set_extent () で緯度経度を指定していました。つまり、データは全部読み込んでいましたが、描画する範囲のみを限定していた、ということになります。

Ⅰ-① Pygrib の範囲指定によるデータ取得

　鉛直方向の描画では Cartopy は使いませんので、 MSM から必要な範囲のみの GRIB2 データを取り出せれば、と思い、Pygrib の説明を読んでみたところ、ありました、ありました。

>>> data, lats, lons = grb.data(lat1=20,lat2=70,lon1=220,lon2=320)

上記の data () を使って、緯度経度を指定、そのデータを取得できます。上記は一例ですが、lat1 = 最小緯度、lat2 = 最大緯度、lon1 = 最小経度、lon2 = 最大経度で囲まれる領域のデータです。

jswhit.github.io

Ⅰ-② 気象庁の MSM のメッシュの確認

　では、5 km は緯度経度にすると何度位になるのでしょう。地球 1 周 40000 km として計算でもいいのですが、もっときちんとした数字があるだろうと思い、気象庁の配信資料に関する技術資料を確認します。

www.data.jma.go.jp

　例えば No.500 で確認すると、緯度 0.0625 度、経度 0.05 度とあります。

Ⅰ-③ MSM GPV から一メッシュのデータを取り出すプログラム

　ということで、データの格子点が一つは確実に入り、二つは入らないようにするため、緯度経度幅 0.08 度ということで（単にざっくり決めただけですが）、コードを書いてみます。

import pygrib

# Set parameter 
dat = "0923"　　# 9 月 23 日 12Z
tm ="12"
lat=35.764　　　# 見たい場所の緯度経度
lon=140.392
dev=0.04　　　  # 緯度経度幅 0.08 度

lvls=[1000,975,950,925,900,850,800,700,600,500,400,300,250,200,150,100]　　# 気圧面

filenm = "Z__C_RJTD_2021"+dat+tm+"0000_MSM_GPV_Rjp_L-pall_FH00-15_grib2.bin"
grbs = pygrib.open(filenm)

def get_value(req_data):
    val = grbs.select(forecastTime=0, level=lvl)[req_data]
    val_pos = val.data(lat1=lat-dev, lat2=lat+dev, lon1=lon-dev, lon2=lon+dev)
    Val_pos=float(val_pos[0])
    val_set.append(Val_pos)

val_set=[]
for lvl in lvls:    
    req_data = 3　　# 温度
    get_value(req_data)

print(val_set)

　上記は MSM GPV の気圧面データを利用するものなので、地上面や他の GPV ではファイル名やデータの種類はちょっと変わると思います。さて、上記のプログラムを動かしてみると以下のようになりました。無事 16 の基準気圧面の温度予報値が取得できています。

[295.95611572265625, 296.2651062011719, 296.1998291015625, 297.4837951660156, 295.5482177734375, 291.1206359863281, 286.7886962890625, 280.7910461425781, 274.5665283203125, 266.0440979003906, 255.32080078125, 239.46371459960938, 228.833740234375, 217.5061492919922, 211.73158264160156, 205.81198120117188]

　ひとまず目標は達成されました。そして、これは個人的希望ですが、鉛直方向の温度分布を見るのであれば、気温の他に露点が欲しいと思いました。もっとも、MSM GPV に気温データはあるのですが、露点データはありません。ですが、相対湿度データはあるので、ここから、露点を求めたいと思います。

Ⅱ. 露点を計算する（と言っても、ライブラリを使うだけです）

　普段目にする、例えば空港気象情報には、気温と露点が並んで表示されているので、気温を測るように、露点も何らかの手段で難なく測れるものと思っていたのですが、検索をかけてみると、そう簡単というわけではないようです。事実、何らかの計算や、理論式からの換算表などが用いられることが多いようです。そこで、今回理屈の理解はひとまず断念し（笑）、その値を求める方法のみを確認することとします。　そして、こういう時こそ、python ではライブラリの出番です。

Ⅱ-① Metpy により露点を求める

　前回書いたように、Metpy というライブラリは気象関連の計算描画がいろいろ出来るライブラリなので、その簡易版をわざわざ自分で作るのは、本来であれば、単にプログラミングの練習程度にしかなりません。ですが、今のところは描画の段階でなぜかエラーが出るので、今回は、その計算機能部分だけ拝借することとします。

　こういうライブラリであれば当然なのでしょうが、気温と相対湿度から露点を計算する機能があるのはさすが素晴らしい限りです。

unidata.github.io

　Metpy の面白い所は、物理量に対する姿勢が厳格というか、物理量に次元（単位）を付けて計算するところです。実にもっともな話ではあるのですが、電卓の延長のように計算しようとすると、単位分だけコードが長く、ややくどい感は拭えません。おそらく、次元をコードに組み込めた方が利点のある人がユーザーの多くを占めるのでしょうね。

　次元をコードに組み込むことについて書かれたチュートリアルもバッチリあります。

unidata.github.io

Ⅱ-② 計算してみましょう

　では、あるまだ穏やかめな暑い夏の日の露点を求めてみましょう。

　気温約 32 ℃、相対湿度 85 ％ の時の露点を求めるコードは以下の様になります。

from metpy.units import units
import metpy.calc as mpcalc

dew_pos = mpcalc.dewpoint_from_relative_humidity(305 * units.kelvin, 0.85)
dewp=dew_pos/units.kelvin -273.15

print(dewp, type(dewp))

　上の結果は以下の通りです。

29.01570497604706 dimensionless <class 'pint.quantity.build_quantity_class.<locals>.Quantity'>

　露点約 29 ℃。湿度の暑い夏は、数度気温が下がったところで、大気がすぐ飽和するので、ベタベタ感がなくなることがないことが分かります。日中多少気温が上がったところで、秋の空気はやはり素敵ですね。

　コード中、units.kelvin を掛けたり割ったりしている点が、先に触れた通り、温度という物理量の次元部分です。露点 29.015.... は無次元（dimensionless）、数値だけですが、長いクラス名が付けられています。細かいことは分かりませんが（笑）、Ⅰ-③ で求めた温度のリスト val_set と同等に扱えます。

　前回の記事を書いた時には、コードの忘備録を残そうと意気込んでいたのですが、少し時間が経っただけで、何を残しておこうと思ったのか、そもそも何で忘備録を作ろうと思ったのかすら見事に怪しくなりました。

　あとは取り出した GRIB2 データを図示するだけなので、描画ライブラリに詳しい方の方が色々とできるかと思いますが、自分でも使い勝手の良い図をちょっと模索中なので、また後日、続きのまとめを書くかもしれません。

　まとめをほとんど作らない人生を送っている身には、忘備録を作ることは、やはりなかなか骨が折れます。

　本日も最後までお付き合いどうもありがとうございました。

温度の鉛直分布の表示を試みる - その 1 - 思い立ったら、まず結果編

何もいらない幸せ

　朝方は曇っていました。関東では、北東の風が吹くと、鹿島灘からの湿った空気が流れ込み、低い雲が広がるとよく言われます。

　もっとも、昼前にはほとんど雲も消え、澄んだ秋の青空が広がりました。テレビ越しではありますが、駅伝のある出雲でもきれいな青空が広がっていました。何もしなくても幸せな気分になってしまう季節の到来です。

　ところで、低い雲の上は実際にはどうなっていたのでしょう。先日の記事『雲の底 ～ おそらく大気の不連続面 ～』の中でもリンクを張りましたが、雲は低、中、高層雲の三層に分類します。もっとも、実際の高さがどの位かが伝えられることはそれほど多くはありません。私たちは普通、地面の高さで生活しているので、雲の高さが問題となることはさほどないからでしょう。

　でも、低い雲が切れ始め、その向こうに青空が広がっているのを目にすると、低い雲が広がっていたその時、雲の上はどうなっていたのだろう、と思う時があります。今日はそんな一日でした。おそらく平和で暇な休日の一時だったからでしょう（笑）。

気象庁の高層気象観測

　地上の気象観測に比べれば、その数は圧倒的に少ないですが、高層の気象観測も全国 16 ヵ所で一日二回行われす。

www.jma.go.jp

　日本全土を程よく網羅し、気球に観測機器を付けて上空に放球するため、自分の住んでいる場所近くのデータが手に入ることはあまりなさそうです。が、いつも使わせて頂いている数値予報から自分の見たい地点のデータを取り出せば、その時点の状況は窺い知れます。

　基準気圧面における予報値が広大な領域にわたって小さなメッシュで整えられた数値予報から、わざわざ一地点のデータのみを取り出すというのは何とももったいない話ですが、高層気象観測の行われる地点以外でエマグラムのようなデータを見てみたい、というのは、確かに時々思うことでした。

MetPy、まだ直っていなかった、"name 'ctables' is not defined"

　python はライブラリが本当に充実しているので、検索をかければ、それらしいことが出来るものはたいてい見つかります。もっとも、ひとたびやりたいことが解決すると、細かい文法はどんどん忘れ、挙句の果ては、何で自分がこういうコードを書いたのかすらよく分からない blackbox と化します（まあたいていは、examples をいじり倒しただけですが…）。

　ということで、MetPy というライブラリにエマグラムの親戚、Skew-T 図が描けた覚えがあったので、気象庁の MSM 予報から Pygrib で一地点の気温を取り出し、描画を試みました。すると、

　　”name 'ctables' is not defined”

というエラーが発生しました。これ、何か前に見たなあ、と思っていたらそのはず、一年ほど前に、Basemap から Cartopy に引っ越しをした時のあれこれの記事の一つ、『MetPy と Cartopy の微妙な関係』で書いたエラーでした。まだ直っていなかったんですね。検索をかけても、その後の進展もあまりないようです。あまり使われない機能なのでしょうか…。

　昨年記事を書いた時には、ライブラリであるものをわざわざ自分で作るのも面倒かな、と思ったのですが、今回は予報値を見てみたい、という欲望が勝りました。ということで、エマグラムとまでは言えませんが、気温、露点温度と気圧の関係図を作ることとしました。

ひとまず結果から

　普段は、変数部分だけいじって、エラーや Warning が出ない限りは、動いてくれるのであればコードを見返すことは、残念ですが全くありません。その結果、細かい使い方を、あきれるほどきれいに忘れてしまうことが、この度改めて分かりました。ということで、自分のリマインドのためにも、まとめの記事を書こうと思ったのですが、長くなりそうなので、次回にすることにしました（笑）。

　ひとまず今日の結果から。使用データは、2021 年 9 月 23 日 12Z (21L) の気象庁 MSM 予報、場所は成田空港、先日の記事『気象予報をもう一度見直す - 成田空港の霧』で書いた通り、予報にない霧が急激に発生した時点の予報初期値です。

　エマグラムとまでは言えませんが、気温と露点の高度変化です。描画は好みで細部をいくらでもいじれそうですが、ひとまず今日はこの位で終了します。

　確かに、地表付近でのみ雲が発生しやすい状況であったことが窺えます。

　明日の祝日は、すでに夏休み中に使われていたんですね。楽しみは先が良いのか、後に取っておくべきなのか。正直、楽しみの種類によりますね。

　本日も最後までお付き合いどうもありがとうございました。

（改題しました、15Oct.2021）

（旧）簡易エマグラムを作る - その 1 - 思い立ったら、結果編

（現）温度の鉛直分布表示を試みる - その 1 - 思い立ったら、まず結果編

年降水量偏差と太陽黒点数の関係 ( は多分ほぼない )

降水量と太陽活動

　先日の記事で、素人考えからすると、太陽活動の強弱で気候の変化が説明できる部分があるのではないかと書きました。

　無線の世界では『サイクル 24 が終わりかけで…』等と言って、太陽活動によって通信に影響が出ることがかつて取りざたされていました。太陽活動周期はおよそ 11 年周期とよく耳にしますが、同じように、気象庁の年降水量偏差と太陽活動には、何か相関が見出せるのでしょうか。　

www.data.jma.go.jp

そもそも太陽活動は、どのような指標を使って考えるのか

　経済活動や災害の規模等を考える場合もそうですが、それを表わす指標や基準は色々選ぶことができます。また、選んだ指標や基準の結果への影響の程度も一様ではありません。そもそも太陽の活動周期が約 11 年というのは、何を根拠とするものだったのでしょう。

　太陽活動による効果（Wikipedia へ）は、磁場、放射、電波、爆発、宇宙線等があります。太陽については、記録された長さから、黒点の数、F10.7 太陽電波フラックス（太陽からの電波のピークに近い波長 10.7 cm の電波）等がよく使われ、黒点の数については、400 年にわたる記録があります。

commons.wikimedia.org

　太陽活動を測ることは、宇宙利用の拡大とも相まって、広範多岐にわたるようで、検索する程度では、残念ながらそのさわりすら把握できませんでした（泣）。

　いずれにせよ、太陽活動周期が約 11 年というのは、黒点数の変動周期で、今回も黒点数の変化を太陽活動周期と考えて、年降水量偏差との関係を調べてみることとしました。

太陽黒点数データ

　太陽の黒点数データも色々な研究機関が公表していますが、今回は国立天文台、三鷹太陽地上観測のデータを使わせて頂きました。1929 年から 2021 年まで日々の数値が公開されており、1929 年からの黒点相対数のグラフも見ることができます。

solarwww.mtk.nao.ac.jp

　　　Copyright 2017 Solar Science Observatory, NAOJ. All rights reserved. (Above)

しばらく眺めてみた結果…

　上のグラフの通り、太陽黒点数は確かにきれいな周期を描いて増減することが分かります。ですが、気象庁の年降水量偏差と比べると…、んー…、

　「時系列変化の状況を見てみると、これはほとんど関連無しかも…」

という残念な予想が見えた気がしました。本当は、pandas を駆使してあれこれやろうかとも思っていたのですが…。ですが、あれこれ調べた半日を一応形にするべく、めげずに可視化と相関を求めてみました。

データの準備と整形

　年降水量の偏差は csv 形式、黒点数の年平均値は txt 形式のファイルでした。対象期間は、気象庁と国立天文台に共通する期間、1929 年から 2020 年としました。そして、年降水量偏差が 1991 年から 2020 年の 30 年平均からの偏差を使っているので、黒点数に関しても、同じ期間の平均値からの偏差を用いました（あまり意味は無さそうですが）。

　ちなみに1991 年から 2020 年の黒点数の年平均個数は、49.77 個でした。

年降水量偏差と黒点数偏差の相関

　年降水量偏差と黒点数偏差の推移を視覚化した結果が以下の通りです。

　同じく、df.corr() により相関係数を求めてみると、0.115019。ほぼ無関係ですね…。やっぱり。

まとめ

　前回の記事を書きながらふと思いついた、年降水量偏差と太陽黒点数偏差の相関は、ほとんどないことが分かりました（泣、残念）。

　太陽活動を黒点数で代表させ、また太陽から届いたエネルギー（正確にはエントロピーでしょうか）が色々変遷して現れる実に多くの結果の一つが降水であることを考えれば、両者に直接的、あるいは有意な関係を見出すのは、ちょっと無理かもしれませんね。

　蛇足ながら、太陽活動と株価の相関や景気の予想を試みた例も多々あることが分かりました。是非はともかく、人間が太陽や月に何かを見出そうとするのは、因果的に世界を理解しようとする、長年にわたり培われた思考パターンの一つなのかもしれませんね。

本日も最後までお付き合いどうもありがとうございました。

日照時間で感じる季節の移り変わり - 2 - 2020 年の軽い EDA

　前回は、国立天文台の「日の出入り」データを Beautiful Soup でお借りし、それを numpy / pandas 形式で整えるところまでやりました。

　必要な日の出入りのデータをお借りすることはできるようになりましたが、そのたびにアクセスするのも先方に迷惑なので（？）、一年分まとめてお借りし、それを自分の computer に保存、その後いじくることとしました。

目次

１．2020 年一年分の日の出入りのデータをまとめて保存する

２．データの整形

３．簡単に見てみる

４．簡単に図示する

まとめ

１．2020 年一年分の日の出入りのデータをお借りし、保存する

　データをまとめ、保存する方針と形式（と言うほどのものではないのですが）は、次のようにしました。

① 2020 年を一つのフォルダへまとめる。

② サイトのページをそのまま保存するイメージ。つまり、ファイル数は月数。

③ 保存ファイル形式は、 csv 形式。

　前回は 1 ページ（つまり 1 か月分）のみ読み込むだけだったのですが、一年分まとめて読み込み、それを保存する、という形にしました。

import urllib.request
from bs4 import BeautifulSoup
import numpy as np
import pandas as pd
import os

yr = "2020"
mnths = ["01","02","03","04","05","06","07","08","09","10","11","12"]
dys = [31,29,31,30,31,30,31,31,30,31,30,31] 　# 2020 only

def gather_srss(i):
    mnth = mnths[i] 
    dy = dys[i]  
    target_url=
         "https://eco.mtk.nao.ac.jp/koyomi/dni/"+yr+"/s13"+mnth+".html"
 
    with urllib.request.urlopen(target_url) as hinodeiri:
        html = hinodeiri.read()
        soup = BeautifulSoup(html, "html.parser")
        srsss = soup.find_all('td')

    A = []
    for srss in srsss:
        a=srss.contents[0]
        A.append(a)

    B = np.array(A).reshape(dy,7)
    for m in range(dy):
        B[m,1] = B[m,1].replace(' ', '0') 
 
    s = pd.DataFrame(B)
    s.columns=['Date', 'SR_time', 'SR_dir', 'Culmination',
             'Culm_alt', 'SS_time', 'SS_dir']
 
    os.chdir(r'C:\　保存するフォルダ名 ')
    s.to_csv("　ファイル名　"+yr+mnth+".csv")
 
for i in range(12):
    gather_srss(i)

２．データの整形

　これで相手のサーバーに迷惑をかけることがなくなりました。

　保存したファイルを読み込み、一年分まとめてデータフレーム（df20）にします。そして、列ごとにデータの型式を揃え、ついでに、日照時間（Daylight_Hours）、太陽の移動角度（Travel_Angle）を付け加えました。

import os
import pandas as pd


os.chdir(r'C:\　保存したフォルダ　')
yr = "2020"
mnths = ["01","02","03","04","05","06","07","08","09","10","11","12"]
df20 = pd.DataFrame()

def read_srss(k):
    del df['Unnamed: 0']
    df['Date'] = df['Date'].replace({' 1':'01', ' 2':'02',' 3':'03',
          ' 4':'04',' 5':'05', ' 6':'06',' 7':'07',
          ' 8':'08',' 9':'09'})
    df['SR_time'] = yr+"/"+mnths[k]+"/"+df['Date']+" "+df['SR_time']
    df['SS_time'] = yr+"/"+mnths[k]+"/"+df['Date']+" "+df['SS_time']
    df['Culmination'] = yr+"/"+mnths[k]+"/"+df['Date']+
          " "+df['Culmination']
 
for k in range(12):
    df = pd.read_csv("SRSS"+yr+mnths[k]+".csv", dtype={'Date':str})
    read_srss(k)
    df20 = df20.append(df)

df20['SR_time'] = pd.to_datetime(df20['SR_time'])
df20['SS_time'] = pd.to_datetime(df20['SS_time'])
df20['Culmination'] = pd.to_datetime(df20['Culmination'])

df20['Daylight_Hours'] = df20['SS_time'] - df20['SR_time']
df20['Travel_Angle'] = df20['SS_dir'].astype('float64') - df20['SR_dir'].astype('float64')

df20.head()

	Date	SR_time	SR_dir	Culmination	Culm_alt	SS_time	SS_dir	Daylight_Hours	Travel_Angle
0	01	2020-01-01 06:50:00	118.1	2020-01-01 11:44:00	31.3	2020-01-01 16:38:00	241.9	0 days 09:48:00	123.8
1	02	2020-01-02 06:51:00	118.0	2020-01-02 11:45:00	31.4	2020-01-02 16:39:00	242.0	0 days 09:48:00	124.0
2	03	2020-01-03 06:51:00	117.9	2020-01-03 11:45:00	31.5	2020-01-03 16:40:00	242.1	0 days 09:49:00	124.2
3	04	2020-01-04 06:51:00	117.8	2020-01-04 11:46:00	31.6	2020-01-04 16:40:00	242.3	0 days 09:49:00	124.5
4	05	2020-01-05 06:51:00	117.7	2020-01-05 11:46:00	31.7	2020-01-05 16:41:00	242.4	0 days 09:50:00	124.7

　データフレーム自体は、前回から目新しさはありませんが、データ型は以下のようにしました。

df20.dtypes

Date                     object
SR_time                datetime64[ns]
SR_dir                   float64
Culmination         datetime64[ns]
Culm_alt               float64
SS_time                datetime64[ns]
SS_dir                   float64
Daylight_Hours    timedelta64[ns]
Travel_Angle        float64
dtype: object

３．簡単に見てみる

 　jupyter notebook のセルでは、それほど長いと感じなかったのですが、ブログにコピペすると、なかなか長いですね。ということで、さっそく 2020 年の日の出入りのデータをざっくり見てみることとします。

df20.describe()

	SR_dir	Culm_alt	SS_dir	Daylight_Hours	Travel_Angle
count	366.000000	366.000000	366.000000	366	366.000000
mean	88.934699	54.686612	271.067486	0 days 12:10:55.573770491	182.132787
std	20.536070	16.482305	20.537343	0 days 01:40:01.575416347	41.073001
min	60.000000	30.900000	241.400000	0 days 09:45:00	122.800000
25%	68.725000	38.375000	250.850000	0 days 10:33:30	141.575000
50%	88.600000	55.000000	271.450000	0 days 12:12:30	182.850000
75%	109.150000	70.975000	291.325000	0 days 13:47:45	222.550000
max	118.600000	77.800000	300.000000	0 days 14:35:00	240.000000

　簡単に見るだけのつもりが、これだけでほぼすべてが語られていました。

　東京の緯度は、北緯 35°41' であることを念頭に見てみると。

① 2020 年は 366 日で、

② 一年平均すると、太陽はほぼ真東 (北から 88.9° の位置) から昇り、真西 (271°) に沈み（要するに、春分、秋分の日辺り）、

③ 春分秋分時の南中高度 (Culm_alt) は 54.6°（90° ー 東京の北緯 35°41'）で、

④ 日中の時間は 12 時間 10 分。

⑤ 南中高度が最大になるのは（おそらく夏至の時で）、77.8° (90°ー(東京の北緯ー北回帰線 23.4°))

⑥ 南中高度が最小になるのは（おそらく冬至の時で）、30.9° (90°ー(東京の北緯＋北回帰線 23.4°))

　まさに教科書通りというか、あるいは、教科書が過不足なく言い表しているというべきでしょうか。

４．簡単に図示する

describe() だけでも十分のような気もしますが、とどめにいくつか図示してみました。

　左から、一年を通じた南中高度の変化、日照時間と年中高度の関係、日照時間と太陽の移動角度の関係を図示しました。ここまできれいに精度よく記述できるからこそ、今年の節分のように、補正も精確にできるのでしょうね。

まとめ

　日の出入りに関しては、怖ろしい精度で記述が可能であることが窺えました。また、データを整えさえすれば、あっという間にまとめてくれる、python やそのライブラリにもびっくりです。

　日照時間の変化は、十分過ぎる位感じることができましたが、これに天気などのデータを加えると、季節の移り変わりを感じることとなるのでしょうか…。おいおい試してみたいと思います。

　今日は、寒さが緩んだというよりも、暖かい一日でした。2 月の第一週目で最高気温が 16 ℃位になるのは驚きです。

　本日も最後までお付き合い、ありがとうございました。

日照時間で感じる季節の移り変わり - 1 - データ収集、スクレイピングまで

節分、冬最後の日

　暦の上では最後の冬の日となりました。2 月 2 日が節分になるのは、何でも 1897 年以来 124 年ぶりだそうです。

　今日の東京は、昼間は 15 ℃位まで気温が上がり、体を動かすと、長袖一枚でも汗ばんでくる暖かな陽気でした。もっとも、朝9時近くまでは雨、そして、昼間の暖かさとは裏腹に、夕方からは信じられない位冷たい風が強まり、どんどん気温が落ちて行きました。要するに、寒冷前線の通過の影響が大きく出た一日でした。ということは、どうやら明日は寒くなりそうです。

　春遠からじと言っても、寒い日は少なからずまだありそうです。けれど、今年もそうでしたが、1 月も半ば位になってくると、日差しの力強さが増しているのを感じます。特に家から急に外出した時など、日差しの印象と気温の差に驚かされます。

　いつもは天気図と合わせて、地上付近や高層の風や気温の分布などの変化の予報を眺めていたのですが、もう少し原始的に、日照時間とか、太陽の高さがどのように変化しているのかに目を向けて、季節の変化を振り返ってみることとしました。

目次

節分、冬最後の日

１．日照時間のデータ集め

２．さあ、Beautiful Soup の出番です…

３．次は、Numpy の出番です

４．最後は、pandas です

まとめ

１．日照時間のデータ集め

　こうした情報も、いつも通り気象庁にお世話になるのかな、と思っていたところ、検索をかけてみると、どうやら、国立天文台の方が担当 (？) のようです。

eco.mtk.nao.ac.jp

　上のサイトでは、各都道府県の県庁所在地における、日の出、南中、日の入り時刻、太陽の方位、高度変化が月別に、3年分見ることができます。

　ご覧頂くと分かりますが、日々のデータの見たいところを眺めながら何となく納得…、でも良いのですが、これだけ詳細にまとまっているので、お借りして、Visualize したら面白いだろうな、ということで、久しぶりにいつもと違うデータをいじらせて頂くこととしました。

２．さあ、Beautiful Soup の出番です…

　まずは、ページのソースを拝見しましょう。

　例えば、東京都の 2021 年 2 月の『日の出入り』ページで、右クリックして「ページのソースを表示(V) Ctrl+U」を表示させてみました。

<td id="m0201"> 1</td><td> 6:41</td> <td>110.6</td> <td>11:55</td> <td> 37.3</td> <td>17:08</td> <td>249.6</td></tr> <tr><td> 2</td> <td> 6:40</td> <td>110.2</td> ……

のような感じで、時刻データ等は <td> タグで挟まれているので、ここをお借りすることとします。

import urllib.request
from bs4 import BeautifulSoup

hinodeiri = "https://eco.mtk.nao.ac.jp/koyomi/dni/2021/s1302.html"

with urllib.request.urlopen(hinodeiri) as deiri:
    html = deiri.read()
    soup = BeautifulSoup(html, "html.parser")
    srsss = soup.find_all('td')

print(srsss)

　これで、先ほど表示させたページのソースの <td> タグだけが取れていれば、データ収集ひとまず完了です。結果は以下の通りです。

[<td id="m0201"> 1</td>, <td> 6:41</td>, <td>110.6</td>, <td>11:55</td>, <td> 37.3</td>, <td>17:08</td>, <td>249.6</td>, <td> 2</td>, <td> 6:40</td>, <td>110.2</td>, <td>11:55</td>, <td> 37.6</td>, <td>17:09</td>, <td>249.9</td>, <td> 3</td>, <td> 6:40</td>, ……

良さそうです。続いて、タグ内の数値だけお借りしましょう。

A = []
for srss in srsss:
    a=srss.contents[0]
    A.append(a)

[' 1', ' 6:41', '110.6', '11:55', ' 37.3', '17:08', '249.6', ' 2', ' 6:40', '110.2', '11:55', ' 37.6', '17:09', '249.9', ' 3', ' 6:40', '109.9', '11:55', ' 37.9', '17:10', '250.3', ' 4', ' 6:39', '109.5', '11:55', ' 38.2', '17:12', '250.7', ' 5', ' 6:38', '109.1', ……

無事にリストになってくれました。

　『日の出入り』ページをから分かるように、データは、

　2月の 28 日 ( 行 ) × 「日付、日の出、方位、南中、高度、日の入り、方位」の7 列

の配列となっています。先の予定は特に決めていませんが、ひとまず配列の形を整えることとしました。

３．次は、Numpy の出番です

import numpy as np

B = np.array(A).reshape(28,7)

といっても、reshape するだけですが、格段にサイトの表に近づきました。

[[' 1' ' 6:41' '110.6' '11:55' ' 37.3' '17:08' '249.6']
 [' 2' ' 6:40' '110.2' '11:55' ' 37.6' '17:09' '249.9']
 [' 3' ' 6:40' '109.9' '11:55' ' 37.9' '17:10' '250.3']
……

　よく見ると、日の出の時刻が、時間の方が '6' のみになっています。6 が一文字なのか、空白付きなのか、結果からははっきり分かりませんが、あとで、pandas を使う時に、エラーになりそうな感じがプンプンします。

　ということで、時間も二桁表示に変えておくこととしました。

for m in range(28):
    B[m,1] = B[m,1].replace(' ', '0')

[[' 1' '06:41' '110.6' '11:55' ' 37.3' '17:08' '249.6']
 [' 2' '06:40' '110.2' '11:55' ' 37.6' '17:09' '249.9']
 [' 3' '06:40' '109.9' '11:55' ' 37.9' '17:10' '250.3']
……

　これならデータフレームで datetime を使う時にも、何の問題もなく扱えそうです。

４．最後は、pandas です

　まだ、どんな感じでデータをいじるか、何も決めていないのですが、ひとまずデータフレームにして、準備を整えます。

import pandas as pd

s = pd.DataFrame(B)
s.columns=['Date', 'SR_time', 'SR_dir', 'Culmination', 'Culm_alt', 'SS_time', 'SS_dir']
s.head()

	Date	SR_time	SR_dir	Culmination	Culm_alt	SS_time	SS_dir
0	1	06:41	110.6	11:55	37.3	17:08	249.6
1	2	06:40	110.2	11:55	37.6	17:09	249.9
2	3	06:40	109.9	11:55	37.9	17:10	250.3
3	4	06:39	109.5	11:55	38.2	17:12	250.7
4	5	06:38	109.1	11:55	38.5	17:13	251.0

まとめ

　国陸天文台の『日の出入り』ページのデータを Beautiful Soup でお借りして、Numpy、Pandas で整えるところまで到達しました。

　こちらを使わせて頂き、視覚化、時系列化等、いろいろ学ばせて頂こうと思案中です。

　今年の恵方は「南南東」だそうです。家族が買ってきた恵方巻のラベルを見て驚いたのですが、「清水寺祈祷済」と書いてありました。恵方巻は、おいしいだけでなく、実は（？）、神聖な習慣なのでしょうか…。

　本日も最後までお付き合い、どうもありがとうございました。

『伝道の書』の動詞の Lemmatize を試みる

　前回はテキストの品詞分解を試してみました。動詞には時制があり、形容詞や副詞には比較級があり、と分類すると、品詞のタグが英文の場合、36 にも上ることは、ある意味当然なのかもしれません。しかし、辞書で単語を調べる時のように、各単語をその基本形に帰着させ、基本形の使用頻度を見ることが出来れば、そのテキストの語彙の傾向が抽出できるのではないか、と期待したのですが…。

目次

１．動詞の基本形化

　１．１　基本形化の機能 WordNetLemmatizer

　１．２　Lemmatize された動詞の品詞タグの変化

　１．３　Lemmatize された動詞の中身

２．まとめと所感

１．動詞の基本形化

　基本形化することですっきり整頓出来そうなトークンとして、まず動詞が思い浮かびました。動詞の品詞タグは 6 つに分かれていますが、時制や人称、活用を統一して把握することが出来れば、品詞タグのばらつきがグッと減るのではないか、と思い、試してみたのですが…。

１．１　基本形化の機能 WordNet Lemmatizer

　公式ページによると、nltk には WordNet Lemmatizer という機能があり、これを使うことで、基本形化が出来るそうです。ただし、WordNet にない単語はそのまま返す、とあります。

www.nltk.org

　WordNet の辞書の中身がどのようなものかは分かりませんが、動詞を対象としてやってみることとします。Lemmatizer は以下の形で使います。

from nltk.stem import WordNetLemmatizer
wnl = WordNetLemmatizer()
wnl.lemmatize('基本形化したい言葉')

 　Lemmatize した動詞群の品詞タグがどう変化するのか、比べてみました。

１．２　Lemmatize された動詞の品詞タグの変化

　下図中、左は『伝道の書』をトークン化した単語のうち、動詞系とタグ付けされたものの数、右はトークン化した単語のうち、動詞系とタグ付けされたものを lemmatize し、それについて再度品詞のタグ付けを行ったものの数を表します。 

 　トークン化しただけでは、動詞過去形（VBD）と区分されるものが最多で、次に現在形（VBP）、原形（VB）、過去分詞（VBN）、3人称単数現在（VBZ）、動名詞または現在分詞（VBG）となっていました。lemmatize すれば、これらが原形もしくは現在形にかなり集約されるのではないかと考えていたのですが、その予想は大きく裏切られ、実際には大半が名詞、あるいは形容詞と品詞タグ付けされ、その上細分化されることとなりました。以下がその code です。

token = nltk.word_tokenize(ecclesia)
stop_words = nltk.corpus.stopwords.words('english')
symbol = ["'", '"', ':', ';', '.', ',', '-', '!', '?', "'s"]
clean_token = [w.lower() for w in token if w.lower() not in stop_words + symbol]
pos = nltk.pos_tag(clean_token)

wrd_VB = []
for i in range(len(pos)):
 if "VB" in pos[i][1]:
  wrd_VB.append(pos[i])

word, hinshi = zip(*wrd_VB)

lem_word = []
wnl = WordNetLemmatizer()
for s in word:
 lem_word.append(wnl.lemmatize(s, pos="v"))

lem_pos = nltk.pos_tag(lem_word)
lem_word, lem_hinshi = zip(*lem_pos)

df = pd.DataFrame({'Word':hinshi})
dg = pd.DataFrame({'Word':lem_hinshi})

fig=plt.figure(figsize=(11,5))
sns.set()
sns.set_style=("darkgrid")

plt.subplot(121)
plt.xlim(0,200)
plt.title('Tokenized Verb Tags')
sns.countplot(y="Word", data=df, palette="magma",
 order=df.Word.value_counts().iloc[:20].index)

plt.subplot(122)
plt.xlim(0,200)
plt.title('Lemmatized Verb Tags')
sns.countplot(y="Word", data=dg, palette="magma",
 order=dg.Word.value_counts().iloc[:20].index)

 　前回試みたように、トークン化した動詞の中には、少なからず違う品詞の単語も混じっていました。それらが Lemmatize の後、元の品詞と判別される可能性は大いにあるように思います。それにしても、500 程度あった動詞タグのうち、280 程度が名詞と形容詞と判別されるのは、ちょっと行き過ぎです。ということで、Lemmatize された動詞の品詞がどう判別されているのか、中身を一部見てみることとしました。

１．３　Lemmatize された動詞の中身

print(lem_pos)

('ecclesiastes', 'NNS'), ('king', 'VBG'), ('saith', 'JJ'), ('labour', 'JJ'), ('passeth', 'NN'), ('cometh', 'NN'), ('ariseth', 'VBZ'), ('hasteth', 'JJ'), ('arise', 'NN'), ('wind', 'NN'), ('accord', 'NN'), ('run', 'VB'), ('sea', 'NN'), ('come', 'VB'), ('satisfy', 'NN'), ('see', 'NN'), ('fill', 'JJ'), ('hear', 'VB'), ('do', 'VBP'), ('do', 'VB'), ('say', 'VB'), ('see', 'VB'), ('remembrance', 'JJ'), ('come', 'VB'), ('come', 'VBP'), ('give', 'JJ'), ……

　語尾が th の過去形は動詞と判別されず、名詞や形容詞になっています。また、satisfy や fill といった単語も動詞と判別されていないようです。一方、king は動名詞と判別されたようです。

　形容詞や副詞についてもやってみたのですが、トークン化後の分類から、動詞同様に細分化され、余計にわからなくなってしまいました…。

２．まとめと所感

　残念ですが、現時点では Lemmatize された動詞の品詞がどうしてこのように変化するのかについて、理由は思いつきません（泣）。Lemmatize した後に pos_tag という操作をすることに問題があるのか、あるいはそもそも古い英語のテキストを対象にしているために誤判別が起きやすいのか、それとも、WordNet という辞書の特徴をもっと理解する必要があるのか。いずれにせよ、ひとまず解決の糸口は思いつかないので、もう少し別のテキストも対象にぼちぼちやることとしました。

　数値化された予測などの正当性が 80 % ということと、言語のような質的データを定量化したものの正当性が 80 % ということには、水準は同じでも、私たちにとっては、やはり意味が異なります。言葉の意味を取り違えても、問題ない場合ももちろんありますが、時としてそれは表現しているものがまったく変わる可能性もあります。言語処理の難しさや面白さはこの点にあるんでしょうね。

　今回は単に達成できなかったことの報告となってしまっただけでした。

　最後までお付き合いありがとうございました。

『伝道の書』を品詞分解する - pos_tag によるカウント -

　先日テキストのトークン化を試してみました。内容によって大きく変わりそうですが、英文は、どういった品詞がどのような割合の下に成り立っているのでしょうか。前回に引き続き、『伝道の書』を品詞分解してみました。

目次

１．nltk の品詞取得

　１－１．品詞取得の command

　１－２．品詞タグの種類

２．『伝道の書』の品詞構成

３．品詞別高使用単語上位 20

４．まとめと課題

１．nltk の品詞取得

１－１．品詞取得の command

 　トークン化したテキストの、単語の品詞は次のように取得できます。ついでに、どんな結果になるのか、print してみました。

pos = nltk.pos_tag(clean_token)　 # clean_token : tokenized text
print (pos)

[('ecclesiastes', 'VBZ'), ('preacher', 'RB'), ('1:1', 'CD'), ('words', 'NNS'), ('preacher', 'RB'), ('son', 'NN'), ('david', 'NN'), ('king', 'VBG'), ('jerusalem', 'NN'), ('1:2', 'CD'), ('vanity', 'NN'), ('vanities', 'NNS'), ('saith', 'VBP'), ('preacher', 'JJ'),　......

　左が単語、右が品詞を表すタグのリストが得られます。

１－２．品詞タグの種類

　では、品詞タグにはどのようなものがあるのでしょうか？公式ページにもありますが、次の command でタグの一覧や説明が見られます。

nltk.help.upenn_tagset()   # All list of pos tag

nltk.help.upenn_tagset('VBP')   # Description about 'VBP' tag

www.nltk.org

　tag は全部で 36 に分けられています。一覧は、上の command、あるいは例えば以下のページをご参照頂くとして、

www.ling.upenn.edu

大まかに次のように分けてみました。

　　①名詞系（NN、NNS、NNP、NNPS）

　　②動詞系（VB、VBD、VBG、VBN、VBP、VBZ）

　　③形容詞系（JJ、JJR、JJS）

　　④副詞系（RB、RBR、RBS）

　　⑤疑問詞系（WDT、WP、WP$、WRB）

　　⑥その他（接続詞など）

果たしてかなりザックリですが、トークン化した『伝道の書』がどのような品詞比率になっているか、実際に図示してみました。

２．『伝道の書』の品詞構成

　トークン化では、大文字を小文字にして単語を統一し、カンマ等の記号は省きました。先ほどの pos_tag を使って上位 20 のタグを図示したものが下のグラフです。

　一目見て分かるように、圧倒的に名詞（NN）が多く、続いて形容詞（JJ）、次に聖書特有の、文頭につく"1:1" のような基数（CD）が続きます。その後は、副詞（RB）、動詞系（VBD、VBP、VB、VBN、VBZ、VBG）が目立ちます。先日トークン化した際には、一番多く使われていた単語は、"shall"、続いて"man" でしたが、全体としては名詞が一番多いようですね。

上の script はこんな感じです。

import pandas as pd
import nltk
from nltk import word_tokenize,sent_tokenize
import seaborn as sns
import matplotlib.pyplot as plt

with open("Ecclesiastes.txt", mode="r",encoding="ms932", errors="ignore") as file:
 ecclesia = file.read()
 
token = nltk.word_tokenize(ecclesia)
stop_words = nltk.corpus.stopwords.words('english')
symbol = ["'", '"', ':', ';', '.', ',', '-', '!', '?', "'s"]
clean_token = [w.lower() for w in token if w.lower() not in stop_words + symbol]
pos = nltk.pos_tag(clean_token)

word, hinshi = zip(*pos)　　　　　　　　# create dataframe
df = pd.DataFrame({'Word':hinshi})

fig=plt.figure(figsize=(7,7))
sns.set()
sns.set_style=("darkgrid")
sns.countplot(y="Word", data=df, palette="magma",
 order=df.Word.value_counts().iloc[:20].index)

３．品詞別高使用単語上位 20

　やや繰り返しになりますが、単語全体の使用頻度の上位ではなく、品詞別に使用率の高い単語上位 20 をまとめてみました。

　まずは名詞から。

有難い言葉が並びますが、hath は動詞ですし、wise は形容詞ですね。

　続いて形容詞です。

こちらも、sun、god、labour など名詞と思われる単語が混じっています。

　では動詞はどうでしょう。

こちらも、god、labour、king といった名詞が混じっていますね。

　最後に副詞を調べてみます。

　こちらもやはり、数単語は副詞で無さそうなものが混じっています。ちなみに、品詞別の図示は、Dataframe を作った 2 行を以下の script に変えただけです。

wrd_VBP = []
for i in range(len(pos)):
 if "RB" in pos[i][1]:　　　　　　# RB は副詞。欲しい tag をどうぞ
 　wrd_VBP.append(pos[i][0])

df = pd.DataFrame({'Word':wrd_VBP})

４．まとめと課題

　以上、『伝道の書』を品詞分解してみました。大筋はかなり納得、繰り返される言葉から、内容が垣間見られる想いがしますが（ちょっと言い過ぎですね）、どの品詞分解に関しても、違う品詞の単語が混じってしまうということは、トークン化からタグ付けの間にもう一手間必要なようですね。

　ライブラリの使い方にとどまらず、英語の構造や語源自体の理解も必要なようで、やはりなかなか奥が深そうです。

　最後までお付き合いありがとうございました。

『伝道の書』を自然言語処理する - nltk の install、トークン化、使用頻度の plot まで -

　先日の記事を書き、とにかく始めようということで、ライブラリを Install、あれこれいじっていると…、こんなところにも matplotlib が使われているんですね。

目次

１．まずは Install

２．テキストデータの準備、読み込みあれこれ

　２－１．別にデータは何でもいいと思いますが…

　２－２．さあ、読み込みです

　２－３．単語のトークン化

３．トークンのカウントと視覚化

４．トークンのデータ型

５．英単語はやはり横書きがいい

　５－１．matplotlib の barh を使おう

　５－２．seaborn だって使える

６．まとめと野望

１．まずは Install 

　『Python による AI プログラミング入門』の 10 章、自然言語処理を片手に、言われるがままライブラリを Install します。

www.oreilly.co.jp

pip install nltk

　場合によっては、ライブラリの一部はすでに Install されているかもしれません。私は Anaconda 3 / Jupyter notebook を普段使っていますが、一部はすでに Install されているようでした。

　いずれにせよInstall が完了したら、次は nltk のデータセットのダウンロードです。このデータセットは、色々な分析の判断基準となる辞書や物語のようです。Python Shell からでも、Jupyter notebook からでもダウンロード可能です。

import nltk
nltk.download()

　Window が新しく出てくるので "ALL" を選びます。回線の早さにもよると思いますが、わりと時間がかかります（1,2 分ではないです）。

　結構色々なものをダウンロードしていたので公式ページを覗いてみると、リストがありましたが、膨大でした…。必要になったら見返そうということで、今回はスルーしました。

www.nltk.org

　ちなみに、データセットには以下の 9 つのお話（book と言うようです）が入っています。

from nltk.book import *
texts()

とすると、お話の一覧が得られます。

text1: Moby Dick by Herman Melville 1851
text2: Sense and Sensibility by Jane Austen 1811
text3: The Book of Genesis
text4: Inaugural Address Corpus
text5: Chat Corpus
text6: Monty Python and the Holy Grail
text7: Wall Street Journal
text8: Personals Corpus
text9: The Man Who Was Thursday by G . K . Chesterton 1908

ついでに、

sents()

とすると、それぞれの text の書き出しが見られます。

sent1: Call me Ishmael .
sent2: The family of Dashwood had long been settled in Sussex .
sent3: In the beginning God created the heaven and the earth .
sent4: Fellow - Citizens of the Senate and of the House of Representatives :
sent5: I have a problem with people PMing me to lol JOIN
sent6: SCENE 1 : [ wind ] [ clop clop clop ] KING ARTHUR : Whoa there !
sent7: Pierre Vinken , 61 years old , will join the board as a nonexecutive director Nov. 29 .
sent8: 25 SEXY MALE , seeks attrac older single lady , for discreet encounters .
sent9: THE suburb of Saffron Park lay on the sunset side of London , as red and ragged as a cloud of sunset .

　横道にそれました。あと一つ、gensim というライブラリも Install してひとまず準備完了です。

conda(あるいはpip) install gensim

２．テキストデータの準備、読み込みあれこれ

　２－１．別にデータは何でもいいと思いますが…

　何ができるのか、全く全貌は見えていませんが、まずは、色々操作するためのテキストデータが必要そうです。データセットの text 3 に旧約聖書の『創世記』が入っているのを見て、じゃあ私は、同じく旧約聖書から、『伝道の書』と呼ばれる『コヘレトの言葉（Ecclesiastes）』をデータとしようと決めました。

　興味のあるものをインターネットから探して、text ファイル等にすればよいと思います。

　２－２．さあ読み込みです

　何となく予想はしていたのですが、指定しないと文字コードのエラー等出る可能性が高いです。今回はひとまず先に進みたかったので、以下のようにしました。テキストデータは、"ecclesia" という名前にしました。読み込んだデータの一部を print してあります。

with open("Ecclesiastes.txt", mode="r", encoding="ms932", errors="ignore") as file:
 ecclesia = file.read()
 print(ecclesia)

　当たり前ですが、普通のテキストファイルです。

Ecclesiastes
or
The Preacher
1:1 The words of the Preacher, the son of David, king in Jerusalem.
1:2 Vanity of vanities, saith the Preacher, vanity of vanities; all is vanity.
1:3 What profit hath a man of all his labour which he taketh　…

２－３．単語のトークン化

　文章の分析にあたり、単語やカンマ、ピリオドなどに分割し、分割したものをトークン、作業をトークン化というそうです。ひとまずやってみます。

import nltk
from nltk import word_tokenize,sent_tokenize　　　# import library

token = nltk.word_tokenize(ecclesia)　　　　　　　 # トークン化

stop_words = nltk.corpus.stopwords.words('english')
symbol = ["'", '"', ':', ';', '.', ',', '-', '!', '?', "'s"]
clean_token = [w.lower() for w in token if w.lower() not in stop_words + symbol]

print(clean_token)

　トークン化すると、カンマやピリオドも一つのデータとなるため、そういったものを除外し、また、文頭に来た単語は、最初の文字が大文字になるので、それを小文字にする作業をします。そうすれば、テキストデータは意味のある単語に分割されるはずです。上の結果は、次のようになりました。

['ecclesiastes', 'preacher', '1:1', 'words', 'preacher', 'son', 'david', 'king', 'jerusalem', '1:2', 'vanity', 'vanities', 'saith', 'preacher', 'vanity', 'vanities', 'vanity', '1:3', 'profit', 'hath', 'man', 'labour', 'taketh', 'sun', '1:4', 'one', 'generation', 'passeth', 'away',　 …

　確かになっているようです。

３．トークンのカウントと視覚化

　いろいろ解析はあるのでしょうが、各トークンがどのくらい現れているのか、言い換えれば、ある単語がどの位使われているかを数えるコマンドがあり、しかも図示も出来ます。詳しくは、公式サイトをご覧頂くとして、使用頻度上位 20 個を図示してみました。

clean_frequency = nltk.FreqDist(w.lower() for w in token if w.lower() not in stop_words + symbol)
clean_frequency.plot(20, cumulative=False)

　横軸が左から、使用頻度の多いトークン、縦軸がその回数です。便利です。でも、日本語や中国語は、縦書き横書きあっても良いと思うのですが、英語で縦書きはちょっとなあ、というのが否定できない印象だったりもします。

　しかし…、これは何となく（いえ、どう見ても）matplotlib ですよねー。トークンのデータ 型 を見てみることにしました。

４．トークンのデータ型

 　clean_token の type を調べてみると、予想通り、

<class 'list'>

となり、list でした。というか、「２－３．単語のトークン化」で print した clean_token は、よく見ればリストそのものですね。

　リストと分かれば、あとは単なる視覚化の問題です。こうして nltk から Visualization へと脱線して行きました。

５．英単語はやはり横書きがいい

　５－１．matplotlib の barh を使おう

　nltk.FreqDist で試した時同様、使用頻度の高いトークン上位 20 個を横棒グラフにしてみました。やっぱり、英単語は横書きが良いと思うのですが。 

　matplotlib で棒グラフに出来ると分かると、色もちょっといじりたくなってしまいました。

import collections

c = collections.Counter(clean_token)
c20 = c.most_common(20)
c20r =list(reversed(c20))

import matplotlib.pyplot as plt
import numpy as np

sns.set()
sns.set_style=("darkgrid")
fig=plt.figure(figsize=(7,7))

y_pos = range(0, 20)
clrs = ['darkolivegreen', 'limegreen', 'greenyellow', 'darkkhaki']
plt.yticks(y_pos, word, fontweight='bold', rotation=10, fontsize='14')
plt.barh(y_pos, count, color=clrs)
plt.xlabel('\nWord Count', fontweight='bold', color = 'black', fontsize='14', verticalalignment='center')
plt.ylabel('Word\n', fontweight='bold', color='black', fontsize='14', verticalalignment='center')

plt.show()

５－２．seaborn だって使える

　seaborn なら、countplot があるので、上位 20 語の指定はもっと簡単です。

import seaborn as sns
import matplotlib.pyplot as plt
import pandas as pd

df = pd.DataFrame({'Word':clean_token})

fig=plt.figure(figsize=(7,7))
sns.countplot(y="Word", data=df, palette="magma",
 order=df.Word.value_counts().iloc[:20].index)

　もっとも、FreqDist 一行で済んだ code をここまで増やすことに価値を見出すかどうかは、やっぱりお好みでしょうね…。

6．まとめと野望

　今回は、nltk を Install し、トークン化を試みました。そして、トークンのカウントと視覚化で、大きく脱線して行きました。

　英語の NLP でやることもたくさんあるようですが、まずは日本語が扱える環境を整えたいと思います。

　本日も最後までお付き合いありがとうございました。