和歌山サイエンスカフェ「インフィニティ」

米国の航空宇宙局（NASA）の火星探査車
「パーシビアランス(Perseverance)」が無事到着したニュースが話題になっています。
昨日は、同探査機で録音された「火星の風の音」がネットに公開されて
これまた話題になっています。

NASAのサイトでは、二種類の音声がアップロードされています。
１つは探査機のノイズも込みのもの（生音に近いのかな？）と
もう１つは探査機のノイズをフィルター処理して、ノイズを低減したもののようです。
ぜひ、下記のサイトにアクセスして、火星の音を楽しんでください。
https://www.nasa.gov/press-release/nasa-s-mars-perseverance-rover-provides-front-row-seat-to-landing-first-audio

火星はかなたにあるお隣の惑星です。近づいたときでも、約６０００万キロメートルだそうです。
地球よりは少し小ぶりな星だそうですが、大気があり、風が吹き四季もあるそうです。
朝焼けや夕焼けもあるみたいです。
でも大気はとても薄くて、地表の大気は約７ｈＰａ程度だそうです。
地球は約１０００ｈＰａですから相当薄いですね。

NASAのPodcastで解説を聞くと、空気が薄いので音がソフトに聞こえるそうです。
火星の大気中での音の聞こえ方がシミュレートされて、それも聞くことができます。
https://soundcloud.com/nasa?fbclid=IwAR3oHdtvfMMIc-Vx0lcqNu_JVLruQQDJ36RlhaNjE2HsYIQDbgMdoLy190w

はるか遠くの星で聞こえる音を聞くことができるなんて不思議ですね。

録音されたのはDPA4006という無向性のマイクで市販されていたマイクだそうです。
https://www.hibino-intersound.co.jp/dpa_microphones/208.html

その音を聞きながら、スペクトログラムを描いてみたらどうかと思って試してみました。
スペクトログラムとは、音声波形を短時間フーリエ変換したもので、
周波数帯の時間変化をある程度確認できるものです。
（不確定性原理のため時間分解能と周波数分解能はトレードオフになります）

まず、ノイズ込みのもの波形を描くと、

です。
約１８秒のものです。
これのスペクトログラムを描いたものがこちらです。

専門家でもないので分かりませんが、
波形の振幅が大きいくなってるところで確かにスペクトルが大きくなっているように見えます。

特に工夫もせず、ハミング窓という窓関数をかけたものなので、最適な窓関数や
オーバーラップなども変えれば違う見え方がするかもしれません。

でも、こうやって最新の宇宙探査の一端に触れられるのはとても好奇心が刺激されますね。

一方のノイズ除去の音声の波形とスペクトログラムはこちらです。



たしかにノイズ成分が除去されているように見えます。
１６kHzあたりの波打つ線がなんなのか気になりますね。

火星探査の新たなニュースも楽しみです。


  

先週、インテルのデプスカメラRealsense D455が届きました。
注文してからなかなか発送されず、気をもんでいましたが、
デンバーを発ち、インディアナポリス、アンカレッジを経由して
日本に届いたようです。


横幅１５ｃｍほどの小さな箱が１個だけでした。

中身はこんな感じ。




基本的には赤外線カメラ２つを用いて、
カメラからの距離を測るようになっています。

これを使って、職場にAR砂場のような地形教育できるコンテンツを
おけないかと思っていました。

AR砂場ってのは、米国のUC DavisのオープンソースARsandboxや、
日本ではオープンソースのMagic-sandを使ったART Teknikaさんの「砂場マッピング」が
有名です。

「砂場マッピング」は２０１８年の佐賀県で開催された「ボクの土木展」で話題だったそうです。
また、寒地土木研究所さんはARsandboxを使って展示をされていたりするそうです。

ともかく、私はまだ現物を見たことはないのですが、動画や記事から、とても興味深く感じていました。

調べたところARsandboxはマイクロソフトのKinectにのみ対応。砂場マッピングは昨年D435対応版を
公開されたみたいですが、新しいD455には対応してません。

買うなら、最新の機種を、と思い購入しました。
で、あとは、自分で作れそうだと思い始めたところです。

開発に使う言語はPythonにしました。ずっとJavaでプログラミングしてましたが、
画像処理をするにあたってC++で開発されてるOpenCVに関する情報があまりなかったためです。
Pythonだと、RealsenseのSDKやOpenCVのラッパーに関する情報も多く、
また、等高線等を引くライブラリも充実しています。
Pythonは触ったことがある程度だったので、大学院での勉強もかねてです。

４日ほど、帰宅後にちょこちょと触ってたら、それなりに
リアルタイムで深度情報を使って等高線を引けるようになりました。

実際にプロジェクターを使ってAR砂場的なことができるのか試して、
職場の仲間にフィードバックをもらったところです。
開発としてはまだまだですが、方向性は見えてきました。



様子をYoutubeに上げてみました。
雰囲気分かっていただければ幸いです。

https://youtu.be/4PZBwr44SII






  

地震波形データを取り扱うための処理システムであるWINシステムの
インストールから始めます。
システムはubuntu20.04を想定。

０．開発環境を整える。
sudo apt install
を使って、g++,gcc,gfortran,tkなどを
インストールしておく。
他にもx11関連もインストールしなくてはいけない。libx11-devとか？
コンパイル時にエラーが出るので、その都度、
ググったりしながら対応する必要がある。

１．WINシステムのダウンロード
東京大学地震研究所のサイトにアクセス。
http://wwweic.eri.u-tokyo.ac.jp/WIN/

ダウンロードサイトから
WIN_pkg-3.0.6.tar.gz
（２０２０年７月１２日現在）
を取得。

２．ダウンロードしたファイルを/usr/localにコピーし、展開
sudo cp WIN_pkg-3.0.6.tar.gz
sudo tar zxvf WIN_pkg-3.0.6.tar.gz

３．展開したディレクトリ内に移動して、コンパイル
cd WIN_pkg-3.0.6
sudo ./configure --prefix=/usr/local/win
sudo make
sudo make install

４．できたファイルにパスを通す
vi .bashrc
などで設定を行う
export PATH=$PAHT:"/usr/local/bin"
を追加

５．テストをする

  

防災科学技術研究所が配信している
高感度地震観測網Hi-netってのがある。
https://www.hinet.bosai.go.jp/?LANG=ja
地震学の見習い坊主なので、詳しくとはよく分からない。

とりあえず、防災科技研が独自に拡張したWIN32形式っていうので
データが配信されている。

これを扱うために、東京大学地震研究所が開発したWIN形式というものに
変換する必要があるあしい。
そして、WINシステムというツール群で扱うのが、
とりあえず入門者がやってみるべきことらしい。

ということで、WINシステムのインストールや
防災科技研が配信しているwin32toolsを使ってみようという試みを始める。  

2017年初頭のモバイルノート（ウルトラブック）な
パソコンを比較検討中。
それなりのCPUパワーと、１６GB、５１２GB（できれば１TB）のSSD
で、１．３キロくらいより軽くて、コンパクトな電源アダプターなのが
必要条件だとする。
お金を出せれば候補も増えるので、２０万円以下くらいが候補。

１．MacBookPro　１３インチ（タッチバーなし）
　　メモリを１６GB、SSDを５１２GBにカスタマイズして、
　　１７９，８００円（税別、学割適用）
　　http://www.apple.com/jp_edu_1460/shop/buy-mac/macbook-pro?product=MLL42J/A&step=config#
　　いまなら、8,500円分のApple Storeギフトカードのプレゼント付き。

　　良い点
　　・英語キーボードが選べる。
　　・Retinaディスプレイになる。
　　・タッチバーなしモデルは、コスパは良好
　　悪い点
　　・CPUがskylake。秋くらいに Kabylakeになるとも言われている。）
・type-C に対応した、SDカードリーダとHDMIソケットのアダプターを買う必要が
　　　自分の場合ある。
　　

２．東芝 dynabook V82/B (東芝ダイレクトだと VZ72/B)
Core i7 7500U で８GBの512GBのSSD(Serial ATA)
２１４、０００円（税別） （東芝ダイレクトの会員になると割引あり）
https://dynabook.com/2in1-mobile-notebook-tablet/v-series/v82b-2016-winter-model-onyx-metallic-12-5-inch-high-end-convertible-pv82bmp-nja.html

　良い点
　・国内メーカーという点
　・いま所有してるR632は、非常に良いPCだった
　・ノングレアな画面で、タッチ対応
　・CPUはKabylake
　・HDMIもVGA接続も可能なアダプター付き（電源供給も可能）

　悪い点
　・英語キーボードなし
　・メモリ８GB設定なし
　・SSDに１TB設定なし。
　・ちょっとベゼル幅が大きめ

３．HP spectre x360（パフォーマンスモデル）
　　CPUはCore i7 7500U、メモリは16GBで、1TBのSSD（NVMe接続）
　　１６１，８００円（税別、学割適用、電話申し込みの場合）
HPの学割サイト　http://jp.ext.hp.com/campaign/personal/notebooks/2017student_pc_collection/
http://jp.ext.hp.com/directplus/personal/?jumpid=ps_74m86jxecn

　良い点
　・CPUはKabylake
　・１６GBも１TBも満たす。しかも、NVMe接続
　・コストパフォーマンつ高め
　・２月２８日までに申し込むと、Office Businessが無料

　悪い点
　・やや重め（１．３ｋｇ）
　・英語キーボードなし
　・電源アダプターのコードが太め
　・HPにやや悪い印象


他にもDELL(XPS13)とかLenovo(Carbon)なども候補にあがるかも。
隠し玉はXiaomiですが、一応、法的な問題もあるので。

spectrex360とXPS13との比較動画


spectre360とMacBookPro13との比較動画

  

windows7ノートに
Latexをインストールした。

Texインストーラ３を利用
http://www.math.sci.hokudai.ac.jp/~abenori/soft/abtexinst.html

参考にしたサイト
http://did2memo.net/2016/04/24/easy-latex-install-windows-10-2016-04/  

ここ数ヶ月、流体の数値計算の勉強を続けている。
２次元の河床変動プログラムを作るのが、当面の最終目標。

河床変動の計算の基礎式は、
１：流体の質量保存則（連続の式）
２；流体の運動量保存則（ナビエストークス方程式）
３；流砂の連続式
４；流砂量式
である。これら４つの方程式を解くことで
流れの場を一般的には求めることができる。

これらの式は複雑に絡み合い、特に２の運動量保存則は、
非線形の方程式（移流方程式）であり、一般に解析的に解くことが
できない。
そこで、数値計算が多用される。

数値計算においては、上記の方程式を離散化して、
解く方法がとられている。

その離散化の方法の違いにより、大きく
１：　格子法
２：　粒子法
に区別される
（他に個別要素法があるが、どう分類されるのかはまだ、勉強不足）

最近は２の粒子法も注目を集めているが、
計算範囲が大きい場合や、時間スケールを長くするには、
計算パワーが必要なことから、まだまだ一般的ではない。
粒子法については
『粒子法入門 流体シミュレーションの基礎から並列計算と可視化まで C/C++ソースコード付』
が詳しくて分かりやすい。

このため、一般には１の格子法が採用される。
格子法は、さらに、
A：差分法
B：有限体積法
C；有限要素法
などに分類される。

差分法は、微分方程式をテイラー展開を利用して離散化する。
また、有限要素法は、重み付きの関数を用いて離散化する。
（有限要素法も勉強不足）
一方で、有限体積法は、コントロールボリュームと呼ばれる
有限な体積の要素に積分することで、その収支を用いて離散化する。
有限体積法では、コントロールボリュームの質量収支等を考慮するので、
質量が保存されやすいという特があり、この点が他の計算方法との違いと
言われている。

連続体としての質量の保存を重視する流体力学の数値計算では、
有限体積法が人気を集めているようである。
他の工学分野（構造力学等）では、有限要素法が一般的なようである。

以上のような理由から、有限体積法を用いて、
数値計算（数値シミュレーション）を進めていくこととした。


なお、使用言語は、
Javaを使っていくことにする（途中で変更するかもしれない）
もともと数値計算でよく用いられているFortranも使っていたが、
Androidアプリ開発のため、Javaを勉強し始めたことと、
オブジェクト指向を理解するためにも、Javaをプログラミング言語として
学んでいくことは、今後に役立つと考えたためである。


有限体積法の勉強に用いた参考文献を以下に示す。








  

GMT5.3.1を使い、国土地理院の基盤地図情報サービスの
１０mメッシュのDEMから、地形の陰影図や鳥瞰図を作る。

まずはデータのダウンロード
http://www.gsi.go.jp/kiban/
（今回は紀伊半島の一部のデータを利用）

登録してダウンロードし解凍してできる
FG-GML-5035-37-dem10b-20161001.xml
を利用します。

GMTで扱えるgrd形式に変換する必要がある。
grd形式に変換するスクリプトを
産業技術研究所の野田さんが公開されているので
利用させていただく。ありがとうございます。
https://staff.aist.go.jp/a.noda/programs/jpgis2grd/jpgis2grd.html

このスクリプトを動かすとgrd形式に変換できるが
一部環境に合わせて変更する。
（commandを「gmt xyz2grd」とか出来上がるファイル名などの変数を）

$ perl jpgis2grd.pl
を実行
しかし、エラーで止まる。
途中のどこかに３列目の数値がないところがあるようだ。
（今回は 843751行目）

このスクリプトのglobal optionに
my $outfile2 = "out2.xyz"
を追記。

grdに変換するあたりを

print STDOUT "Making a grd file.\n";
# my $command="gmt xyz2grd $outfile -R$xmin/$xmax/$ymin/$ymax -I$xint/$yint -G$grdfile -N-9999 -Dlongitude/latitude/height/1/0/$opts{title}/$opts{remark} ";
my $command_sed="sed -e '843751d' $outfile > $outfile2";
my $command="gmt xyz2grd $outfile2 -R$xmin/$xmax/$ymin/$ymax -I$xint/$yint -G$grdfile -di-9999 -Dlongitude/latitude/height/1/0/$opts{title}/$opts{remark} ";
if ($opts{verbose}) { print STDERR "Running \"$command\"\n"; };
system("$command_sed");
system("$command");
system("gmt grdinfo $grdfile");

と一部修正。
（perlに詳しくないので、あんまりいいやり方ではなさそうですが）

再度、
$ perl jpgis2grd.pl -v -n
を実行
（実行状況や中間ファイルを残すのに　-v -nオプションをつけた）

すると、out.grdファイルができるので
とりあえず、これを
FG-GML-5035-37-dem10b-20161001.grd
と名前を変更。

次はgmtで陰影図を描く
以下のスクリプトで
#!/bin/bash
#
# DEM10mメッシュのgrdファイルから陰影図を作る。
# 2016/12/12
#
# perl jpgis2grd.pl で出力される範囲をコピー
# xmin=135.87500000, xmax=136.00000000, ymin=33.58333333, ymax=33.66666667

grd_file=FG-GML-5035-37-dem10b-20161001.grd
shade_file=5035-37.shade
eps_file=5035-37_inei.eps
cpt_file=color.cpt
ctb=-1000/1000/50 #color table min/max/interval
region=135.875/136.00/33.5833333/33.6666667

gmt makecpt -Cglobe -T$ctb -Z > $cpt_file

gmt gmtset FORMAT_GEO_MAP D

gmt grdgradient $grd_file -A90 -G$shade_file -Ne0.6 -V
gmt grdimage $grd_file -R$region -JM15 -C$cpt_file -I$shade_file -X4 -Y4 -K > $eps_file
gmt psscale -D7.5/-1/12/0.3h -C$cpt_file -Ba200f100:"altitude(m)": -I -V -K -O >> $eps_file
gmt psbasemap -R$region -JM15 -Ba0.05f0.025g0.01 -O -V >> $eps_file

こんな図ができる



次は鳥瞰図
#!/bin/bash
#
# DEM10mメッシュのgrdファイルから陰影図を作る。
# 2016/12/12
#
# perl jpgis2grd.pl で出力される範囲をコピー
# xmin=135.87500000, xmax=136.00000000, ymin=33.58333333, ymax=33.66666667

grd_file=FG-GML-5035-37-dem10b-20161001.grd
shade_file=5035-37.shade
d3_file=5035-37_3d.eps
cpt_file=color.cpt
ctb=-1000/1000/50 #color table min/max/interval
region=135.875/136.00/33.5833333/33.6666667
boundary_cap=wSnEZ
gmt makecpt -Cglobe -T$ctb -Z > $cpt_file
gmt gmtset FORMAT_GEO_MAP D
gmt grdgradient $grd_file -A90 -G$shade_file -Ne0.6 -V
gmt grdview $grd_file -Jm130 -Jz0.002 -C$cpt_file -I$shade_file -p150/25 -Ba0.05f0.025g0.01 -Bza400 -B$boundary_cap -Qi300 -N0+ggray -V > $d3_file

を実行すると、鳥瞰図ができる。



角度とかは適宜設定する。

  

久しぶりに綺麗なグラフとかを描きたくて
GMTを触ってみた

バージョンが５に上がっていて、いろいろ
変わったところが多い。
インストールはubuntuの場合cmakeでしたり
するようになっていたり。

プログラミングもコマンドからモジュールという
概念に変わっているようだ。
例えば
gmt pscoast ほにゃらら
と入力する。

より日本語の情報が少ないので
できたことはアップしてみる。

今回は紀伊半島周辺の海岸線図に
県境と、距離スケールと方位を示す矢印を追記する。
矢印、というかベクトルを書くのに手間取った。

gmtのバージョン5以降は
インストールできているとして、
県境のデータをダウンロード。
こちらを参考に
http://www.geocities.jp/ne_o_t/otenki.htm

後は、観測地点なんかも記載できたらいいので、
std_name.dataと名付けたファイル
　　NachiRiver 33.668 135.904
と、std＿point.dataと名付けたファイル
　　135.904 33.668 LM 0.0 12p,Helvetica,255/0/0 NachiRiver
を用意。

今回はmacで作業してるので、あんまりgmtのファイルが
どこにあるかわからないので、shareと名付けたファイルに
県境のデータ ken.txtを置いて、以下のようなシェルスクリプトを書いて実行。

#!/bin/bash
# kenzakai_wakayama.sh

rangeSW=134.5/136.5/33/35
scale=M10c
afg=a1f0.5g0
ruler=136/33.3/33.2/50k
pen1=0.6p,black #[width,color]
pen2=0.4p,50
g_color=200
s_color=255

point_data=stn_point.dat
name_data=stn_name.dat

fig=Wakayama.eps

gmt gmtset FORMAT_GEO_MAP +DF
gmt set MAP_VECTOR_SHAPE 0

gmt psbasemap -J$scale -R$rangeSW -B$afg -X3c -Y3c -P -V -K > $fig
gmt pscoast -J -R -B$afg -Lf145/30/35/500 -Dh -P -W$pen1 -G$g_color -S$s_color -L$ruler -V -O -K >> $fig
gmt psxy ../share/ken.txt -J -R -W$pen2 -O -K >> $fig
gmt psxy $point_data -i2,1 -J -R -Sc0.2c -G255/0/0 -V -O -K >> $fig #地点の位置
gmt pstext $name_data -J -R$rangeSW -D0.5/0 -F+j+a+f -O -K >> $fig #地点の名前 -Dはシフト

#方向矢印の記述
#矢印
gmt psxy -R -J -SV0.4+s+e+a45 -W1 -G0 -O -K << EOF >> $fig #片矢印
135,33.3,135,33.65 #始点、終点
EOF
#横棒
gmt psxy -R -J -SV0.5+s -W0.8 -G0 -O -K << EOF >> $fig #線分を引く
134.93,33.45,135.07,33.45 #始点、終点
EOF
#Nの文字
gmt pstext -R$rangeSW -J -D0/0.1 -F+j+a+f -O << EOF >> $fig
135 33.65 BC 0.0 12p,Helvetica,black N
EOF

#白抜きのpngの作成
#convert -density 100 test2.ps test2.png
# mogrify -trim -density 300 -bordercolor 'transparent' -border 10x10 -format png *.eps
#mogrify -trim -density 300x300 -transparent white -format png $fig


以上。
最後でコメントアウトしてるけど、
imagemagicを使って、白のところを透明にしたpngファイルも
作成できるようにしている。

こんな図ができると思う。（クリックして拡大）


ベクトルを書くのに、psxyの-SVオプションが結構戸惑った。
WANtaroHP (Programing on Mac)さん
http://kbtkt.web.fc2.com/sub_0_gmt5.html
を参考にさせてもらった。ありがとうございます。  

礫の移動を表示するため
QGISに、
WKTのテキストデータを
レイヤー追加して、ラインを引いた。

矢印表示にできるように
ラインのプロパティを以下のように設定。



こんな感じです。（赤の矢印）


画像をクリックしたら、大きく表示されます。  

３年ほど使ってるMac mini(Late2012)の
ハードディスクをSSDに交換してみました。
ずっと思ってたけど、ちょっと難しく思っていたので。

このために買ったのは、

と、
交換で出てくるHDDのための外付けケース

と、
ちょっと特殊なドライバーが必要だということで


SSDが値上がり気味ってことなので、
慌てて買いました。

件のMacmini

一年ほど前にメモリは８GBに交換済み。

裏返すと

黒い蓋をヒネって開ける。


購入したドライバーでビスを外して、
ごそごそとファンを外す。




金属の金網みたいものを外す。この時
wifiのアンテナ線も外します。


すると、黒いカバーに覆われたHDDが見える。
ごそごそと引き抜く。
自分のMacminiは、下側（ちゃんと置くと）にあったので、取りやすい方だった。
上側にあるものもあるらしい。そうなると、作業がさらに増えるみたい。



HDDについてる部品を外して、SSDに装着して、
手順を遡って、元に戻すだけ。


でも、wifiのアンテナの接続部分とか、
fanの電源端子が小さくて、つけるのは結構手間取った。
ピンセットがあると便利。

SSDには予めSierraを外付けのHDDケースを使って、
インストールおいたので、すぐに起動できた。
MacはOSを簡単にダウンロードして、外付け
のドライブにインストールできるのがすごいな。

この交換作業は色々ネットに情報があるので、
今更ですが。
Mac miniはLate2012に根強い人気があるらしい。  

３年ぶりにスマホを買い替えました。
Nexus5からの買い替えで、
Asusのzenfone3(ZE520KL)
にしました。


電話はガラケーなので、
スマホはocnのモバイルワン（月約１０００円コース）
の通信線専用です。
simフリー版しか興味ないところです。
microSDも入れたままなので、
DSDSにも興味はないです。

なるべくコスパがいいものをってことで
zenfone3がmotoGplusと悩みましたが、
後者は電子コンパスなしってことで諦めました。

並行輸入したら安いみたいですが、
あんまり文句を言うところがないのは嫌なので
適技の問題もあるかもしれないので、
NTTXstoreで購入。
同時に、モバイルoneも買い替えました。

いろいろレビューはあるので、そちらを見てもらったら
いいですが、
かなり満足です。
指紋認証もとてもスムーズです。

ただ、lightroomモバイルを使っての
raw撮影（dng撮影）は一応できますが、
カラーマップがだめらしくて、使い物になりません。

あと、android6でありながら、
マルチユーザー（マルチアカウント）が使えません。
これは結構いたかった。

でも、その他は、十分だと思います。

USBのtypeCなので、
変換コネクタ

保護フィルム

ケース

も買いました。



  

河川シュミレーションソフトのiRICを
ひさしぶりに動かしてみました。
http://i-ric.org/ja/

まだあまり慣れてないので、
河床変動も計算できるＮａｙｓ２ＤＨソルバーを
チュートリアルで練習してみた。

ものの数時間で基本的なところは分かるみたい。
面白いですよ。

河の形を与えて、上流からの流量とかを
設定するだけ。





１１月２６日の和歌山サイエンスカフェはまだまだ
参加者募集中です！！

よろしくお願いします。
詳しくは以下のリンクを！
http://wscinfinity.ikora.tv/e1234831.html  

QGISでラインを引くにはWKTという形で記述した
CSVファイルを用意するとよい。

こんな感じのファイル

ID WKT hoge
1 LINESTRING(-9634.725 -260090.329,-9634.715 -260090.363) hogehoge
2 LINESTRING(-9637.913 -260090.775,-9637.916 -260090.777) hogehogehoge

これを
「レイヤー」
「レイヤーの追加」
「デリミテッドテキストファイルの追加」
を選び、
WKT(Well Known Text)形式を選択する。  

気象庁が出している土壌雨量指数を
読み込むプログラムをC#の勉強を
兼ねて作ってみた。

全然、わからなかったので、
「ロケッこが行く」のブログ
http://blog.syo-ko.com/?eid=1228
を大変参考にさせてもらった。
ありがとうございます。

grib2形式のバイナリーファイルでできている。
データ部はランレングスという形で表現されていて
この読み方に結構時間を取られた。

ちなみに、気象庁が出している技術資料は
http://www.data.jma.go.jp/add/suishin/jyouhou/pdf/302.pdf
から取れる。
でも、タンクの数が違ったり、まだまだ謎が多い。
自分の理解不足か、よくわからない。

ちょっとづつC#の勉強しながらなので、
全然クラスの使い方とか、ファイルの読み込み方とか
統一とれてない上に、動作確認のためのコンソールへの
出力とかをコメントアウトのまま残してる。

ネットで検索すると、土壌雨量指数の読み込みプログラムが、
異常な値段で売られてたりするので、
少しでも役に立てればうれしい。

ただ、実行は自己責任でお願いします。

できたテキストファイルをRで描画すると

それっぽいので、たぶんうまくいけてる。

で、ソースは以下のとおり。

using System;
using System.IO;

class DojoUryouANAL
{


static void Main()
{

int FileYear = 2011;　//入力年の指定
int[] nMonth = new int[] { 31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31 };
nMonth[2 - 1] = 28 + (1 / (FileYear % 4 + 1)) * (1 - 1 / (FileYear % 100 + 1)) + (1 / (FileYear % 400 + 1)); //２月のうるう年判別

int nCount = 0;
int FileMonth = 1;
int FileDay = 1;
int FileHour = 0;
int FileMinute = 0;
int FileSecond = 0;

string InputFile1 = "InputFiles¥¥Z__C_RJTD_";
string InputFIle2 = "_SRF_GPV_Gll5km_P-swi_ANAL_grib2.bin";

string[] InputFileNameList = new string[18000];

while (FileMonth < 13)
{
while (FileDay < nMonth[FileMonth - 1] + 1)
{
while (FileHour < 24)
{
for (int minIndex = 0; minIndex < 2; minIndex++)
{
if (minIndex == 0)
{
FileMinute = 0;
InputFileNameList[nCount] = InputFile1 + FileYear.ToString("D4") + FileMonth.ToString("D2") + FileDay.ToString("D2") + FileHour.ToString("D2") + FileMinute.ToString("D2") + FileSecond.ToString("D2") + InputFIle2;
//Console.WriteLine("InputFile = " + InputFileName);
// Console.WriteLine("Year={0},Month={1},Day={2},Hour={3},Minute={4},Second={5}", Year, Month, Day, Hour, Minute, Second);
nCount++;
}
else
{
FileMinute = 30;
InputFileNameList[nCount] = InputFile1 + FileYear.ToString("D4") + FileMonth.ToString("D2") + FileDay.ToString("D2") + FileHour.ToString("D2") + FileMinute.ToString("D2") + FileSecond.ToString("D2") + InputFIle2;
//Console.WriteLine("InputFile = " + InputFileName);
//Console.WriteLine("Year={0},Month={1},Day={2},Hour={3},Minute={4},Second={5}", Year, Month, Day, Hour, Minute, Second);
nCount++;
}
}
FileHour++;
}
FileHour = 0;
FileDay++;
}

FileDay = 1;

FileMonth++;

}


for (int sCount = 0; sCount < nCount; sCount++)
{

string InputFileName = InputFileNameList[sCount];
//string InputFileName = "Z__C_RJTD_20110903170000_SRF_GPV_Gll5km_P-swi_ANAL_grib2.bin";
//string OutputFileName = "output.txt";
string FileHeadS = InputFileName.Substring(11, 10);
string FileYearS = InputFileName.Substring(21, 4);
string FileMonthS = InputFileName.Substring(25, 2);
string FileDayS = InputFileName.Substring(27, 2);
string FileHourS = InputFileName.Substring(29, 2);
string FileMinuteS = InputFileName.Substring(31, 2);
string FileSecondS = InputFileName.Substring(33, 2);


//string OutputFileName = "Dojosisu_ANAL_" + FileYearS + FileMonthS + FileDayS + FileHourS + FileMinuteS + ".txt";
//Console.WriteLine("OutputFileName = " + OutputFileName);

FileStream dataIN; //入力データファイル
StreamWriter fstr_out; //出力データファイル

//ファイルの読み込み（はじめ）
try
{
dataIN = new FileStream(InputFileName, FileMode.Open, FileAccess.Read);
}
catch (IOException exc)
{
Console.WriteLine(exc.Message);
return;
}

int fileSize = (int)dataIN.Length; //オクテット数
byte[] rByte = new byte[fileSize]; //rByte 全オクテットの保持配列

Console.WriteLine("File Size: " + fileSize);

int readSize;
int remain = fileSize;
int bufPos = 0; //読み込み開始ポイント

while (remain > 0)
{
readSize = dataIN.Read(rByte, bufPos, Math.Min(1, remain));

bufPos += readSize;
remain -= readSize;
}
//ファイルの読み込み（終わり）


//ファイルを閉じる。
dataIN.Dispose();




//////解釈の始まり//////
//１〜４オクテット
// for (int i = 0;i < 4; i++){
// Console.Write((char)rByte[i]+",");
// }

//９〜１６オクテット
// long rValue = ByteToInt(rByte, 8, 8);
long rValue = ByteToSingedInt(rByte, 8, 8);
Console.Write(rValue + ",");



Console.WriteLine();
Console.WriteLine("第１節");
//第１節１〜４オクテット ※節の長さ
rValue = ByteToSingedInt(rByte, 16 + 1 - 1, 4);
Console.WriteLine(rValue);
int lenFirstchp = (int)rValue;
//第１節５オクテット
rValue = ByteToSingedInt(rByte, 16 + 5 - 1, 1);
Console.WriteLine(rValue);
//第１節６〜７オクテット
rValue = ByteToSingedInt(rByte, 16 + 6 - 1, 2);
Console.WriteLine(rValue);
//第１節１３〜１４オクテット
rValue = ByteToSingedInt(rByte, 16 + 13 - 1, 2);
Console.WriteLine(rValue);
//第１節１５オクテット
rValue = ByteToSingedInt(rByte, 16 + 15 - 1, 1);
Console.WriteLine(rValue);
//第１節16オクテット
rValue = ByteToSingedInt(rByte, 16 + 16 - 1, 1);
Console.WriteLine(rValue);
//第１節17オクテット
rValue = ByteToSingedInt(rByte, 16 + 17 - 1, 1);
Console.WriteLine(rValue);
//第１節18オクテット
rValue = ByteToSingedInt(rByte, 16 + 18 - 1, 1);
Console.WriteLine(rValue);
//第１節19オクテット
rValue = ByteToSingedInt(rByte, 16 + 19 - 1, 1);
Console.WriteLine(rValue);
//第１節２０オクテット
rValue = ByteToSingedInt(rByte, 16 + 20 - 1, 1);
Console.WriteLine(rValue);


Console.WriteLine();
Console.WriteLine("第３節");

//第３節１〜４オクテット　※節の長さ
rValue = ByteToSingedInt(rByte, 16 + lenFirstchp + 1 - 1, 4);
Console.WriteLine(rValue);
int lenThirdchp = (int)rValue;
//第３節５オクテット　※節番号
rValue = ByteToSingedInt(rByte, 16 + lenFirstchp + 5 - 1, 4);
Console.WriteLine(rValue);
//第３節３１〜３４オクテット
rValue = ByteToSingedInt(rByte, 16 + lenFirstchp + 31 - 1, 4);
Console.WriteLine(rValue);
//第３節３５〜３８オクテット
rValue = ByteToSingedInt(rByte, 16 + lenFirstchp + 35 - 1, 4);
Console.WriteLine(rValue);
//第３節４７〜５０オクテット
rValue = ByteToSingedInt(rByte, 16 + lenFirstchp + 47 - 1, 4);
Console.WriteLine(rValue / 1000000.0);
//第３節５１〜５４オクテット
rValue = ByteToSingedInt(rByte, 16 + lenFirstchp + 51 - 1, 4);
Console.WriteLine(rValue / 1000000.0);
//第３節５６〜５９オクテット
rValue = ByteToSingedInt(rByte, 16 + lenFirstchp + 56 - 1, 4);
Console.WriteLine(rValue / 1000000.0);
//第３節６０〜６３オクテット
rValue = ByteToSingedInt(rByte, 16 + lenFirstchp + 60 - 1, 4);
Console.WriteLine(rValue / 1000000.0);
//第３節６４〜６７オクテット
rValue = ByteToSingedInt(rByte, 16 + lenFirstchp + 64 - 1, 4);
Console.WriteLine(rValue / 1000000.0);
//第３節６８〜７１オクテット
rValue = ByteToSingedInt(rByte, 16 + lenFirstchp + 68 - 1, 4);
Console.WriteLine(rValue / 1000000.0);



int lenNow = 0;

for (int NumTank = 0; NumTank < 4; NumTank++)
{



Console.WriteLine();
Console.WriteLine("第４節");



//第４節１〜４オクテット　※節の長さ
rValue = ByteToSingedInt(rByte, 16 + lenFirstchp + lenThirdchp + lenNow + 1 - 1, 4);
Console.WriteLine(rValue);
int lenForthchp = (int)rValue;

//第４節5オクテット
rValue = ByteToSingedInt(rByte, 16 + lenFirstchp + lenThirdchp + lenNow + 5 - 1, 1);
Console.WriteLine(rValue);

Console.WriteLine();
Console.WriteLine("第５節");

//第５節１〜４オクテット　※節の長さ
rValue = ByteToSingedInt(rByte, 16 + lenFirstchp + lenThirdchp + lenForthchp + lenNow + 1 - 1, 4);
Console.WriteLine(rValue);
int lenFivechp = (int)rValue;

//第５節５オクテット
rValue = ByteToSingedInt(rByte, 16 + lenFirstchp + lenThirdchp + lenForthchp + lenNow + 5 - 1, 1);
Console.WriteLine(rValue);
//第５節６−９オクテット
rValue = ByteToSingedInt(rByte, 16 + lenFirstchp + lenThirdchp + lenForthchp + lenNow + 6 - 1, 4);
Console.WriteLine(rValue);

//第５節１０−１１オクテット
rValue = ByteToSingedInt(rByte, 16 + lenFirstchp + lenThirdchp + lenForthchp + lenNow + 10 - 1, 2);
Console.WriteLine(rValue);

//第５節１２オクテット
rValue = ByteToSingedInt(rByte, 16 + lenFirstchp + lenThirdchp + lenForthchp + lenNow + 12 - 1, 1);
Console.WriteLine("nbit: " + rValue);
int nbit = (int)rValue;

//第５節１３−１４オクテット
rValue = ByteToSingedInt(rByte, 16 + lenFirstchp + lenThirdchp + lenForthchp + lenNow + 13 - 1, 2);
Console.WriteLine("maxv: " + rValue);
int maxv = (int)rValue;

//第５節１５−１６オクテット
rValue = ByteToSingedInt(rByte, 16 + lenFirstchp + lenThirdchp + lenForthchp + lenNow + 15 - 1, 2);
Console.WriteLine("max_M: " + rValue);
int max_M = (int)rValue;

//第５節１７オクテット
rValue = ByteToSingedInt(rByte, 16 + lenFirstchp + lenThirdchp + lenForthchp + lenNow + 17 - 1, 1);
Console.WriteLine(rValue);

int[] LevelValue = new int[max_M + 1];
for (int level = 1; level < max_M + 1; level++)
{
//第５節１８−オクテット レベルの対応値
rValue = ByteToSingedInt(rByte, 16 + lenFirstchp + lenThirdchp + lenForthchp + lenNow + 18 + (level - 1) * 2 - 1, 2);
LevelValue[level] = (int)rValue / 10;
Console.WriteLine("level: " + level + " Value: " + rValue);
Console.WriteLine("level: " + level + " LevelValue " + LevelValue[level]);
}
LevelValue[0] = 9999; ///欠測地

Console.WriteLine();
Console.WriteLine("第６節");


//第６節１−４オクテット ※節の長さ
rValue = ByteToSingedInt(rByte, 16 + lenFirstchp + lenThirdchp + lenForthchp + lenFivechp + lenNow + 1 - 1, 4);
Console.WriteLine(rValue);
int lenSixchp = (int)rValue;

//第６節５オクテット
rValue = ByteToSingedInt(rByte, 16 + lenFirstchp + lenThirdchp + lenForthchp + lenFivechp + lenNow + 5 - 1, 1);
Console.WriteLine(rValue);
//第６節６オクテット
rValue = ByteToInt(rByte, 16 + lenFirstchp + lenThirdchp + lenForthchp + lenFivechp + lenNow + 6 - 1, 1);
Console.WriteLine(rValue);


Console.WriteLine();
Console.WriteLine("第７節");


//第７節１−４オクテット
rValue = ByteToSingedInt(rByte, 16 + lenFirstchp + lenThirdchp + lenForthchp + lenFivechp + lenSixchp + lenNow + 1 - 1, 4);
Console.WriteLine("length of this chapter: " + rValue);
int n_data = (int)rValue - 6 + 1;
int lenSevenchp = (int)rValue;

//第７節５オクテット
rValue = ByteToSingedInt(rByte, 16 + lenFirstchp + lenThirdchp + lenForthchp + lenFivechp + lenSixchp + lenNow + 5 - 1, 1);
Console.WriteLine("number of chapter: " + rValue);



//出力ファイル名の操作
string OutputFileName = "OutputFIles¥¥Dojosisu_ANAL_" + "TANK" + NumTank.ToString() + "_" + FileYearS + FileMonthS + FileDayS + FileHourS + FileMinuteS + "UTC.txt";
Console.WriteLine("OutputFileName = " + OutputFileName);


int j;
int p = -1, m = 1, n = 0, k = 0;
int v = 0;
int l = ipow(2, nbit) - 1 - maxv;

// n_data = 30;

try
{
fstr_out = new StreamWriter(OutputFileName); //土壌雨量指数の出力ファイル
}
catch (IOException exc)
{
Console.WriteLine(exc.Message);
return;
}

for (int iRun = 1; iRun < n_data + 1; iRun++)
{
// Console.WriteLine("iRun: " + iRun);

//第７節６−オクテット
// rValue = ByteToSingedInt(rByte, 16 + 21 + 72 + 58 + lenFivechp + 6 + iRun + 5 - 1, 1);
// Console.WriteLine("iRun: " + iRun + " Value: " + rValue);
// Console.WriteLine("iRun: " + iRun + " Value: " + rByte[16 + 21 + 72 + 58 + lenFivechp + 6 + iRun + 5 - 1]);
// Console.WriteLine("iRun: " + iRun + " Value: " + rByte[16 + 21 + 72 + 82 + lenFivechp + 6 + iRun + 5 - 1].ToString("X"));

v = (int)ByteToInt(rByte, 16 + lenFirstchp + lenThirdchp + lenForthchp + lenFivechp + lenSixchp + lenNow + iRun + 5 - 1, 1);
// Console.Write(v + ",");
if (v <= maxv)
{
if (p >= 0)
{
// Console.WriteLine(" ,m: " + m);

for (j = 0; j < m; j++)
{
k++;
// Console.Write(p+" , ");
try
{
fstr_out.Write(LevelValue[p] + " , ");
}
catch (IOException exc)
{
Console.WriteLine(exc.Message);
break;
}
}
// Console.WriteLine("p: " + p + " ,m: " + m );
p = v;
m = 1;
n = 0;
}
p = v;
m = 1;
n = 0;

}
else
{
m += ipow(l, n) * (v - maxv - 1);
// Console.WriteLine(ipow(l, n) * (v - maxv - 1)+","+m);
n++;
}

}

for (j = 0; j < m; j++)
{
k++;
if (k != 286720)
{
// Console.Write(p);
try
{
fstr_out.Write(LevelValue[p] + " , ");
}
catch (IOException exc)
{
Console.WriteLine(exc.Message);
break;
}
}
else
{
// Console.Write(p+" , ");
try
{
fstr_out.Write(LevelValue[p] + " ");
}
catch (IOException exc)
{
Console.WriteLine(exc.Message);
break;
}
}
}
Console.WriteLine("lenNow:" + lenNow);
lenNow += lenForthchp + lenFivechp + lenSixchp + lenSevenchp;
Console.WriteLine("lenNow:" + lenNow);
fstr_out.Close();
Console.WriteLine("k: " + k); //286,720=512x560

}

//第８節１−４オクテット ※-1
rValue = ByteToInt(rByte, 16 + lenFirstchp + lenThirdchp + lenNow + 0 - 1, 1);
Console.WriteLine("¥n -1");
Console.WriteLine("Message: " + rValue);


//第８節１−４オクテット ※1
rValue = ByteToInt(rByte, 16 + lenFirstchp + lenThirdchp + lenNow + 1 - 1, 1);
Console.WriteLine("¥nEighth Chapter 1");
Console.WriteLine("Message: " + rValue);

//第８節１−４オクテット ※2
rValue = ByteToInt(rByte, 16 + lenFirstchp + lenThirdchp + lenNow + 2 - 1, 1);
Console.WriteLine("¥nEighth Chapter 2");
Console.WriteLine("Message: " + rValue);

//第８節１−４オクテット ※3
rValue = ByteToInt(rByte, 16 + lenFirstchp + lenThirdchp + lenNow + 3 - 1, 1);
Console.WriteLine("¥nEighth Chapter 3");
Console.WriteLine("Message: " + rValue);

//第８節１−４オクテット ※4
rValue = ByteToInt(rByte, 16 + lenFirstchp + lenThirdchp + lenNow + 4 - 1, 1);
Console.WriteLine("¥nEighth Chapter 4");
Console.WriteLine("Message: " + rValue);

}
}
//Mainメソッド　終了


//整数の累乗を計算するメソッド
static int ipow(int i, int j)
{
int k, l;
for (k = 0, l = 1; k < j; k++) l *= i;
return (l);
}


//バイトを、整数に変換するメソッド
static long ByteToInt(byte[] source, int Start, int Length)
{
long rVal = 0;
int intB;
int intE = Start + Length;

for (intB = Start; intB < intE; intB++)
{
rVal *= 256;
rVal += Convert.ToInt64(source[intB]);

//Console.WriteLine("intB: " + intB + " ,souce: " + Convert.ToInt64(source[intB]) + " ,rval: " + rVal);
}

return rVal;

}


//バイトを、符号付き整数に変換するメソッド
static long ByteToSingedInt(byte[] source, int Start, int Length)
{
long rVal = 0;
int intB;
int intE = Start + Length;
bool bSng = true;

//先頭オクテットの確認
if ((source[Start] & 0x80) == 0x80)
{
// Console.WriteLine("負数である");
bSng = true;
source[Start] = (byte)(source[Start] & 0x7F);
}
else
{
bSng = false;
// Console.WriteLine("正数である");
}


for (intB = Start; intB < intE; intB++)
{
rVal *= 256;
rVal += Convert.ToInt64(source[intB]);

//Console.WriteLine("intB: " + intB + " ,souce: " + Convert.ToInt64(source[intB]) + " ,rval: " + rVal);

}

if (bSng == true)
{
rVal *= -1;
}

return rVal;

}


}

  

今日から「まいにちスペイン語」を
開始しました。NHKのラジオ講座です。

放送は朝の７時１５分から３０分
もしくはお昼の２時４５分から３時です。

４年前に一度半年間聞いていましたが、
それっきりで、今度お世話になる研究室に
スペイン語圏の方がいらっしゃるってことで、
これはチャンスだと学び直すことにしました。

英語は、学校教育で読み書きから入りましたが、
スペイン語は、耳としゃべりで学習してみたいと思ってます。
スペイン語はなんだか陽気で、鬱な気分も吹っ飛び感じで
たのしいですよぉー！

学生時代はロシア語、ドイツ語もかじりましたが、中途半端だったので、
スペイン語を選んでみました。



kidle本だと、かなり安く手に入るみたいです。
  

和歌山大学や近畿地方整備局が主催の「防災ジオツアー」に
参加してきました。

平成２３年に甚大な土砂災害のあった那智川流域（和歌山県那智勝浦町）
のジオサイトや砂防堰堤を巡りながら、防災に関わる知識を深めようという
趣旨のイベントでした。
参加費１０００円で、４０名程度の参加者でした。


初めは、河川の氾濫により多くの被害が出た井関地区。
災害の怖さを忘れないための碑で、被害や状況の説明を受けました。
ここには当時の様子を伝える浸水深３．５５ｍを記したポールも立っています。
水流の凄まじさに思いを馳せつつ、参加者全員で黙祷を捧げました。


次は大きな土石流が発生した金谷山の崩壊地（源頭部）が望める場所に。

崩壊地付近には、上部に火成岩（花崗斑岩）、下部に堆積岩（熊野層群）の境界が
見えます。そこは、小さな滝になっていました。
この境界付近が大きな土石流の要因になっているのでは、と考えられています。
この付近では、学生時代の友人が研究していた玉ねぎの皮が向けるように、
割れていく花崗岩の球状風化や、
もっと過去の土石流で流れ、堆積した大きな石がたくさん見られました。

続いて、那智の滝に移動して、みんなで記念撮影。
那智の滝の麓近くでは、大正時代からの水力発電があるそうです。
先日の雨で水量も多く、迫力ある滝の様子に見入っていました。
那智の滝も、地質境界付近に位置し、滝となっている部分は
火成岩からなっているそうです。

大門坂の駐車場でお弁当を食べて尻剣谷に入りました。
ここも２３年の災害で土石流が発生した渓流です。

この渓流に分け入ると、



一度溶けて固まったような形の石や、人が掘ったと思われる穴、赤い土など
他のところとちがって、少し自然っぽくないものに出くわします。
というのも、那智地域では、昔から銅の採掘・精錬が行われていたからだそうです。
その名残が、落ちている石や地形に現れています。

さらに深く分け入ると、

石を積み上げてできた精錬所跡に出くわします。
熊野の産業の歴史の一端を垣間見ることができる貴重な場所です。

最後は、この渓流に作られた砂防堰堤を見学。

なかなか間近で見ることがない砂防堰堤。
鋼製のスリットで、大きな岩を止めることで土石流を止めることを狙ったものです。
おっきなジャングルジムみたいです。

近畿地方整備局の職員さんが丁寧に砂防堰堤の構造や工法、
効果や計画論を教えてくれました。
大人気で、質問攻めに遭われていました。
みなさん、砂防堰堤に関心が高いようでした。


前日の天候と打って変わって、とても穏やかで温かい一日でした。
普段は見えないところに入って、地形を詳しく見たり、自然の石を眺めたり
その違いに目を向けることで文化や歴史を感じることができました。
  

「夏休み子ども科学電話相談」

出張での車中、NHKラジオで電話相談の番宣を聞きました。
昔は現場に行く車中で、いつも楽しみに聞いていました。

子どもたちの質問に、真摯な姿勢で答える
専門家のみなさんの姿勢が印象的な番組。

質問の着想にも、「わからないことはわからない」と
はっきり答えたりもする、その回答もいつも素敵。

いままで聞いた中で、一番印象的な問題は、
「セミはなぜ７年土の中にいるの？」という趣旨の
ものに、
「７年いると確認されてるのは１種類だけで、
　他は類推だ」というもの。
「セミが本当に何年土のなかにいるのかを
　調べるのは、すべての土を調べて、まず
　他に卵がないことを確かめた上で、待つ
　必要がある」
だから、なかなか実証するのが大変。

しかも、
「アロエの木の下で育つセミは、
　４年くらいで成虫になる」というお話にはびっくり！

ちょっとうろおぼえなので、間違ってるかもしれませんが、
新鮮な発見に溢れる「夏休み子ども科学電話相談」は
永遠に続いて欲しいな。

http://www.nhk.or.jp/radiosp/kodomoq/  

国立天文台の天体シミュレーションソフト「Mitaka」が
日本経済新聞のHPに紹介されています。

和歌山サイエンスカフェでも何度か登場してくれた、このソフト。

すでに８０万回以上、ダウンロードされているそうです！

http://www.nikkei.com/article/DGXMZO89692440U5A720C1000000/?dg=1

夏休みの自由研究にもってこいですね。  

切目川ダムと初対面！

ダム本体工事中は立ち寄る機会もなかったので、
初対面でした。
直前のカーブを曲がるときは、ドキドキしましたが、
無事初対面を果たしました。

いままで完成予想図のイメージしかなかったので、
意外に優美な造形美に魅了されました。

たまたまオリフィスゲートから水が出ているところで、
下流面にあるオリフィスの上の覆いの尖ったフォルムが
また魅力的！

しかも、上流の流入口には、浮き上がった浮遊砂が
流線を描いていて、こちらも芸術的！

まさに水理学の教科書に出てくるような体験をできる
切目川ダムでした。