[Python][OpenGL][pyqtgraph] 複数直方体をヌルヌル動かすだけ

複数直方体を動かすサンプルって意外と見つからなかったので書いてみた:
Cubes are moving with pyqtgraph opengl (github.com)

以下、お時間ある方向けの解説(という名の雑記)をつらつらと。

最初に出会ったのはpygameだった

Python+OpenGLで立方体/直方体動かすか！と思ってググったとある日。

私はキューブ回すおじさんのチュートリアルページに飛ばされた:
Python Programming Tutorials

このおじさんの説明だけだと「edges?? vertices?? ハァ??」と困惑すること間違いなしなので、
後日見つけたこちらのページもついでにどうぞ:
Advanced OpenGL in Python with PyGame and PyOpenGL - Stack Abuse

どうやら8個の頂点をvertices, どこ繋ぐかという二点の組み合わせをedgesと呼んでいるらしい。

このサンプルを回すとたしかにcubeが回転する。素晴らしい！

しかしこのサンプルはpygameを使っている。
pygameだとグリッドラインを引いたりxyz軸表示したりするのがことごとく自作しないといけない。
たぶん。わかりません。

ということで私は早々にpygameに別れを告げ、pyqtgraphの道を歩むのであった。

PyQtGraphは控えめに言って素晴らしい

pyqtgraph。ああ偉大なるpyqtgraphよ。

pyqtgraphの偉大さについては既に解説されているのでこちらで:
Python: PyQtGraphでリアルタイムプロット／アニメーション（サンプルコード・実装デモ付き） - Wizard Notes (wizard-notes.com)

Gridやxyz軸がそれぞれ二行追加すれば表示できるのです。偉大。

早速本題の直方体作成ですが、答えは全てこちらのanswerにあります:
pyqt5 - Plot cube using pyqtgraph in python - Stack Overflow

GLMeshItemという関数を使う都合上、先程出てきたedgesではなくfacesを定義します。
faces!? これは三点の三角形の表面という意味のようです。
6面*2 = 12個の三角形を定義すれば完了。

今回のサンプルではx,y,zは適当、length, width, heightは固定にしました。
この辺りは用途に応じて入力値で上書きしてしまえばいいかと。

update関数内でxを毎サイクル-0.2だけ動かしています。 

完・成！

おわりに

ここまで紹介してアレなんですが、リアルタイム処理ではこのサンプルはお蔵入りになりそうです。。

ros2で入力与えてヌルヌル動かすぜ！という目論見でしたが、rclpyがやけに遅くて色々試しても10~12Hz程度でしか動きませんでした。。だがこれはrclpyの性能限界であってpyqtgraphさんは悪くない！

というわけで大人しくcppで書き始めました。需要があればそのうち上げます。

おしまい。

[Python] 警視庁の事故データを地図上にプロットしてみた

 今年7月に公開された2019年の交通事故データを使い、日本地図の上にプロットしました。

https://nyoshimura.github.io/JapanTrafficAccident/JapanTrafficAccident.html

1. 使ったデータ
2. データの前処理
3. 地図作成

について書きます。

1. 使ったデータ

警視庁のオープンデータを使いました。

• 出典:警察庁ウェブサイト(https://www.npa.go.jp/rules/index.html)
• データ元: (https://www.npa.go.jp/publications/statistics/koutsuu/opendata/index_opendata.html)

日本語なのでencodingはcp932を使う必要があります。

df = pd.read_csv("./honhyo_2019.csv",encoding='cp932')

2019年以前のデータも含まれていますが、データの90%以上が2019年のものです。

2. データの前処理

元のデータだと緯度経度が度分秒表記かつ整数化のため1000掛けてあるので、通常の10進法緯度経度に変換しました。 (参考: https://www.wingfield.gr.jp/archives/2687)

def dms2dec(latlong_raw):
    var_degree = int(latlong_raw/1e7)
    var_minutes = int(latlong_raw/1e5 - var_degree*100) 
    var_seconds = (latlong_raw/1e3 - var_degree*1e4 - var_minutes*100)
    decimal_latlong = var_degree + var_minutes/60 + var_seconds/3600
    
    return decimal_latlong

3. 地図作成

foliumを使ってopenstreetmap上にプロットしました。

(参考: https://qiita.com/Kumanuron-1910/items/12ce7aa02922927de2f4)

dataframe、 色、地図、凡例名を渡してレイヤー生成する関数を作成してみました。

def addLatLong2map(df_in, color_in, folium_map, legend_name):
    lgd_txt = '<span style="color: {col};">{txt}</span>'
    FeatureGroup_in = folium.FeatureGroup(name= lgd_txt.format( txt= legend_name, col= color_in))
    for index_df in range(len(df_in)):
        folium.CircleMarker([dms2dec(df_in["地点　緯度（北緯）"][index_df]), 
                             dms2dec(df_in["地点　経度（東経）"][index_df])],
                             radius=2,color=color_in,
                             fill=True,fill_opacity=0.9).add_to(FeatureGroup_in)
    folium_map.add_child(FeatureGroup_in)

# create map
folium_map = folium.Map(location=[35.415377, 139.595271], zoom_start=11)

# for loop for lat/long
df_person = df[df["事故類型"]==1].reset_index()
df_nodeath = df_person[df_person['死者数'] == 0].reset_index()
df_death = df_person[df_person['死者数'] > 0].reset_index()

start_time = time.time()

addLatLong2map(df_nodeath,'#ffa500',folium_map, 'pedestrian accident')
print("--- %s seconds for accident loop ---" % (time.time() - start_time))

addLatLong2map(df_death,'#ff0000',folium_map, 'pedestrian death')

# turn on layer control
folium_map.add_child(folium.map.LayerControl())

print("--- %s seconds for all loop ---" % (time.time() - start_time))
    
# save map
folium_map.save('./docs/JapanTrafficAccident.html')

ソースコードはnotebookですがこちら:

https://github.com/nyoshimura/JapanTrafficAccident/blob/main/FirstAnalysis.ipynb

Github pagesをホストに使ってhtmlファイルも保存しました:

https://nyoshimura.github.io/JapanTrafficAccident/JapanTrafficAccident.html

今後はこの地図上に色々情報付加していく所存。

おしまい。

[covid19] 最近の専門家提言資料を一通り読んだらテレビで教えてくれないことばかりだった話

テレビで「新しい生活様式」ばっかりやるけど、厚生労働省のHPに上がっている提言一通り読んだら知らないことばっかり書いてあってビビったので思わずまとめてみた: https://github.com/nyoshimura/covid19-JP

• 実効再生産数は数値目標として微妙
• 政府と専門家が提言内で暗に示唆している数値目標
• 緊急事態の解除条件
• PCR検査数増やしても緊急事態解除はできない

実効再生産数は数値目標として微妙

5/1のスッキリで皆さん知りたがっていた実効再生産数 (1人の感染者が何人に移すか)。

ドイツが毎日公表しているとかで謎の信頼を得ているこの数字ですが、直近20日間については実際の再生産数より低く見積もる可能性を日本の専門家達は指摘しています。

(5/1専門家提言より抜粋)
 

この数字は感染爆発を事前に検知するために活用されるべきなのに、直近20日間は低く見積もるんですよ。微妙すぎる。

日本の専門家達はこの数字より倍化時間(何日で感染者数が倍になるか)が指標として使えると提言内で示唆しています。

政府と専門家が提言内で暗に示唆している数値目標

テレビでは中々教えてくれないですけど、専門家提言資料にはたくさん数字が出てきます:

• 人工呼吸器やECMO(人工肺)より病床数が不足 (使用率38.2%, 東京100%オーバー)
• 軽症者の療養施設に関して39都道府県で13000室が利用可能だが地域別に見ると不足
• 空床数見える化に参加している医療機関は46%, 拡充必須
• 倍化時間3日を下回るとオーバーシュート傾向, 要緊急事態宣言
• 新規感染者数が100人を超えているとコントロールできない

いくつかの数字の根拠が曖昧なため、首相が宣言できるレベルではないです。

ただ、実効再生産数ではなく倍化時間を指標として緊急事態宣言を出したことでオーバーシュートを防げた。
 陽性率ではなく新規感染者数を追うことで感染爆発リスクを回避してきた。
他国やメディアに惑わされず専門家が事実ベースで対策を講じていることは間違いないと思います。

緊急事態の解除条件

これも実は専門家提言に曖昧だが書いてある:

(5/1専門家提言より抜粋)
 

数字は適当ですが、例えば下記4つが専門家提言から読み取れる日本の指標(KPI)です。

指標	条件	(地域ごとの)緊急事態
新規感染者数	e.g. 50人以上の場合	発動
倍化時間	e.g. 3日以下の場合	発動
病床使用率	e.g. 50%以下の場合	解除
空床可視化参加率	e.g. 80%以上の場合	解除

病床使用率については政府資料よりhttps://www.stopcovid19.jp/が分かりやすい:

PCR検査数増やしても緊急事態解除はできない

専門家提言資料にはPCR検査数と陽性率についても数字が出ています。 しかし検査数を増やして陽性率が多少増減したところで、

• 新規感染者が多い限り感染経路特定が追いつかない
• 病床数が足りない限り医療崩壊は防げない

という観点から、緊急事態解除には何ら影響を及ぼさないです。本質ではないのです。 

まとめ

正直想像以上に専門家提言はよくまとまっていましたが、残念なのは経済的観点がこれまで一切考慮されていないところですね。医療と疫学の観点しか提言には含まれていないため、失業数も数値目標や出口戦略の一つとして必要かと。日本では自殺数と失業数の相関が極めて高いので:

おしまい。

日本におけるCOVID19の陽性率, 重症率, 死亡率を可視化してみた

毎日テレビで数ばかり報道されても困るので。
割合を調べようと思ったら結構誰も出していないので自分で出しました。
データは東洋経済オンラインさんの物を流用しました : https://github.com/nyoshimura/covid19

陽性率の重要性

テレビでは人数ばかり報道されていますが, 100人検査した日と1000人検査した日では陽性になる人数も10倍変わりますよね。なので割合が知りたい。

あと, 4人検査して3人陽性だったら75%!!というのは母数が少なすぎですよね。なので累計で割合を出しました。

• PCR_Positive_rate(陽性率) = 陽性の人数 / PCR検査した人数
• PCR_Severe_rate(重症率) = 重症の人数 / 陽性の人数
• PCR_Death_rate(死亡率) = 死者数 / 陽性の人数

3月は重症率がかなり大きい。これは検査数の少なさが一つの要因です。
3月上旬は検査数を抑えていたため正確な陽性率が分からず, 重症化してから検査するケースが今に比べて多かったと推察されます。

一方4月に入ってからは検査数を強化しているため, 11%まで陽性率が上昇しています。
この陽性率が日本の「現状」であり, 韓国の2.6%などと比較すると「高い」と言えます。

また死亡率は4月中旬から緩やかに上昇傾向です。
感染者の増加率が一定であるにも関わらず死亡率が上昇しているのは望ましくありません。

病床数の変動, 呼吸器の稼働率, 医療従事者の稼働率などのデータも見ないと断言はできませんが, 医療が追いついていない, 疲弊している可能性が読み取れます。

今後どうなるのか

最近だと堀江さんが分かりやすい対談をアップしてました:

新型コロナウイルスについて4月末時点での状況を専門家に聞きました(前編)https://t.co/H7yrfN9uhK pic.twitter.com/v29JLAtinK

— 堀江貴文(Takafumi Horie) (@takapon_jp) April 29, 2020

ワクチン, 免疫, 重症化を防ぐ特効薬, と多方面から対策を検討している。
それでもワーストでは2年後まで終息しないそうです。

堀江さん一派の言う「緊急事態宣言延長なんてバカバカしい」という主張はもっともだと思う一方, 医療崩壊を起こすことによる死者数もまた計り知れない。

経済的弱者を捨てるのか, 患者を捨てるのか。
意見は割れるでしょうが, 前者は政府が「救える」, 後者は誰も「救えない」。
陽性率は依然10%と高く, 死亡率は上昇傾向です。

おしまい

自動運転・電動化・カーシェアが流行るのエネルギー問題に帰着する説

すべてがエネルギー問題になる。

仕事がなくなるのは自動運転のせいではなくカーシェアのせい

本日公開されたニュースによると、GMが管理職を25%, 従業員を15%減らし5つも工場を減らすらしいです。

GMが大規模リストラ発表　「電気自動車や自動運転にフォーカスするため」https://t.co/uo974l7HO3 pic.twitter.com/4fCUVfOmhX

— ITmedia NEWS (@itmedia_news) 2018年11月27日

もちろんGMだけの問題ではありません。
自動車業界ではこれから同じことが各社で起きる可能性があります。
めちゃめちゃ職を失う人が出るわけですが、これは自動運転のせいなのでしょうか。

クルマの「完全自動運転」でなくなる職業・モノ全リスト https://t.co/YjQgJmZcNj #現代ビジネス

— 現代ビジネス (@gendai_biz) 2017年12月12日

こんな記事流行りましたね。
たしかに今回の変革は「自動運転にフォーカスするため」と書いています。
が、だったら5つも工場を閉じるのは話が飛びすぎです。

自動運転によってカーシェアが普及する
→車の台数が減るから工場を減らす
→ビジネスモデルが変わるから既存ビジネスの社員減らす

なわけです。

でもカーシェアなんてそんな普及しないっしょ、と思っている方のために次で深堀りします。

カーシェアはエネルギー問題の一つの解決策

前章で「自動運転によってカーシェアが普及する」と書きましたが、厳密には違います。
そもそも今大して自動運転存在しないのにカーシェア普及し始めてますし。

月額基本料がかからないのはいいですね

カーシェア、月会費無料が続々　長期間借りっぱなしも：日経https://t.co/QtsM0T3oyx …

>2018年の会員数は前年比21.6％増の約132万人、車両数も同19.4％増の2万9200台。自家用乗用車の6000万台超はもちろん、レンタカーの60万台超とはまだ大きな差がある。 pic.twitter.com/sUOXoomeGC

— 世界四季報（セカ報） (@4ki4) 2018年11月19日

カーシェアが普及すると、

1. ユーザーが便利: 維持費なく新しい車に必要なだけ乗れる
2. ユーザーが快適: 車の総数が減るので渋滞も減るし駐車場もガラガラになる
3. 国が嬉しい: 渋滞、排ガスによる経済的損失が減る←重要

そう、エネルギー問題解決に貢献するのです。国的にはこれ大事。

貧しい人というのは、多くを持たないもののことではない。多くを求めるものをいうのだ。―ムヒカ元ウルグアイ大統領 pic.twitter.com/AkYMFoWCBv

— Liberal Japan (@LiberalJapan) 2018年11月27日

消費を減らしゴミを減らす。

消費を減らすための一つの答えがカーシェアであり、それを加速するのが自動運転。
そしてゴミを減らすためにガソリン脱却・電動化が求められる。

この清く、正しい考えにより多くの人が失業するわけです。

仕事は増える

とはいえ、消費を減らしてゴミを減らすのは簡単ではないです。

例えばご飯。
ファミレスではメニュー通り料理が提供されますが、毎日大量の残飯が捨てられています。
客の空腹度が数値化できれば料理の量をカスタマイズできます。

例えば電気。
誰もいないオフィスの電気は自動で消せるし、昼の消費量ピーク時は電気代を上げて昼休みをシフトさせたり。

「ケンブリッジTOUモデル」が遂に結果出してる！！（ちょこっと手伝ってました）

日本初となるスマートメーターデータを活用したTOU料金の実証実験の成果発表　～AIによる数値予測通り平均9.2%の電気代削減と7.7%の電力消費量削減に成功～ https://t.co/D6XawcQEr0

— ysmrnbt (@ysmrnbt) 2018年11月27日

モノづくりよりシステムづくりが増えるだけで、仕事は増えると思いますよ。
という説を提唱する者でした。

おしまい。

[Python] いまさら聞けない「ローパスフィルタ」の基礎

フーリエとかラプラスなんて知らないよ。。
周波数領域？ふーん。なるほどー。(移動平均でいいじゃん。)
そんな方向けにローパスフィルタについてまとめてみます。

まずはフーリエ変換の雰囲気をつかむ

ローパスフィルタを知るにはまず「フーリエ変換」を(雰囲気だけでも)知っておく必要があります。
ネットで調べれば下の図のような説明は多く見かけるでしょう。
しかし物事を理解するには歴史、背景、目的みたいなものを多少でも知ることが大事だと思うのです。

フーリエ変換 = 時系列信号を周波数成分(sin波)に分解する変換

ですが、そもそもなぜこんな変換を昔の人はしたかったのか。
その起源は波動にあります。

例えば音。子供しか聞こえないモスキート音は「周波数が高い」。
あるいは光。紫外線は「波長が短い」。
ラジオなど無線にも「周波数帯」なるものがある。

この世のいろんな信号は「波」に変換され「周波数」なるもので表現されています。
この世のあらゆる信号は波に変換できる(仮)。その変換こそがフーリエ変換なわけで。


上のgifは少し前にtwitterで回ってました。
いろんな周期の波を足せばどんな信号も表現できるのです。

こうして信号をフーリエ変換すると、どの周波数がどれだけ含まれているか、というヒストグラムが出ます。

左を時間領域、右を周波数領域、と表現したりもします。
ではなぜこんな変換をするのか。答えは「ノイズ」です。

一番わかりやすいのは音。イヤホンやラジオの音質は大事ですよね。
あるいは光。UVカットサングラス、はどうやってカットしてるんでしょう？？

欲しい信号だけ通して、いらない信号はカットする。それが「フィルタ」であり、
例えば高い周波数のノイズをカットするのが「ローパスフィルタ」なわけです。

次に伝達関数の雰囲気をつかむ

このフーリエ変換、制御的にとても嬉しい点が一つあります。
それは、周波数領域は時間に依存しない表現ができること。
つまり定常状態を表現できるわけです！

例えば車のハンドルを右に切った時。
車がどれくらい遅れて、どれくらい右に動くかを一つの数式で表現できます。

これって実はめちゃめちゃすごいことで、
本当はハンドルを切る時間とか量によって時間領域の信号は毎回違うわけです。
でもそんなの関係ねぇ！という物理的なモデル構築がフーリエ様によって可能になったわけです。

一つの数式(=モデル=伝達関数)で、ハンドルから車までを表現できる。
フーリエ変換のおかげでモデルベース制御なるものが生まれたわけですね。
(雰囲気説明なので多少の誤りはお許しを)

ローパスフィルタの伝達関数

ローパスフィルタはwikipediaを見ての通り、電気系が起源です。
電気信号はとにかくノイズ(高い周波数の不要な信号)が多いので、それを除去したいと。
以下にローパスフィルタの伝達関数を示します。

この図のτ(タウ)を時定数とか呼びます。

時定数の逆数をcut-off frequency(遮断周波数)とか呼びます。

とかそんなことはどうでもよくて。

大事なのは上図の下側、等価変換です。

伝達関数を書けても実装できなければつまらない。
実装しやすくするために伝達関数を変換することも時には必要なのです。

上図を式にするとoutput = input / (1+sτ)
output * (1+sτ) = input
output = (input - output) / sτ

という具合に等価変換できます。

1/sは時間領域では積分に相当するので、これで実装はかなり簡単になりました。

1/sが積分？なぜ？？という方はぜひラプラス変換をもっと勉強してみましょう。

そんなわけでローパスフィルタをpythonで実装すると...

一行で完結。

output = (input - output) / sτ

cut_f = 1/τなので cut_f * (input - output)を積分しています。

積分ってそうやって実装するの？という方は下図を見て納得しましょう。

(もちろん他にも実装方法はあります)

おわりに

知人に「ローパスとは」って聞かれて参考リンクでも送ろうと思ってググったら、けっこう敷居の高い記事が多いんですね。社会人になってから知らない分野を勉強するのってやっぱ大変なんだなぁ。自分もちゃんと手計算で等価変換とかしたのは院生になってからだし。ちなみにハイパスはsτ/(1+sτ)でローパスと足すと1になるんですよ。

おしまい。

Pythonで緯度経度をXY座標に変換する方法はいくつかあるらしい

異なる手法で試して値が一致しなくて困惑したのでまとめておきます。

緯度経度をXY座標に変換する

このテーマ自体は調べたら下記二つ分かりやすい記事がありました。

• Pythonで緯度経度を平面直角座標(XY座標)に変換する
• pyproj

しかし以下のStack Exchangeを試して混乱の時間が始まったのです...
Convert GPS coordinates to Web Mercator EPSG: 3857 using python/pyproj

(メルカトル法で計算した値とpyprojで出てくる値が違うぞ...)

ただ結論から言うと、違ってOKのようですね。
まずはメルカトル法から見ていきます。

Spherical Pseudo-Mercator projection (長い)

Open Street Mapのwikiにpythonコード載ってました。すばらし。
ただこれだと不完全で、上記のStack ExchangeやBokehブログでは地球の半径を基にスケーリングしてます。

何をやっているかというと、
地球を完全な球と仮定して、地球を風呂敷のように広げてX, Yに変換しているイメージです。
これで世界中どこでもX, Yが得られますが、とても簡単で大雑把な仮定なので精度がイマイチというウワサ。

EPSG-based projection

一方最初の記事のEPSGコードに基づく変換は、地球を細かく区切ってエリアごとに風呂敷を広げるイメージ。

そのため、

• 精度は前のやつより良い
• エリアごとにEPSGコードが違う
• 始点(X=Y=0)の位置もエリアごとに違う

というわけでした。そりゃ値違うわ。

さいごに

pythonのコードはGithubで公開してます。
(わざわざレポジトリ作るまでもなかったか。。。)
おしまい。