bosai.go.jp からユーザ登録のメールが届いたので早速、震源データをダウンロードし解析してみます。震源データは気象庁一元化震源リスト（http://www.hinet.bosai.go.jp/eq_inf/regUsers/jma/）を使っています。ダウンロードするにはユーザ登録しなくてはならないようだが、登録は簡単でした。

　このページでは最大一週間のデータしかダウンロードできません。まあ一週間だけでも 2000 件以上あるので、まずはこれを使ってみようかと思う。データ部分をコピーしテキストエディタにペースト。それからテキストに変換し保存します。

			
f$ = "aaaa_bb_cc_dd_ee_ff_gg_hhhhh_iiiiiiii_jjjjj_kkkkkkkk_lllll_mmmmmmm_nnnnn_ooooo";

				a = loadFixedTextData("E200209.txt",f$);

				msize(a);
			

　年月日時分秒をユリウス日に変換し、度数分布を計算。そして折線図を描くと地震頻度のヒストグラムができる（折線図のヒストグラムというのも変か？）。

t = jst2mjd(a[#,1|6]);
t[#,1] = t[#,1] - t[1,1]+1;
f = frequency(t[#,1],50);
changeCopyBitsMode(1);
parametricPlot(f[#,1],f[#,2]/(f[2,1]-f[1,1]))+0;

　このグラフは９月1日から７日まで観測された地震のヒストグラムです。日本では一日あたり 300〜400 回の地震が観測されているようです。この一週間、何となく地震頻度が減少しているような印象を受ける。

　長い期間でグラフを描くと月２回の周期で増減を繰り返しているので月の引力の影響だと思う←素人の想像だから真に受けないで！

　それではさっそく震源の位置（緯度／経度）、震源の深さ、マグニチュードからここ一週間の震源の散布図を三次元表示してみよう。三次元表示には Craig Kloeden 氏が開発した Rotater 5.0b1 を利用している。 CalendarMemo の mc では rotater 関数から利用できる。もし Rotater を持っていなければダウンロードして下さい。

b = change(a,|9,11,13,15|);
c = selectIf(b,b[#,1]>25 and b[#,1]<50 and b[#,2]>125 and b[#,2]<150);
c[#,1|3] = nzScore(c[#,1|3]);
m = int(c[#,4])+1;
rotater(c[#,1],c[#,2],-c[#,3]/4,m,0);

　驚いた事に、たった一週間の震源情報の散布図なのに日本の形が見えるようだ。すばらしい、日本って地震が多い国ですね。さて散布図を回転させ海溝付近のマントルが日本に潜り込んでいる状態をみてみましょう。

　これは日本を大東島の方から見た画像です。深い所の震源が日本列島の下に潜り込んでいるようすが判ります。教科書の通りだ！

　たった１週間の震源データなのに素人でもこれだけの事が判るという日本の地震観測体制ってすばらしい。それに一般に公開してくれるというのも素晴らしい。情報公開っていいですね〜

　もっと詳しく調べるため８月分のデータをダウンロードし、解析してみます。これから行うデータ処理は CalendarMemo の mc (matrix calculator) を使ったらこんなことができるという例です。

　処理の結果は間違いありません。でも私自身、地震に関しては高校生程度の知識なので解釈はたぶん間違えているでしょう。

f$ = "aaaa_bb_cc_dd_ee_ff_gg_hhhhh_iiiiiiii_jjjjj_kkkkkkkk_lllll_mmmmmmm_nnnnn_ooooo";
a = loadFixedTextData("EQ200208.txt",f$);
msize(a);
10736.0000
.OK.

一か月では地震観測回数は１万件をこえるんですね！すごい

t = jst2mjd(a[#,1|6]);
t[#,1] = t[#,1] - t[1,1]+1;
f = frequency(t[#,1],|31,1,61|);
changeCopyBitsMode(1);
parametricPlot(f[#,1],f[#,2]*2)+0;

一か月間の地震の発生頻度のグラフです。横軸が８月の日付け、縦軸が発生頻度（件／日）。一日を午前午後の２回計算しています。どうも午前と午後で上下しているようです。

tt = changeMN(f[1|60,1],2);
ff = changeMN(f[1|60,2],2);
m = nSum(ff);
m;
5527.0000 4919.0000
.OK.
m[1,1]/m[1,2]*100;
112.3602
.OK.

parametricPlot(tt[#,1],ff)+0;

　上記グラフの黒線が午前に発生した地震の頻度で赤線が午後に発生した地震の頻度です。８月は午前中が 5527 回、午後が 4949 回で、午前中の方が午後より１２％程発生頻度が多いようです。12％って偶然にしては多すぎるような気がしませんか？

　午前と午後の違いは何かと言えば太陽の位置くらいしかないでしょう。だから太陽の重力の影響？、それとも気温の影響？

　上記グラフを見ると地震発生頻度はおおよそ 250回／日〜500回／日 程度の変動をしています。これが震源の深さに関係あるのでしょうか。もし地震が他の影響を受けず独立して発生しているとするなら月内変動は震源の深さによって影響を受けないはずです。

　震源の位置を１Kmづつ区切って震源の深さと月内変動の関係を調べてみましょう。ここで月内変動の指標には毎日の発生頻度を一か月間まとめた標準偏差とします。

result = size(0,3);
for(s=0;s<50;s=s+1){
　aa = selectIf(a,a[#,13]>=s and a[#,13]<s+1);
　t = jst2mjd(aa[#,1|6]);
　t[#,1] = t[#,1] - t[1,1]+1;
　f = frequency(t[#,1],|31,1,31|);
　bb = |s,msize(t),sd(f[#,2])|;
　result = cat(result,bb);
　bb;
}
parametricPlot(result[#,1], result[#,3])+0;

　横軸が震源の深さ(Km) で縦軸が月内変動を示す標準偏差です。地下 10Km で変動がピークになり、地下 15Km より深くなると月内変動が小さくなっています。これは月内変動を示す地震は地下 15Km より浅いところで発生したものであるといえるのではないでしょうか。地下 15Km 付近を境に何か違うんでしょうね。

　それでは確認のため震源を地下 15Km で区切って上と下での地震発生頻度のグラフを描いてみます。

s = 15;
fm = |0,32,10,0,400,50|;
aa = selectIf(a,a[#,13]>=s);
t = jst2mjd(aa[#,1|6]);
t[#,1] = t[#,1] - t[1,1]+1;
f1 = frequency(t[#,1],|31,1,31|);
aa = selectIf(a,a[#,13]<s);
t = jst2mjd(aa[#,1|6]);
t[#,1] = t[#,1] - t[1,1]+1;
f2 = frequency(t[#,1],|31,1,31|);
x = addP(f1[#,2],f2[#,2]);
parametricPlot(f1[#,1],x)+0;

　横軸が８月の日付け、縦軸が発生頻度（件／日）です。赤線が震源の深さが 15Km より浅いところの地震発生頻度。黒が震源の深さが 15Km より深いところの地震発生頻度です。震源の深さが 15Km より深いところの地震の発生頻度のばらつきが少ないようです。

　これって何なんでしょうね？もう少し長期間のデータがあればもう少しましな処理ができるのですが。ちなみに８月の満月は２５日で１５日は半月でした。だからと言って月の影響？

　８月分の震源の分布を描いてみます。

b = change(a,|9,11,13,15|);
c = selectIf(b,b[#,1]>25 and b[#,1]<50 and b[#,2]>125 and b[#,2]<150);
pointMode(1);
setPlotWindowSize(400,500);
parametricPointPlot(c[#,2],c[#,1],|125,150,5,25,50,5|)+0;

　横軸が経度、縦軸が緯度の単なる散布図です。この散布図から震源に構造がみえるかどうか見てみましょう。

c1 = series(25,50,0.1)*10;
c2 = series(125,150,0.1)*10;
d = cross_count(c1,c2,c[#,1|2]*10);
d;
e = log(d+1);
e = flipV(e);
densityPlot(flipV(d),|10,0|);
densityPlot(e,|1,0|);
densityPlot(smooth(e),|1,0|);

　0.1 度範囲に発生した地震の回数を濃淡図で表示したものです。細かすぎるようで二枚目以降は対数変換してから表示しています。日本列島を東西に分けるような縦線が見えます。この画像からはこの程度しか判らないようです。

g = smooth(smooth(smooth(e)));
densityPlot(g,|1,0|);

frameLineDraw(0);
x = hokan(|0,250||0,250|,250,251)/250-0.5;
y = hokan(|0,0||251,251|,250,251)/251-0.5;
setParametricPlot3D(|450,450,0,60,10,9|);
RGB = popRGBColor(1);
x# = parametricPlot3DColor(x,y,g/10);
RGB = popRGBColor(0);
y# = parametricPlot3DColor(g,x,y,g/10,|1,0|);
x#*1+y#*0.3;

　まあ頻度を３Dで見てもあまりぱっとしないようです。でも四国付近と九州、いや関西より南はけっこう地震がありますね！　じゃ最後は rotater を使って関東のプレートの潜り込みをみてみますか

c = selectIf(b,b[#,1]>32 and b[#,1]<38 and b[#,2]>137 and b[#,2]<143);
c[#,1|3] = nzScore(c[#,1|3]);
m = int(c[#,4])+1;
rotater(c[#,1],c[#,2],-c[#,3]/4,m,0);

　以上のように、 CalendarMemo のデータ処理機能を利用するとインターネットで公開されているデータを使って色々な処理を試してみることができます。今回は気象庁で公開している震源データを使ってみました。情報公開制度というのはいいもんです。でも情報公開せよと叫んでも受け取る側にその能力がないと無駄でしょう。 CalendarMemo を使って受け取る能力を養いましょう（〜ん、最後はしかっり宣伝だ！）。