玉川研究室で、CReSS を使うには。その２ 地形データ篇
最初に戻る
研究室には、２万５千分の１地形図から作られた 50m メッシュ標高データと
呼ばれるデータがある。(例えば、192.168.1.12 の
/home/tama/Map や /work/Map 以下や、学生部屋の本棚の中)
これを変換して、CReSS で使えるようにする。
この作業は、計算を実行する並列計算機ではなく、
たとえば、192.168.1.12 など。他のPCで行なって下さい。
並列計算機は、あくまでも計算用として使いたいので、GMT などのソフトを
整備していません。
まず、計算領域を決定する。そして、その範囲をカバーするように地図を選択する。
この時、 国土地理院発行地図一覧図 (日本地図センター) を、参照しつつ、
CDROM にある zumei.csv を調べるか、
あるいは、CDROMについている、標準地域メッシュコード一覧図を見て、使用するファイル名を
求める。
以下、Unix 系OS上での作業を仮定して書かれている。
Windows からなら、ダウンロードしたファイルを、GMT や 地図データが入っていて
なおかつ Unix系OS (Linux 等) の動作しているPCに(WinSCP などで)転送して、
ターミナル(putty など)上で操作する。
個々の地図データファイルの名前を指定して作成する方法
(実は、ちょっと領域が大きくなると指定するのがえらく面倒臭くなる)
ブロック毎に指定する楽な方法はこちら
(こちらを推奨するのですぐ下を読まずに、クリックしてずっと下へ飛んで行ってください。)
処理用に玉川が書いたソフト ForCReSS_ter.c を使ってデータを変換する。
このプログラムは、「データファイルの並びの前に、
X方向のファイル数 Y方向のファイル数 データファイルのあるディレクトリ 出力ファイル名
を書いた１行を加えたコマンドファイル」を作って、標準入力から食わせることで
処理を行う。
(この際、データファイル名が sea であれば、海面の高度 (とりあえず -0.01とした) で
埋める。)
たとえば、高山(543712.mem)を中心に東西南北10個ずつファイルを取り、縦横21個ずつ使うとすると
(実は、後の計算の実行のところで、これでは領域が足りない事が分かった。
縦横２４個とって(右下に多い)作成したデータを実際の計算に使った。 ここを参照 )
21 21 /home/tama/Map/Alldata/ Takayama21.ter
sea sea sea sea sea sea 553636.mem 553637.mem 553730.mem sea sea 553733.mem 553734.mem 553735.mem 553736.mem 553737.mem 553830.mem 553831.mem 553832.mem 553833.mem 553834.mem
sea sea sea sea sea sea 553626.mem 553627.mem 553720.mem sea sea 553723.mem 553724.mem 553725.mem 553726.mem 553727.mem 553820.mem 553821.mem 553822.mem 553823.mem 553824.mem
sea sea sea sea sea 553615.mem 553616.mem 553617.mem 553710.mem 553711.mem 553712.mem 553713.mem 553714.mem 553715.mem 553716.mem 553717.mem 553810.mem 553811.mem 553812.mem 553813.mem 553814.mem
sea sea sea sea sea 553605.mem 553606.mem 553607.mem 553700.mem 553701.mem 553702.mem 553703.mem 553704.mem 553705.mem 553706.mem 553707.mem 553800.mem 553801.mem 553802.mem 553803.mem 553804.mem
sea sea sea sea 543674.mem 543675.mem 543676.mem 543677.mem 543770.mem 543771.mem 543772.mem 543773.mem 543774.mem 543775.mem 543776.mem 543777.mem 543870.mem 543871.mem 543872.mem 543873.mem 543874.mem
sea sea sea 543663.mem 543664.mem 543665.mem 543666.mem 543667.mem 543760.mem 543761.mem 543762.mem 543763.mem 543764.mem 543765.mem 543766.mem 543767.mem 543860.mem 543861.mem 543862.mem 543863.mem 543864.mem
sea sea sea 543653.mem 543654.mem 543655.mem 543656.mem 543657.mem 543750.mem 543751.mem 543752.mem 543753.mem 543754.mem 543755.mem 543756.mem 543757.mem 543850.mem 543851.mem 543852.mem 543853.mem 543854.mem
sea sea 543642.mem 543643.mem 543644.mem 543645.mem 543646.mem 543647.mem 543740.mem 543741.mem 543742.mem 543743.mem 543744.mem 543745.mem 543746.mem 543747.mem 543840.mem 543841.mem 543842.mem 543843.mem 543844.mem
543630.mem 543631.mem 543632.mem 543633.mem 543634.mem 543635.mem 543636.mem 543637.mem 543730.mem 543731.mem 543732.mem 543733.mem 543734.mem 543735.mem 543736.mem 543737.mem 543830.mem 543831.mem 543832.mem 543833.mem 543834.mem
543620.mem 543621.mem 543622.mem 543623.mem 543624.mem 543625.mem 543626.mem 543627.mem 543720.mem 543721.mem 543722.mem 543723.mem 543724.mem 543725.mem 543726.mem 543727.mem 543820.mem 543821.mem 543822.mem 543823.mem 543824.mem
543610.mem 543611.mem 543612.mem 543613.mem 543614.mem 543615.mem 543616.mem 543617.mem 543710.mem 543711.mem 543712.mem 543713.mem 543714.mem 543715.mem 543716.mem 543717.mem 543810.mem 543811.mem 543812.mem 543813.mem 543814.mem
543600.mem 543601.mem 543602.mem 543603.mem 543604.mem 543605.mem 543606.mem 543607.mem 543700.mem 543701.mem 543702.mem 543703.mem 543704.mem 543705.mem 543706.mem 543707.mem 543800.mem 543801.mem 543802.mem 543803.mem 543804.mem
533670.mem 533671.mem 533672.mem 533673.mem 533674.mem 533675.mem 533676.mem 533677.mem 533770.mem 533771.mem 533772.mem 533773.mem 533774.mem 533775.mem 533776.mem 533777.mem 533870.mem 533871.mem 533872.mem 533873.mem 533874.mem
533660.mem 533661.mem 533662.mem 533663.mem 533664.mem 533665.mem 533666.mem 533667.mem 533760.mem 533761.mem 533762.mem 533763.mem 533764.mem 533765.mem 533766.mem 533767.mem 533860.mem 533861.mem 533862.mem 533863.mem 533864.mem
533650.mem 533651.mem 533652.mem 533653.mem 533654.mem 533655.mem 533656.mem 533657.mem 533750.mem 533751.mem 533752.mem 533753.mem 533754.mem 533755.mem 533756.mem 533757.mem 533850.mem 533851.mem 533852.mem 533853.mem 533854.mem
533640.mem 533641.mem 543642.mem 533643.mem 533644.mem 533645.mem 533646.mem 533647.mem 533740.mem 533741.mem 533742.mem 533743.mem 533744.mem 533745.mem 533746.mem 533747.mem 533840.mem 533841.mem 533842.mem 533843.mem 533844.mem
533630.mem 533631.mem 533632.mem 533633.mem 533634.mem 533635.mem 533636.mem 533637.mem 533730.mem 533731.mem 533732.mem 533733.mem 533734.mem 533735.mem 533736.mem 533737.mem 533830.mem 533831.mem 533832.mem 533833.mem 533834.mem
533620.mem 533621.mem 533622.mem 533623.mem 533624.mem 533625.mem 533626.mem 533627.mem 533720.mem 533721.mem 533722.mem 533723.mem 533724.mem 533725.mem 533726.mem 533727.mem 533820.mem 533821.mem 533822.mem 533823.mem 533824.mem
533610.mem 533611.mem 533612.mem 533613.mem 533614.mem 533615.mem 533616.mem 533617.mem 533710.mem 533711.mem 533712.mem 533713.mem 533714.mem 533715.mem 533716.mem 533717.mem 533810.mem 533811.mem 533812.mem 533813.mem 533814.mem
533600.mem 533601.mem 533602.mem 533603.mem 533604.mem 533605.mem 533606.mem 533607.mem 533700.mem 533701.mem 533702.mem 533703.mem 533704.mem 533705.mem 533706.mem 533707.mem 533800.mem 533801.mem 533802.mem 533803.mem 533804.mem
523670.mem 523671.mem 523672.mem 523673.mem 523674.mem 523675.mem 523676.mem 523677.mem 523770.mem 523771.mem 523772.mem 523773.mem 523774.mem 523775.mem 523776.mem 523777.mem 523870.mem 523871.mem 523872.mem 523873.mem 523874.mem
というような内容のファイルを作って(ここでは SAMPLEDATA とする)、
cc -o ForCReSS_ter ForCReSS_ter.c
unlimit
./ForCReSS_ter < SAMPLEDATA |& tee LOG
とすると、画面にメッセージが出力され(標準エラー出力)、
同じものが LOG にも
書かれる。(コマンド tee の機能、詳しくは man tee せよ )
こんな情報が末尾にでるはず
ここで使う C 言語で書かれたプログラム ForCReSS_ter.c は、ここからダウンロードする 。
また、２行目のコマンド unlimit は、ユーザに割り当てられているメモリ使用量の
上限を外すコマンドで、今回のように大きなメモリを使用するソフトを起動する際に
用いると良い。
for CReSS information ( 1 , 1 ) (center of lowest left grid) =
( 135.997451 , 35.253403 ) in deg.
data number 4200 x 4200
dx 0.000625 dy 0.000417 deg.
Upper Right = ( 138.622451, 37.004803)
その結果、4200 x 4200 の標高データが、南西の隅から東へ、
次に南から２番目の列のデータがやはり西から東へ、と一番北の列まで続けて、
4byte浮動小数点データで出力されたデータが、Takayama21.ter に出力される。
(もしも、データが大きすぎるエラーになったら、プログラムの冒頭部の #define を
いじるべし。また実行時の Segmentation Error の類は、配列確保の問題なので、
unlimit を実行してから試してみる(上述)などの対応をすべし)
ここで、その出力結果を確認してみる。
簡単には od コマンドを使って od -f Takayama11.ter | more とでもすれば、
充分数字が確認できるが、
ここでは、その結果を確認するためのソフト、check_map.c もあるので、
それを使って確認をしてみる。
check_map.c はここ
cc check_map.c -o check_map
./check_map
xnum ynum =
4200 4200
xpos ypos (lowest left in degree)=
135.997451 35.253403
dx dy in deg.
0.000625 0.000417
input file name=Takayama21.ter
output file name=Takayama21.txt
(ここで、途中、対話的にやりとりするので、上の例には、ソフトの出力と
キーボードからの入力が混在している)
その結果できた Takayama21.txt は、X Y Z のデータの並ぶテキストデータであるので、
GMT を使った以下のようなスクリプトを実行すると、
チェック用のスクリプト(GMT)はここ
set ROPT = 135.997451/138.622451/35.253403/37.004803
set TXTDATA = Takayama21.txt
set GRDDATA = Takayama21.grd
set FIG = Takayama21.ps
set CPTFILE = CPT.cpt
set TOPT = 0/3000/500
set JOPT = M10
gmtset BASEMAP_TYPE plain
xyz2grd ${TXTDATA} -G${GRDDATA} -I0.000625/0.000417 -R${ROPT}
makecpt -Crainbow -T${TOPT} > ${CPTFILE}
psbasemap -Ba1f0.2/a1f0.2:."Check for Takayama 21x21": -J${JOPT} -R${ROPT} -K > ${FIG}
grdimage ${GRDDATA} -R -C${CPTFILE} -J${JOPT} -O -K >> ${FIG}
grdcontour ${GRDDATA} -R -C${CPTFILE} -J${JOPT} -O >> ${FIG}
のような図が描け、まぁちゃんと変換できているんだろうなぁと確認できる。
(実際には、この後 convert -rotate 90 Takayama21.ps Takayama21.png として作成した図が、
ここにある)
(2007.11)
上の変換プログラム ForCReSS_ter.c は、ファイルの指定のために
多くのファイル名を書かねばならないので、とても面倒である。
そこで、４桁の地図名 (5437 とか言うようなもの, 多分 第一次地域区画と言う) 単位で指定し、
更に、左下と右上だけを指定して読み込むように改造したものを作成した。
もし、その地図データが無ければ、海だということにして、上と同様の
海の処理をしてある。
プログラムは、ここ ForCReSS_ter_2.c
ダウンロードした後、
cc ForCReSS_ter_2.c
./a.out 地図データのあるディレクトリ名 左下図面番号 右上図面番号 出力ファイル
と指定するだけで、出力ファイルに CReSS の地形処理プログラム terrain が
期待する左下からべたっとバイナリ書きの標高データが出来上がる。
例えば、
./a.out /work/Map/Alldata 5034 5135 test-small.bin
: : :
for CReSS information ( 1 , 1 ) (center of lowest left grid) =
( 133.997646 , 33.336931 ) in deg.
data number 3200 x 3200
dx 0.000625 dy 0.000417 deg.
Upper Right = ( 135.997646, 34.671331)
の様に実行した場合、紀伊水道の南側の標高データが作成される。
チェックの様子は、上の例と同様であるが、以下のようなものである。
cc check_map.c -o check_map
./check_map
xnum ynum =
3200 3200
xpos ypos (lowest left in degree)=
133.997646 33.336931
dx dy in deg.
0.000625 0.000417
input file name=test-small.bin
output file name=test-small.txt
できたテキストファイル test-small.txt を
set ROPT = 133.997646/135.997646/33.336931/34.671331
set IOPT = 0.000625/0.000417
set TXTDATA = test-small.txt
set GRDDATA = tmp.grd
set FIG = test-small.ps
set CPTFILE = CPT.cpt
set TOPT = -7000/7000/500
set JOPT = M10
gmtset BASEMAP_TYPE plain
xyz2grd ${TXTDATA} -G${GRDDATA} -I${IOPT} -R${ROPT}
makecpt -Ctopo -T${TOPT} > ${CPTFILE}
psbasemap -Ba1f0.2/a1f0.2:."Check": -J${JOPT} -R${ROPT} -K > ${FIG}
grdimage ${GRDDATA} -R -C${CPTFILE} -J${JOPT} -O -K >> ${FIG}
grdcontour ${GRDDATA} -R -C${CPTFILE} -J${JOPT} -O >> ${FIG}
の様なスクリプトで処理すると
のような図ができ、ちゃんと動いている事が
分かる。
とは言うものの、大きな領域で作成すると、図化することもできない
(もしかしたら、CReSS の terrain すら動かないかも) 巨大なデータに
なっていることがある。
この様な場合、実際には 50m などという解像度が
不要なことが普通である。(必要なら処理できるコンピュータ環境で仕事しているはず)。
そこで、データを m × n 個で平均する処理を行うことにする。
そのためのソフトは、 data_reduce.c である。
使い方は、
cc data_reduce.c -o data_reduce
./data_reduce filename m n m_ave n_ave outputfile
ここで、
- filename が、ForCReSS_ter_2.c などで作られた CReSS 用 binary 標高データファイル。
- m, n がそのファイルの東西、南北方向のグリッド数。
- m_ave n_ave が、それぞれの方向の平均するグリッド数。
- outputfile が、平均した結果を出力するファイルの名前。
である。(2007 Nov. 追加)
ソフトのダウンロードは、ここ
GMT に関しては
とりあえず研究室の関連ページ参照