2018/11/09

JAGURSのNaN

JAGURSの計算結果には陸域(水が届かなかったメッシュ)のデータは含まれません.一方で,出力は2次元のメッシュデータなので何かの値が入っています.どうやら,それが -10000000000 のようです.これを無視してしまうと(というか,autoscaleで描画すると)以下のようになってしまいます.
これは,海に値が入っていないわけでも,陸が黒いわけでもなくて,データなしの -10000000000 が負に大きすぎて黒くなっているに過ぎません.なので,陸域の点を使って透明にしてあげて,色をつけてあげます.
右の2番めのアイコンを使うと,画像の点を指定して透明にする値を入力できます.
これをApplyすると,背景の地図が見えるようになります.
背景(地形図)が見えるようになりましたが,レンジが適切ではないので波が表示されないので,調整します.

これで見た目も良くなりました.

直交座標系での計算

JAGURSで直交座標系(CARTESIAN)で計算するには,Makefileで
    #CARTESIAN=ON
を有効にするために#のコメントアウトを外します.

従って,球座標(JAGURSのデフォルト)で作成したバイナリとは別になりますので時々で使い分ける場合はファイル名を変えるなり,ディレクトリを別にするなりしておく必要があります.

OpenTSUNAMI (杞憂プロジェクト)はぢめます.詳しくは後日.

2018/11/05

JAGURS on OSX (Brief summary in en)

I succeeded to compile and run JAGURS on Mac OSX in the following steps.

1. Install libraries via MacPorts

$ port install netcdf +dap +netcdf4
$ port install netcdf-fortran +gcc8
$ port install proj
$ port install fftw-3 +gfortran
$ port install szip

2. Install compiler from HPC Mac OS X compiler.

3. Download JAGURS from GitHub

4. Modify Makefile using CC=/usr/local/bin/gcc and FC=/usr/local/bin/gfortran. You should use -cpp option in gfortran command line to use preprocessor and -fopenmp option to use OpenMP directive. OR download Makefile.MacOSX.0500 and save as Makefile in src/ directory (diff patch).

日本語はこちら

gdalでnetCDFをFITSに変換する

GEBCOを紹介したら,天文な方からFITSで無い?と言われたので,変換してみた.

地震津波の分野では,NetCDFがよく使われています.もちろん,リモセンの時にも見てたので色々な分野で使われていると思います.一方,MacPortsのデフォルトのgdalではサポートしていないので,variants で指定する必要があります.

1年前にgdalでnetCDFを使うために,以下のようにしていました.
port install gdal +expat +netcdf +geos +hdf4 +hdf5 +xerces

さらにFITS も扱えるようにするため,さらにproj4を忘れていたようなので,以下のように作り直し(正確には新しいMacBookで新規インストール)ました.
port install gdal +expat +netcdf +geos +hdf4 +hdf5 +xerces +cfitsio +proj4

% gdalinfo GEBCO_2014_2D.nc | head -1
Driver: netCDF/Network Common Data Format
こんな感じで読めたのを確認しました.

% gdal_translate GEBCO_2014_2D.nc -of FITS GEBCO_2014_2D.fits
で変換.どうやって確認しようかなと思ったら,MacPotsに ds9 が入ってました!これは驚き.

ds9 @7.6b8_1 (science, graphics)
    SAOImage DS9 astronomical imaging and visualization application

% foreach i ({0..23})
dd if=GEBCO_2014_2D.fits iseek=$i bs=80 count=1 2>/dev/null && echo
end
SIMPLE  =                    T / file does conform to FITS standard
BITPIX  =                   16 / number of bits per data pixel
NAXIS   =                    2 / number of data axes
NAXIS1  =                43200 / length of data axis 1
NAXIS2  =                21600 / length of data axis 2
EXTEND  =                    T / FITS dataset may contain extensions
COMMENT   FITS (Flexible Image Transport System) format is defined in 'Astronomy
COMMENT   and Astrophysics', volume 376, page 359; bibcode: 2001A&A...376..359H
LONGSTRN= 'OGIP 1.0'           / The HEASARC Long String Convention may be used.
COMMENT   This FITS file may contain long string keyword values that are
COMMENT   continued over multiple keywords.  The HEASARC convention uses the &
COMMENT   character at the end of each substring which is then continued
COMMENT   on the next keyword which has the name CONTINUE.
HIERARCH LAT#AXIS = 'Y       '
HIERARCH LON#AXIS = 'X       '
END

ヘッダも出来てた.

が,,,180度回転してるな...

ImageMagickのconvert で -rotate 180 で回転は出来たものの,出力されたFITSは正の値しか持ってなかった.一度値を変換して,ほげほげすれば出来そうな気もするけど,一旦おしまい.

(追記 11.8)
ファイルがでか過ぎてマカリで開けなかったという連絡が来たので,半分に解像度を落とす.
 % grdsample  GEBCO_2014_2D.nc -GGEBCO_2014_2D_1min.nc -I1m
 % gdal_translate GEBCO_2014_2D_1min.nc -of FITS GEBCO_2014_2D_1min.fits

これで,サイズは1/4と思ったら,Float32になっちゃってて半分にしかならなかった...


2018/10/27

ガウス分布初期波源計算

JAGURSでは,Gaussian分布を海底地殻変動(初期水位)として計算することが出来ます.が,,,マニュアルには球座標使用の場合しか書かれていないので少し補足です.

tsun.par に
init_disp_gaussian =1
とすることで使えます.で,gaussianというファイル名でパラメータを用意します.が,単位に注意が必要です.

球座標の場合:
&gaussian h0=波高 lon_o=経度[度] lat_o=緯度[度] L=半値幅[km] /

直交座標の場合:
&gaussian h0=波高 lon_o=X座標 lat_o=Y座標 L=半値幅[座標系の単位] /

のように,m(メートル)が単位になっているUTM等を利用している場合,Lで指定する半値幅はmが単位になります.一方で,球座標の場合,座標は「度」単位,半値幅はkm単位で固定です.

なお,波高に負を指定することで凹んだ初期水位を表現することも可能です.実装は mod_init_disp_gaussian.f90 に入っていますので必要な関数に書き直しちゃうのも簡単かなと.

JAGURS V0407以前では,直交座標の場合には波高(h0)が指定できない(1 mに固定)ので注意が必要です.

2018/10/23

JAGURS on Mac 2

だいぶ時間も経ったし,JAGURSも新しくなったので手順の見直しも兼ねて再整理しておきます.以前の記事はこちら

基本的な環境は以前と同じ.



$ port install netcdf +dap +netcdf4
$ port install netcdf-fortran +gcc8
$ port install proj
$ port install fftw-3 +gfortran
$ port install szip

修正後Makefile
Makefileの差分

というわけで,v0501になっても,手順の変更はなしでした.

せっかくなのでバイナリも配ろうかと思いましたが,どうやらMac OSXではstatic linkが出来ないとのこと...諦めます(手順を教えて頂ける方がいらっしゃいましたらやります).これで出来る??? https://code.i-harness.com/ja/q/ce413

2018/10/04

QGIS 3.x

新しもの好きではないので積極的に何かを食べてみようとは思わないのですが,QGISはそろそろ 3.x にしておいた方が無難かなと思ったのでやってみた.

https://www.qgis.org/ja/site/forusers/download.html
インストーラーも整備されているので何のこともないのですが,Pythonがやっかいでした.日本人が書いていると思うのですが,分かりにくい日本語...Python 3.6 が必要なのですが(私はMacPorts で3.6を入れている),このQGIS3は,https://www.python.org/ 本家で配布されているバイナリじゃなきゃ動きませんよ.ということでした.こういう基盤ツールが複数入るのは後々の罠になるので嫌なのですが,そうも言っていられないので,最新 3.6.6 はあるものの, 3.6.5 で動作確認したと書かれているので(どうせQGIS以外では使わないので) 3.6.5 のバイナリをインストールしました.それ以外は,書いてある通りでOK.

最低限の整備として,背景地図.

Open Street Map の場合:
XYZ Tiles -> OpenStreetMap -> 選択したレイヤをキャンパスに追加 で終わり...


地理院地図の場合:
XYZ Tiles -> 新しい接続 に
標準地図 https://cyberjapandata.gsi.go.jp/xyz/std/{z}/{x}/{y}.png
淡色地図 https://cyberjapandata.gsi.go.jp/xyz/pale/{z}/{x}/{y}.png
を追加.
参考)QGIS3で地理院地図を追加する方法
地理院地図)https://maps.gsi.go.jp/development/ichiran.html

困ったな.津波の正負の波を入れるのに,赤ー青をよく使っていましたが,「反転」のチェックボックスがなくなった...困ったな...と思っていたら,カラーテーブルを選ぶプルダウンの中に,"invert color" があった...


あまりにもQGIS 3.x の情報が少ないので,そういう検索に掛かり易いようにしておこうかなと思います.

2018/07/27

QGISでライン表示

x,yが書かれたファイルをQGISで読み込む(レイヤ→レイヤの追加→デリミティッドテキストレイヤの追加)と,ポイントデータとして読み込まれます.


これを線で結ぶには,[Points to path] という操作が必要だそうです.
プロセッシングからツールボックスを出して,QGISジオアルゴリズム→ベクタ作成ツール→ポイントをパスへをダブルクリックすると,下のようなメニューが開きます.
この時,グループフィールドというのが指定できないと動いてくれないようです.例えば四角形1つだったとしても,x,y に加えて,例えば全部 0 とかの数字が入ったカラムを用意しておいて,それをグループフィールドとして指定しないとうまく動いてくれないようです.
ポリゴンにしたい場合は,さらに lines to polygons (ベクタ→ジオメトリツール→ラインをポリゴンへ)で変換出来ます.
それぞれのステップ(points to lines, lines to polygons)で,元のレイヤーはそのまま残るので,必要に応じて消すなりなんなり.


QGISで点つなぎを参考にさせて頂きました.

2018/05/29

2018/03/24

Keybord settings

Macで新しいキーボードを手に入れてから(Magick Keyboard の雰囲気が随分変わってしまったので古いマシンも書い直した)最初にやるのはこれ↓

Caps Lock を ^Control に割り当てる.

キーボードを英語配列にしたのは,Happy Hacking を使い始めた頃かな…その前が思い出せない...学校や自宅のPC-9801は当然日本語配列だったはずだし,その後, Mac LC 475を中古で買った頃も何も考えてなかったはずだから日本語配列だったと思われる.たぶん,その後のデスクトップはLinux boxで作る世界に入ったから英語にしたのかも.その原因は,研究で使ってた Solaris 端末の影響か?その後,VAIO noteで英語配列が選べる時に英語にしたことで,完全US配列化した.それが,505だったのか,SRだったのか思い出せない.移行,ずっとUS配列で来たので "A" の左に Ctrl が居てくれないと手がついて行けないのよね.

一方で,キーピッチやキーストロークには一切の趣向は持っていません.どんなキーボードでも基本的には指が合わせてくれます.

おまけは余ったMagick keyboardをiPad miniに繋いでみたの図


2018/03/19

VMWare windowsの画面が真っ黒

なかなか困りました.
VMWare fusionでWindowsを動かしていますが,ひょんなことから(いつからか不明)画面が真っ黒なままなかなか操作できない.マウスを持っていくと,ぐるぐるマークが回っているので何かは動いている模様.仮想マシンメニューで再起動しようがシャットダウンしてから起動しようが,Windowsの起動画面っぽいところまでは見えてるけど,すぐに真っ黒.常時真っ黒かというと,一瞬ちらっとログインプロンプトが見えて利する状態.

で,何で解決出来たかというと,設定ウィンドウの中の「ディスプレイ」で,「Retinaディスプレイのフル解像度を使用」のチェックを外した所,無事に表示されました.その後チェックをONにしても動いているようなので,このオプション自体が悪さをしているといういわけでもなさそうです(これまでも使えてたし).何らかのリセットがかかって無事に認識されたのでしょう.


2018/03/17

Mewのマスターパスワード(Google2段階認証)

Mewでマスターパスワードを設定すれば,個々のIMAPやPOPのパスワードを覚えなくても良いというのは知ってはいたのですが,なんとなく使っていませんでした.

が,Googleの認証を2段階認証にしてしまうと,アプリパスワードが必要になって,こいつがとてもじゃないけど覚えられたものじゃない...8桁の不規則パスワードならなんとか覚えてましたが,12桁は無理.「ほとんどの場合、アプリ パスワードはアプリや端末ごとに一度入力するだけなので、覚える必要はありません。」とか手順書に書かれている通り,パスワードはアプリが記憶するものともう世の中(Google)が決めてしまっているんですね.Mewで起動するたびに入力する時代は終わったか.

というわけで,mewでもパスワードを覚えてもらうために,mew-use-master-passwdを設定しました.GnuPGは,MacPortsでgnupg2を入れました.Mew側はGnuPGのバージョンは複数ケアしてくれてるようです.
% port installed gnupg2
The following ports are currently installed:
  gnupg2 @2.2.5_0+pinentry_mac (active)


2018/03/14

天文学辞典

日本天文学会から天文学辞典βがウェブで公開されていました.公開当初はいろいろな不備があると思われますが,,,ということですが,網羅性は凄いですね.

2018/01/09

地形モデル

遠地津波計算をするためにGEBCOのデータを変換する手順は以前書きましたので,今回は近地津波の計算をするための直交座標系での地形データの作成手順を書いておきます.

元データは内閣府の南海トラフの巨大地震モデル検討会で作成された地形データがG空間情報センターからダウンロード出来ます.無償にはなっていますが,ダウンロードするためにはアカウントを作成する必要があります.地形データ以外の粗度や堤防,初期水位データ等もダウンロード出来ます..

例として,JAGURS で使える形式に変換してみます.
手順は大きく分けで3つ.(1)まずはデータをダウンロードして,(2)まともなGISフォーマットに変換して,(3)最後に必要なフォーマットへの変換をします.
地形データは,ここからダウンロードできます.利用規約はもちろん確認が必要ですが,「計算入力データ(地形,粗度,堤防データ)の説明」を確認して下さい.ここに何系がどの領域に対応しているか書かれています.例えば四国を計算する場合は第04系になるのでその地形データをダウンロードします.ダウンロードしたデータにはメタ情報が何も書かれていないので,別途「計算範囲設定」を確認する必要があります.下記の ncols, nrows, xllcorner, yllcorner の値は,エクセルのメッシュ個数(X方向,Y方向),南西端の位置(X座標,Y座標)に対応します.

データは 10f8.2 で記述してあるということなので,Fortranで読み込むのは簡単ですが, f8.2 の場合,標高(陸域)が1000mを超える場合に,符号も含めると -1000.00 となり8桁になってしまうため,前の数字とスペースなしで連結されてしまいます.Fortranなら読めるでしょって言ってしまえばそれまでなのですが,一般的なアスキーフォーマットと考えた場合,スペースで区切られていないと使えません.というわけで, 10f8.2read して, 10f9.2write します.さらに,アスキーフォーマットとしてGISデータになるようにヘッダを付与します.プログラム(reformat.f)も掲載しておきます.
     ncols 720
     nrows 540
     xllcorner -730200
     yllcorner -857700
     cellsize 2430
     NODATA_value -9999

ヘッダの中身は,さほど難しくなく,ncolsはX方向のメッシュ数,nrowsはY方向のメッシュ数,xllcorner は南西端のX座標,yllcorner は南西端のY座標,cellsize は解像度(上記例は 2430m),NODATA_value は不要ですがお決まりの -9999 を設定しています.

reformat.f:
      PROGRAM reformat
      real*8 data(10)
      character(len=128) :: fname
      read(*,'(a)') fname
      open(10,file=fname)
      do
         read(10,'(10f8.2)',end=100) data
         write(*,'(10f9.2)') data
      end do
 100  END

作成したアスキーファイル(bathy.SD01.asc)をJAGURS用のcf形式に変換します.
% grdreformat bathy.SD01.asc=gd bathy.SD01.grd=cf
これで,ネスティングしない場合は計算できるようになりました.ネスティングする場合は,メッシュの細かい調整が必要なので後日掲載します(一応,手元では動くデータは作れましたがこれで良いか確認します).

日本海講側は,中央防災会議のものが公開されています(古いです).