PyROOT を使う時の注意点 (書きかけ)

Python ROOT

PyROOT (Python + ROOT) を使う時の注意点をいくつか。適宜更新予定です。公式の PyROOT の manual が参考になります。

1. Memory の開放

1.1. C++ 内部で new したものの扱い

Python はメモリ管理を気にしなくて良いと言われますが、気をつけなくてはいけない場合も存在します。

#include "TGraph.h"

TGraph* newgraph()
{
  return new TGraph(1000);
}

import ROOT
import resource

ROOT.gROOT.ProcessLine('.L mem.C+')

for i in range(10000):
    graph = ROOT.newgraph() # graph という変数が新しい TGraph で上書きするだけ
    if i%1000 == 0:
        print i, resource.getrusage(resource.RUSAGE_SELF)[2]

for i in range(10000):
    graph = ROOT.newgraph()
    del graph # 試しに del を実行してみる
    if i%1000 == 0:
        print i, resource.getrusage(resource.RUSAGE_SELF)[2]

for i in range(10000):
    graph = ROOT.newgraph()
    graph.IsA().Destructor(graph) # 明示的に TObject の destructor を呼び出す
    if i%1000 == 0:
        print i, resource.getrusage(resource.RUSAGE_SELF)[2]

上記のような mem.C と mem.py を用意して、mem.py を実行してみましょう。どのような場合に new された TGraph が解放されるかが分かります。

$ python mem.py
0 42803200
1000 59510784
2000 76222464
3000 92925952
4000 109633536
5000 126341120
6000 143048704
7000 159764480
8000 176472064
9000 193183744
0 209895424
1000 226611200
2000 243314688
3000 260022272
4000 276729856
5000 293437440
6000 310153216
7000 326860800
8000 343564288
9000 360280064
0 378089472
1000 378109952
2000 378109952
3000 378109952
4000 378109952
5000 378109952
6000 378109952
7000 378109952
8000 378126336
9000 378126336

3 つ目の方法、明示的に TGraph::~TGraph() を呼び出した場合以外は、memory leak することが分かります。これは考えると当たり前なのですが、うっかり「Python は memory 管理しなくて大丈夫」と思って作業すると大変なことになります。

Python 側では new された TGraph が C++ 内部で使われ続けているのか判断できないため、Python が自動的には delete をしてくれません。

参照 http://wlav.web.cern.ch/wlav/pyroot/memory.html

2. 参照

例 1. TGraph::GetPoint(Int_t i, Double_t& x, Double_t& y)
import ROOT
import array

graph = ROOT.TGraph()
graph.SetPoint(0, 10, 20)

x = array.array('d', [0])
y = array.array('d', [0])

graph.GetPoint(0, x, y)
print x[0], y[0] # 10.0 20.0

3. ポインタ

例 1. TGraph.GetX()

TGraph::GetX() では、member 変数である Double_t* fX という配列の pointer を返します。PyROOT の中ではこれを ROOT.PyDoubleBuffer という type に変換しますが、中身は C の pointer のようなものです。

そのため、以下の例のように配列内の要素に直接演算を行ったり、配列の要素を取り出すことができます。また、当然 C の配列ですので、segmentation fault を引き起こす可能性があります。

import ROOT

graph = ROOT.TGraph()
graph.SetPoint(0, 10, 20)

print graph.GetX()[0] # -> 10.0

graph.GetX()[0] += 30

print graph.GetX()[0] # -> 40.0

graph.GetX()[10000000] # -> seg fault

4. 関数呼び出し速度の向上

PyROOT は CINT や ACliC で実行する C++ な ROOT script に比べると処理速度が遅くなる場合があります。そのひとつに、関数呼び出し時の overhead の存在があります。

これは PyROOT の問題と言うより、Python の仕様なのですが、次のように code を工夫することで、若干ですが処理速度が上がるはずです。

hist = ROOT.TH2D('hist', '', 100, -5, 5, 100, -5, 5)
x = ROOT.gRandom.Gaus() # 別に乱数じゃなくても良いが
y = ROOT.gRandom.Gaus()

for i in range(100000):
    hist.Fill(x, y) # これが単純なやり方

hist = ROOT.TH2D('hist', '', 100, -5, 5, 100, -5, 5)
x = ROOT.gRandom.Gaus()
y = ROOT.gRandom.Gaus()
fill = hist.Fill # 別名で呼べるようにする

for i in xrange(1000000):
    fill(x, y) # hist.Fill を探しに行かない分、速度が改善する

time で測定すると、10% くらい早くなったのが分かります。

python tmp.py  15.17s user 5.42s system 109% cpu 18.829 total
python tmp.py  13.86s user 4.56s system 108% cpu 16.914 total

Hilo に住んでた人の偏った Hawaii 島 (Big Island) 情報

観光

某所から Hawaii 島 (Big Island) の情報についての問い合わせがあり、返信していたら予想以上に長文になったので、もったいないので個人的なやり取りのメール本文を一部削除して転載。知らない誰かの役に立つこともあるでしょう。

        • -

僕が Big Island に住んでいた時は、遊ぶ余裕がなくてあんまり観光していないのですが、知っている範囲のみで。西側はほとんど知りません。書き出したら楽しくなってきて、えらい長文になってしまいました。Maui にも行かれるようでしたら、追加情報をお送りします。何かご質問あれば、ご遠慮なくお尋ね下さい。

1. Hilo

僕の住んでいたのは Hilo (ヒロ) という街で、Hilo 空港 (ITO) から車で 10 分くらいの場所です。SFO から直行便もあります。宿はそこそこありますが、恐らくこの街に泊まる人は、Mauna Kea に行きたい人 (観光 or 研究)、Univ. Hawaii (UH) に用事のある人、島を一周する時に休憩として 1 泊する人、diving する人、などだと思います。

Hilo の街自体に面白いものがあるか? と問われれば、特にありません。雨も多く、また街も結構寂れてるので、僕は住んでいてあまり楽しくありませんでした。雨が嫌いなので。

Hilo の観光 spot はあまり無いのですが、挙げるとすれば

Big Island Candy

クッキー屋さん。観光客が多い。味は結構美味しいほうですが、包装代がかさむのか、少し値段が高いです。何かお土産を Hilo で買うならここ。弁当箱くらいのしっかりした作りの箱に入った詰め合わせがお薦め。柄は時期によって変わるので、気に入った柄の箱であれば、食べ終わった後にも使えます。
http://g.co/maps/6d5hc

Bear's Coffee

Cafe + Bagel 屋さん。息抜きのときは、僕はここでいつも salmon bagel を食べていました。半年前くらいに San Fransisco で食べた salmon bagel より美味しかった記憶が。地球の歩き方に載っているので、日本人客が多いらしい。
http://g.co/maps/zscdw

Hilo Bay Hostel

Hilo には高級ホテルや motel はなく、そこらへんの宿も別に面白みはないのですが、ここの hostel ならのんびりとした雰囲気です。○○さんご夫妻の好みによってはお薦め。うちの夫婦は気に入りました。ただし、shower は共同。ゴキブリもいました。土曜か日曜かにすぐ側で farmer's market をやるので、そこで南国っぽい果物を買ってきて、hostel の台所で切って食べたり。
http://g.co/maps/59j8q

Harper Car and Truck Rental

観光地ではなく、レンタカー屋です。Mauna Kea にご自分で運転されて登るときには、ここで四駆を借りるのがお薦めです。なぜかというと、他のレンタカー屋だと Mauna Kea や、そこに至る Saddle Road (Route 200) を走行する場合に保険が効かないからです。事故多発地帯なので (天文関係者も必ずやる)、保険の効くレンタカーじゃないと危険です。
http://g.co/maps/886rp

Coconut Island Park

Hilo は溶岩っぽい海岸であまり綺麗な海の印象ではないのですが、Coconut Island は唯一、海沿いで落ち着ける場所と言うか、まあ、何にもない小さな島なのですが。ものっすごい狭い beach もあります。運が良ければウミガメに会えます。遭遇率 50% くらい。
http://g.co/maps/rzk26

King Kamehameha

いわゆるカメハメハ大王の像です。Honolulu、Hilo に加えて、Las Vegas にもありますね。全米で三箇所だそうです。まあ、この写真で見る以上でも以下でもないです。
http://g.co/maps/vqyzv

Ken's House of Pancakes

ほぼ唯一、夜遅くも Hilo でやってる店です。変な時間に Hilo に戻ってくると、もうこれしかありません。味はまあアメリカのそれです。でも、日本人には観光名所みたいです。キムタクも来たらしいです。「Yokozuna」サイズのものを注文すると、ドラを鳴らしてくれます。これを聞くためだけに、僕は不必要に「Yokozuna」サイズを注文します。
http://g.co/maps/z28k4

Pizza Hawaii & Deli

ここの Seafood Pizza が好きでした。
http://g.co/maps/wafn8

Imiloa Astronomy Center of Hawaii

最近できた天文系の博物館です。僕は行ったことがありません。UH の裏ですが、UH 自体は Stanford に比べて、観光しても面白くないです。
http://g.co/maps/6du38

2. Hilo 周辺

Mauna Kea + すばる望遠鏡

西側からも行けますが、Hilo を拠点にしたほうが近いです。すばる望遠鏡は一般向けの施設公開をしていますので、国立天文台のページから申し込みできます。業者を手配して連れて行って貰うか、自力での運転になります。夜間に運転する場合、車のヘッドライトを消すなどして運転する必要があると思いますので (観測用の光害防止)、初心者には危険だと思います。天体観測ツアーなど企画している業者もあるので、夜間に登頂する場合はツアーがお薦めです。
http://www.naoj.org/Information/Tour/Summit/j_index.html

僕自身は、国立天文台の教授の方の運転で Kea に二度登ったことがありますが、自力の運転は未経験です。Saddle Road は起伏とカーブの激しい道で、しかも地元民が物凄い飛ばすので危険です。ただ、最近は再舗装などがされて、僕のいた頃より真直ぐで平らになっていると思います。島の一周ではなく横断をする場合は、この道を通る必要があります。

4,000 m 超に車で簡単に行けるのは Big Island くらいしかないので、お薦めします。ただ、短時間で登るために体が慣れず高山病が出やすいです。

Waipio Valley

Hilo 周辺で僕の一番のお薦めです。普通のレンタカーで行くと、渓谷を上の展望台から眺めて「フーン」で終わるのですが、ここで Harper のレンタカーが活躍します。
http://g.co/maps/7cakk

Google Maps だと道が途中までしかありませんが、http://g.co/maps/yv2em ←こんな風に四駆以外は禁止されています。この道を下ってからが Waipio Valley の本番で、短距離ですが楽しいドライブになります。あれ？道が行き止まりになった、と思ったら、そのまま目の前の川をズンズンと進んで行ったりします。どこまで進んで良いのかも謎なのですが、川を渡ると道の続きを発見したりします。

↓こんな感じ。
http://carphotos.cardomain.com/ride_images/3/2367/4081/30917040008_large.jpg

周辺は熱帯雨林なので、Yosemite や Redwood とは違った森が楽しめます。

Mauna Loa

僕が宇宙線用の望遠鏡を建設していた場所です。観光客はほとんど来ません。ここまでの道も、一応舗装はされていますがデコボコです。当然 Harper をお薦めします。面白いか面白くないかでいうとあんまり面白くないですが、Mauna Kea に登る道に比べて溶岩がいっぱいあります。交通量もとても少ないので、気兼ねなく途中で車を停められます。Mauna Kea を正面から眺めるには最適。

↓こんな上まで登らなくても、ずっと同じような景色です。写真に見える変な建築物の中に望遠鏡が収まっています。後ろを振り返ると Mauna Kea が見えます。ここで標高 3,300 m です。RV が一台停まっていますが、この中で寝泊まりをしていました。Google Car がここまで来ているのが凄いですね。
http://g.co/maps/e7qz4

観光客で Mauna Loa に登る人は滅多にいません。行けばかなりのマニアです。

3. その他

Kilauea

僕が行った時は溶岩の流れる向きが変わっていて、最も溶岩に近い場所まで歩いても、見ることができませんでした。事前に溶岩情報を調べて、どこに行けば良いか計画を立てたほうが良いでしょう。

Hawaii 島の最南端 (South Point)

透き通った紺色の海が非常に綺麗な断崖です。飛び込んでいる若者がたくさん。一周する道から少し外れますが、お薦めです。
http://g.co/maps/ktxj9

Green Sand Beach

緑の砂と言われていますが、そんなに緑じゃなかったです。South Point 周辺から徒歩でも行けますが、結構遠い (片道 40 分くらい?)。舗装されていないので、レンタカーだと、多分 Harper の保険じゃないと駄目じゃないかと思います。時間に余裕があれば。
http://g.co/maps/hfskf

Black Sand Beach

僕が行った時は日が暮れて何も見えませんでした。

西側

栄えてるそうですね…。ほとんど何も知りません。

4. 事前学習の映画

Picture Bride

日系移民の映画です。工藤夕貴主演。
http://en.wikipedia.org/wiki/Picture_Bride_(film)

ホノカアボーイ

Waipio Valley の手前の町、Honokaa を舞台にした映画です。吉田カバンの息子の実話か何かが原作。Mauna Loa の溶岩っぽい道が出てきます。
http://ja.wikipedia.org/wiki/%E3%83%9B%E3%83%8E%E3%82%AB%E3%82%A2%E3%83%9C%E3%83%BC%E3%82%A4

フラガール

これは御存知かもしれませんが、一応。まあ、福島の話なので、関係ないと言えば関係ないです。
http://ja.wikipedia.org/wiki/%E3%83%95%E3%83%A9%E3%82%AC%E3%83%BC%E3%83%AB

プレゼンに使う物理の記号を Unicode で書こう

Mac 研究

Keynote や PowerPoint でプレゼン資料を作る時に、物理の記号を書こうとして困るときがあります。h bar などは Unicode に存在するので問題ないですが、例えば演算子を書きたい、反物質を書きたい、平方根 √ の上の棒もちゃんと書きたい、そういう需要が僕にはあります。

もちろん、LaTeXiT で PDF を挿入すればいくらでも好きな数式を貼り付けられますが、文中にちょっとだけ数式を入れたくなると難儀します。

Unicode には合成用の特殊記号が用意されています。最近 Twitter などで見かける AA (絵文字) でも、これを使った複雑なものが存在します。解説はこことか。

さて、いくつかの簡単な式を作ってみます。


▲ 図 1. Keynote + Unicode で作った物理の式。右側は LaTeXiT の出力。

[改訂第7版]LaTeX2ε美文書作成入門

[改訂第7版]LaTeX2ε美文書作成入門

作りかたの解説

反陽子

まず、普通に「p」と入力します。

次に「ことえり」の文字ビューア (Character Viewer) を開き、左側の pane から Unicode を選びます。そうすると、中央上の段のどこかに Unicode で 00000300 の箇所 (Combining Diacritical Marks、合成用発音記号かな?) が出てくるので、それを選択します。


▲ 図 2. ことえりの文字ビューア (英語環境での表示)

0305 に相当する場所に上付きの棒があるので、これを先ほどの「p」に続けて入力します。記号を double click すれば良いです。そうすると、p̅ のように表示されます。この合成用の記号は、1 つ前の文字を修飾してくれるため、このように上に棒をつけたり、色々とできます。

Lucida Grande だと italic が用意されていないので、例えば Myriad Pro*1 で italic にすると、もう少しそれっぽくなります。

演算子

演算子は、p̅ と同様の手法で文字の上にハットを乗せれば書くことができます。Unicode で 0302 を「p」の直後に入力すれば (図 2 で上棒の 3 つ左)、p̂ のように書けます。Myriad Pro だとこのハットが存在しなかったため、図 1 では Myriad Pro での例を載せていませんが、ハットだけ Lucida Grande にするような合わせ技も使うことができます。演算子はさすがに italic にしたほうがそれっぽいでしょう。

平方根

まずは「るーと」を変換して「√」を入力します。そのあとに、p̅ などと同様に上線を書けばいいわけですが、「√」の右端と綺麗に繋がってくれる保証はありません。Font によって、見栄えをましにするために微調整が必要です。

Lucida Grande の 「√s」の場合は、「√」+「s」+「上線」+「半角スペース」+「上線」を入力しています。ちょっとだけ「√」の右上と横棒の左端がずれてしまっています。

Myriad Pro の場合は、「√」+「半角スペース」+「s」+「上線」+「半角スペース」+「上線」という順番で入力しています。これだと、たまたまですが、とても綺麗にルートが描けています。今の職場の同僚が「√s」を使いまくる方々なので、横棒がなくてずっと落ち着かなかったんですが、この方法が広まるとありがたい。

上付き、下付き

Neutrino と陽子の衝突の式のように上付き・下付きがある場合は、Keynote や PowerPoint の機能で普通にやっちゃいましょう。

*1:Myriad Pro は無料で入手可能です。http://d.hatena.ne.jp/oxon/20100808/1281244820

Lion の GCC ではまった

Mac

OS X を Snow Leopard から Lion に移行したところ、WCSLIB の build でこけました。

$ cd wcslib-4.9/
$ ./configure --prefix=/usr/local
$ make -j 4
(略)
gcc -std=gnu99 -DHAVE_CONFIG_H -I. -I.. -g -O2 -c wcserr.c
wcserr.c:160: internal compiler error: Segmentation fault: 11
Please submit a full bug report,
with preprocessed source if appropriate.
See <URL:http://developer.apple.com/bugreporter> for instructions.
make[2]: *** [libwcs-4.9.a(wcserr.o)] Error 1
(略)

普通に configure/make の流れで build しようとすると、GCC が segmentation fault を起こして make できません。

Lion + Xcode 4 になってから GCC が LLVM のもの (/usr/llvm-gcc-4.2/bin/llvm-gcc-4.2 が実体) になったため、GCC と異なる挙動をする、もしくはまだ bug が潜んでいるということのようです。LLVM と GCC については、「Xcode 4 でデフォルトになった LLVM って何?」が簡潔で分かりやすいです。

$ gcc --version
i686-apple-darwin11-llvm-gcc-4.2 (GCC) 4.2.1 (Based on Apple Inc. build 5658) (LLVM build 2336.1.00)
Copyright (C) 2007 Free Software Foundation, Inc.
This is free software; see the source for copying conditions.  There is NO
warranty; not even for MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.

さて、LLVM ではない「普通の」GCC も Xcode 3 を入れると /usr/bin/gcc-4.2 として入るようです。こいつを使えば先の WCSLIB の問題は解決するようですが、自分の環境では Xcode 3 を入れていないので、代わりに clang を使います。

$ cd wcslib-4.9/
$ CC=clang ./configure --prefix=/usr/local
$ make -j 4

これで通ります。もちろん、clang も GCC とは異なるので、今回 WCSLIB ではたまたま上手く通っただけかもしれません。

ちなみに、llvm-gcc が seg fault を起こす、WCSLIB の wcserr.c の該当箇所は以下の通りです。これだけ見ただけでは正常な C の code なので、何が原因かはよく分かりません。

static int wcserr_enabled = 0;

int wcserr_enable(int enable)
{
  return wcserr_enabled = (enable ? 1 : 0);
}

twilog の統計情報から follower 数の変化を数値で取り出す方法

原発

ある Twitter ユーザが Twilog に登録している場合、その人の follower 数の推移などを見ることができます。例えば僕の場合、@AkiraOkumura だと、ここに行けば↓こういうグラフが見られます。この blog に幾つか原発関連の記事を書いたりしたことで、震災以降に大きく follower 数が増えました。


▲ @AkiraOkumura の follower 数の推移。3.11 以降に増加し、一旦落ち着いて、なぜか忘れたけど再度増加。

一人だけのグラフを眺めているだけなら、twilog の生成した図を開けば済みますが、複数の人と比較したりするには、やはり元の数値を落としてきて自分で作図したいのが理系脳。さて、どうやるか。

1. まずは Twilog の統計情報を見る

最初に、http://twilog.org/AkiraOkumura/stats/7-3 のような統計情報を開かなくてはいけません。このようなページを開くことで、Twilog の server 側に統計情報の書かれた file が upload されます。

(このページに行かないと、上記の図も出ない可能性が大。)

2. akiraokumura-7-3.txt を開く

すると、http://twilog.org/chart-dir/a/akiraokumura-7-3.txt という text file が生成されます。この中に必要な情報が詰まっています。ここで注意が必要なのは、必ず http://twilog.org/AkiraOkumura/stats/7-3 を先に開かなくては駄目、ということ。

また、@AkiraOkumura の場合は chart-dir の次の階層が /a になり、また例えば @hayano の場合は /h になります。

3. 適当なソフトで作図する

↓中身は日付と人数の羅列です。

&title= ,30,&
&x_label_style=10,#000000,2,16&
&x_axis_steps=16&
&y_ticks=5,10,9&
&line=2,#ff7f7f&
&values=182,182,183,183,184,183,185,184,183,184,186,187,186,187,187,189,189,191,192,191,190,195,197,197,197,196,196,195,197,196,196,195,195,198,197,200,200,200,201,201,202,202,202,201,201,203,203,203,202,204,204,202,203,203,204,205,205,208,207,207,207,205,206,207,206,206,206,207,207,208,209,210,209,209,210,209,209,208,208,207,210,209,207,207,207,209,208,209,212,213,213,217,223,224,226,228,227,227,228,228,227,227,226,225,224,226,225,226,225,225,228,228,231,231,233,232,232,233,234,234,234,237,238,236,235,248,248,253,255,253,254,252,252,255,255,255,256,257,257,258,259,259,260,263,265,264,263,259,269,269,270,270,270,270,269,274,275,279,282,281,283,284,289,290,292,290,288,288,288,285,286,285,284,285,285,317,317,318,423,450,502,511,620,645,670,690,736,811,832,960,993,1006,1015,1018,1017,1022,1067,1108,1113,1154,1188,1200,1207,1225,1499,1757,1816,1826,1841,1863,1943,1996,2003,2042,2047,2074,2098,2135,2253,2361,2377,2436,2482,2502,2520,2536,2554,2564,2587,2599,2588,2589,2633,2640,2653,2747,2803,2815,2824,2868,2869,2875,2876,2871,2877,2881,2872,2881,2881,2906,2909,2910,2906,2906,2919,2931,2931,2944,2945,2948,2951,2958,3071,3060,3059,3065,3061,3077,3104,3135,3145,3150,3192,3202,3214,3214,3230,3234,3246,3262,3262,3267,3306,3339,3394,3398,3397,3394,3397,3443,3444,3443,3440,3442,3439,3437,3436,3438,3434,3428,3421,3422,3419,3412,3413,3428,3433,3438,3434,3440,3444,3451,3463,3464,3465,3504,3520,3530,3525,3525,3525,3525,3521,3517,3518,3520,3520,3522,3519,3516,3509,3512,3517,3520,3518,3515,3512,3523,3524,3519,3517,3526,3529,3527,3526,3526,3527,3530,3527,3525,3525,3525,3527,3532,3537,3534,3540,3543,3539,3528,3524,3525,3530,3526,3527,3525,3521,3518,3505,3505,3506,3502,3506,3505,3514,3601,3627,3635,3641,3640,3644,3650,3654,3652,3649,3672,3695,3699,3703,3700,3698,3699,3698,3693,3696,3697,3732,3758,3794,3837,3861,3887,3905,3906,3920,3922,3924,3923,3923,3919,3920,3915,3911,3907,3908,3908,3908,3908,3906,3906,3906,3909,3939,3978,3985,3989,3994,4000,3998,3993,4011,4014,4014,4016,4010,4012,4012,4010,4001,3999,3997,3995,3990,3998,3995,3995,3993,3996,4001,4001,4013,4021,4038,4045,4042,4045,4045,4039,4039,4035,4035&
&x_labels=08/14,08/15,08/16,08/17,08/18,08/19,08/20,08/21,08/30,09/01,09/07,09/08,09/09,09/10,09/11,09/12,09/13,09/14,09/15,09/16,09/17,09/18,09/19,09/20,09/21,09/22,09/24,09/25,09/26,09/27,09/28,09/29,09/30,10/01,10/02,10/03,10/05,10/06,10/07,10/08,10/09,10/10,10/11,10/12,10/13,10/14,10/15,10/16,10/17,10/18,10/19,10/21,10/22,10/23,10/24,10/25,10/26,10/27,10/28,10/29,10/30,10/31,11/01,11/02,11/03,11/04,11/05,11/06,11/07,11/08,11/09,11/10,11/11,11/12,11/13,11/14,11/15,11/16,11/17,11/18,11/19,11/24,11/25,11/26,11/29,11/30,12/01,12/02,12/03,12/04,12/05,12/06,12/07,12/08,12/09,12/10,12/11,12/12,12/14,12/15,12/16,12/17,12/18,12/19,12/20,12/21,12/22,12/23,12/24,12/25,01/03,01/04,01/05,01/06,01/07,01/08,01/10,01/11,01/12,01/13,01/14,01/15,01/16,01/17,01/18,01/19,01/20,01/21,01/22,01/23,01/24,01/25,01/26,01/27,01/28,01/29,01/30,01/31,02/01,02/02,02/03,02/04,02/05,02/06,02/07,02/08,02/09,02/10,02/11,02/12,02/13,02/14,02/15,02/16,02/17,02/18,02/19,02/20,02/21,02/22,02/23,02/24,02/25,02/26,03/01,03/02,03/03,03/04,03/05,03/06,03/07,03/08,03/09,03/11,03/12,03/13,03/14,03/15,03/16,03/17,03/18,03/19,03/20,03/21,03/22,03/23,03/24,03/25,03/26,03/27,03/28,03/29,03/30,03/31,04/01,04/02,04/03,04/04,04/05,04/06,04/07,04/08,04/09,04/10,04/11,04/12,04/13,04/14,04/15,04/16,04/17,04/18,04/19,04/20,04/21,04/22,04/23,04/24,04/25,04/26,04/27,04/28,04/29,04/30,05/02,05/03,05/04,05/05,05/06,05/07,05/08,05/09,05/10,05/11,05/12,05/13,05/14,05/15,05/16,05/17,05/18,05/19,05/21,05/22,05/23,05/24,05/25,05/26,05/27,05/28,05/29,05/30,05/31,06/01,06/02,06/03,06/04,06/05,06/06,06/07,06/08,06/09,06/10,06/11,06/12,06/13,06/14,06/15,06/16,06/17,06/18,06/19,06/20,06/21,06/22,06/23,06/24,06/25,06/26,06/27,06/28,06/29,06/30,07/01,07/02,07/05,07/06,07/07,07/08,07/09,07/10,07/11,07/12,07/13,07/14,07/15,07/16,07/17,07/18,07/19,07/20,07/21,07/22,07/23,07/24,07/25,07/26,07/27,07/28,07/29,07/30,07/31,08/01,08/02,08/03,08/04,08/05,08/06,08/07,08/08,08/09,08/10,08/15,08/20,08/21,08/22,08/23,08/24,08/25,08/26,08/27,08/28,08/29,08/30,08/31,09/02,09/03,09/04,09/05,09/06,09/07,09/08,09/09,09/10,09/11,09/12,09/13,09/14,09/15,09/16,09/17,09/18,09/19,09/21,09/24,09/25,09/27,10/02,10/03,10/05,10/06,10/07,10/09,10/10,10/11,10/12,10/13,10/14,10/16,10/17,10/18,10/19,10/20,10/21,10/22,10/23,10/24,10/25,10/26,10/27,10/28,10/29,10/30,10/31,11/01,11/02,11/03,11/04,11/05,11/06,11/07,11/08,11/09,11/10,11/11,11/12,11/13,11/14,11/15,11/16,11/17,11/18,11/19,11/20,11/21,11/22,11/23,11/24,11/25,11/26,11/27,11/28,11/29,11/30,12/01,12/02,12/03,12/04,12/05,12/06,12/07,12/08,12/09,12/10,12/11,12/12,12/13,12/14,12/15,12/16,12/17,12/18,12/19,12/20,12/21,12/22,12/23,12/24,12/25,12/26,12/27,12/28,12/29,12/30,12/31,01/01,01/02,01/03,01/04,01/05,01/06,01/07,01/08,01/09,01/10,01/11,01/12,01/13,01/14,01/15,01/16&
&y_min=0&
&y_max=4500&
&bg_colour=#ffffff&
&x_axis_colour=#808080&
&x_grid_colour=#dddddd&
&y_axis_colour=#808080&
&y_grid_colour=#dddddd&

適当なソフトで整形すれば、複数人の比較ができます。


▲ 原発関連で有名な方々をいくつか。

みんな震災以降に急上昇していますね。@hayano さんの場合、突発的に上昇して、その後ゆるやかに減少。@hayakawayukio さんの場合、訓告事件のときにドカンと増えています。

施設名から経緯度を一括自動取得する方法

原発 Python

今さら感のある話題ですが、経緯度の公開されていない放射線量の自治体公開データから、住所や経緯度を一括で取得するやり方の説明をします。使用言語は Python です。Python 以外の言語でも、同じやり方でできるはず。1600 箇所の可視化、4400 箇所の可視化のときに使った方法です。

福島県が経緯度を公開してくれれば良かったのですが、例えばこの PDFのように、施設名しか公開してくれませんでした。そのときの文句は「4 月 19 日に福島県災害対策本部原子力班と電話した内容」に書きました。

1. 「iタウンページ」で住所を取得

NTT の提供するiタウンページを使用すると、施設名から住所を取得できます。例えば福島市立御山小学校の場合だと、こんな感じで検索すれば住所が表示されます。同名の小学校が存在する可能性もあるので、「福島県」や「福島市立」などを入れておくと正確に取得できます。

2. Google Maps で経緯度を取得

次に、取得した住所から経緯度に直すには Google Maps を使います。御山小学校の場合は住所が「福島県福島市御山字長滝１-１」ですので、これを引数に渡してこんな感じで情報を取ります。そうすると、うまく見つかった場合には「140.4585484, 37.7786954」という数値が取得できます。検索が失敗する場合もあります。

3. 使った Python Script

まずは元の PDF を Acrobat なりで開いた後に少し手を加えて、次のような CSV を作ります。

http://www.pref.fukushima.jp/j/schoolmonitamatome.pdf
,,,,,※測定値（μSv/h）
市町村,No.,区分,調査地点,調査月日,1m高さ,1cm高さ
福島市,1,公立小学校,福島市立御山小学校,2011/04/05,4.9,6.3
福島市,2,公立小学校,福島市立三河台小学校,2011/04/05,3.3,4.1
福島市,3,公立小学校,福島市立松川小学校,2011/04/05,2.3,2.7
福島市,4,公立小学校,福島市立下川崎小学校,2011/04/05,4.2,4.9
福島市,5,公立小学校,福島市立水原小学校,2011/04/05,2.1,2.7
福島市,6,公立小学校,福島市立金谷川小学校,2011/04/05,3.1,3.8

これを、次の Python script で読み込みます*1。

# -*- coding: utf-8 -*-

import codecs
import urllib2
import os

f = codecs.open('school.csv', 'r', 'utf')

for line in f.readlines()[3:]:
    name = line.split(',')[3]
    town = line.split(',')[0]

    if name.find(u'市立') >= 0:
        split_name = name.replace(u'市立', u'市立 ').split(' ')
    elif name.find(u'町立') >= 0:
        split_name = name.replace(u'町立', u'町立 ').split(' ')
    elif name.find(u'村立') >= 0:
        split_name = name.replace(u'村立', u'村立 ').split(' ')
    else:
        split_name = (name, '')

    try:
        url = u'http://itp.ne.jp/result/?kw=' + urllib2.quote(u'福島県'.encode('shift_jis')) + '+' + urllib2.quote(town.encode('shift_jis')) + '+' + urllib2.quote(split_name[0].encode('shift_jis')) + '+' + urllib2.quote(split_name[1].encode('shift_jis'))
    except:
        print name
        continue
    
    fp = urllib2.urlopen(url)
    html = fp.read()

    tmp = open('tmp.txt', 'w')
    tmp.write(html)

    tmp = codecs.open('tmp.txt', 'r', 'shift-jis')

    flag = False

    for tmpline in tmp.readlines():
        if tmpline.find('class="icn"') >= 0:
            address = tmpline.split(u'　')[1].split('<a')[0]
            flag = True
            break

    if not flag:
        print name
        continue

    url = u"http://maps.google.com/maps/geo?q=" + address + u"+" + u"+福島県+日本&output=json&sensor=false"

    fp = urllib2.urlopen(url.encode("utf-8"))
    html = fp.read()

    for v in html.split("\n"):
        if v.find("coordinates") >= 0:
            break

    try:
        lng, lat = v.split("[ ")[1].split(", 0 ]")[0].split(", ")
        lng, lat = float(lng), float(lat)
    except:
        lng, lat = float('nan'), float('nan')

    print "%s\t%s\t%f\t%f" % (name, address, lat, lng)

そうすると、次のような出力が得られます。このやり方で 4400 箇所の 9 割以上は自動で経緯度の取得ができました。残った地点は、有志の皆さんのご協力で手入力で頑張りました。

福島市立御山小学校	福島県福島市御山字長滝１−１	37.778695	140.458548
福島市立三河台小学校	福島県福島市三河南町１７−７	37.757742	140.453800
福島市立松川小学校	福島県福島市松川町字南諏訪原３１−１	37.656167	140.464370
福島市立下川崎小学校	福島県福島市松川町沼袋字戸ノ内８３２−３	37.644727	140.500220
福島市立水原小学校	福島県福島市松川町水原字戸ノ内３１	37.665687	140.420651
福島市立金谷川小学校	福島県福島市松川町浅川字陳場２１	37.683573	140.465394

*1:当時、急いで勢い任せで書いたのを記事用にちょこっと修正しただけなので、もっと綺麗に書けます。

C++ 上では NULL な TObject の扱い

ROOT Python

PyROOT だと C の NULL が扱えないので、NULL な TObject を返すような ROOT の関数 (TFile::Get など) では、混乱する時があります。自分は毎回ググってるような…。

例えばここなんかに丁寧に説明がありますが、自分用の忘備録として blog にも書いておきます。

>>> import ROOT
>>> f = ROOT.TFile("test.root", "recreate")
>>> obj = f.Get("test")
>>> type(obj)
<class 'ROOT.TObject'>
>>> obj
<ROOT.TObject object at 0x0>
>>> obj == None
True
>>> obj is None
False
>>> not obj
True

こんな感じで、NULL の替わりに None が来ると思って is None で調べちゃ駄目。でも == None で調べるのは OK。ややこしい。