いい感じの風分布図の書き方
こんにちは。おこめまんです。
今回は
「地上風分布・風鉛直断面を風矢印を用いて書きたい。」
と思います。
GMTを使用し、grdvectorで矢印を書きます。
風を矢印で書く場合、風速が弱い場合、「陥没」する場合があります。
見栄えが悪いので、凹まないようにしましょう。
①すべての矢印の長さをそろえてしまう。
-Qオプションと-Sオプションで調整します。
-Q<arrowwidth>/<headlength>/<headwidth> と-S<矢印の長さ>
-Sオプションで矢印の長さを指定すると、図に描画される矢印の長さはすべて同じになります。
②すべての矢印の長さと大きさを変えつつうまく描く。
こちらも-Qオプションと-Sオプションで調整します。
-Q<arrowwidth>/<headlength>/<headwidth>n<headlength>と-S<矢印の長さ>
-Qオプションで矢印の形を指定するのと同時に、nをくっつけて、最小の長さを指定します。n<ココ！>の数字をheadlengthと同じ値にすることによって、凹むことはなくなります。※-Sオプションを付けることが必須です。
以上、ご参考ください。
時系列データをGMTで描画する
みなさんこんにちは。おこめまんです。
今回は「時系列グラフ」をGMTで書いてみます。
はじめに
みなさんがこのブログにたどり着いたということは、「時系列ってできるの？」とか「時刻上手くGMTで書けないんだが…」と思っていることだと思います。
気象データを時系列グラフにすることで、時間変化を追うことが出来たり、新しい発見があったりするものです。ぜひ、書けるようにしましょう。
さて、みなさんが書きたい図はこんな感じでしょうか。
この図は、ある地点のアメダスデータの風向・風速の時系列変化を表しています。
風向(右軸)はシンボル、風速(左軸)は実線です。
肝心の「時刻」(横軸)に注目すると、3/27 00:00 ~ 3/30 00:00となっています。
おそらく、皆さんがお困りなのは、この時刻をどう書かせたのかではないでしょうか。
シェルスクリプト
※fig_timeseries.shというシェルスクリプトを作成して、その中身を以下の通りとしました。赤色の部分は各々変更しましょう。
#
#!/bin/bash
#
######################################################################
### 作図日時指定 ###
SYEAR=2000
SMONTH=03
SDAY=27
SHOUR=00
SMINUTE=00
SSECOND=00
EYEAR=2000
EMONTH=03
EDAY=30
EHOUR=00
EMINUTE=00
ESECOND=00
### AMeDASデータ ###
ifile='amedas_200003.txt'
ifn=./${ifile}
### 風速の下限/上限指定 ###
minws=0
maxws=20
### 風向上限/下限指定 ###
minwd=0
maxwd=360
######################################################################
### 日時の結合 ###
starttime=${SYEAR}-${SMONTH}-${SDAY}T${SHOUR}:${SMINUTE}:${SSECOND}
endtime=${SYEAR}-${EMONTH}-${EDAY}T${EHOUR}:${EMINUTE}:${ESECOND}
### 作成する図＜風＞ ###
ofile='wind_'$SYEAR$SMONTH$SDAY'.ps'
mkdir -p ./figures/
ofn=./figures/${ofile}
### GMTのデフォルト変更 ###
#gmtdefaults -D
gmtset ANNOT_FONT_SIZE_PRIMARY 28p
gmtset HEADER_FONT_SIZE 35p
gmtset LABEL_FONT_SIZE 30p
gmtset INPUT_DATE_FORMAT yyyy-mm-dd
gmtset INPUT_CLOCK_FORMAT hh:mm:ss
gmtset PLOT_DATE_FORMAT mm-dd
gmtset PLOT_CLOCK_FORMAT hh:mm
############ 作図 ############
#####1軸＜風速＞ #####
psbasemap -R$starttime/$endtime/$minws/$maxws -Jx2.5i/0.5 -Bs1D::/:: -Ba12h:"JST":/2.5:"Wind Speed [m/s]":WS -X5 -Y5 -K > $ofn
### 風速時系列描画 ###
awk '$9 > 0.0 {print $7" "$9}' $ifn | \
psxy -Jx -R -W8 -K -O >> $ofn
### 2軸＜風向＞ ###
psbasemap -R$starttime/$endtime/$minwd/$maxwd -Jx2.5i/0.02777 -Bs1D::/:: -Ba12h:"JST":/60:"Wind Degree [\260]":ES -K -O >> $ofn
### 風向時系列描画 ###
awk '$10 > 0.0 {print $7" "$10}' $ifn | \
psxy -Jx -R -Sc0.5 -G100 -N -O >> $ofn
ps2raster -P -A -Tg $ofn
シェルスクリプトの説明
※大事なところを説明します！
### 作図日時指定 ###
SYEAR=2000
SMONTH=03
SDAY=27
SHOUR=00
SMINUTE=00
SSECOND=00
EYEAR=2000
EMONTH=03
EDAY=30
EHOUR=00
EMINUTE=00
ESECOND=00
→いきなりですが、この部分で描画したい範囲での時刻を指定します。今回は2000/3/27 00:00:00～2000/3/30 00:00:00の時系列グラフを作成します。いろんな日付の作図を行う際にはいちいち指定するのは大変なので、以前書いたこちらを参照して適宜書き換えましょう。
*****************************************************************************************************
### AMeDASデータ ###
ifile='amedas_200003.txt'
ifn=./${ifile}
→ここで作図するinputfileを指定しました。中身ですが、2000/3/27 00:00:00～2000/3/30 00:00:00のデータが1行ずつ記述されているものとします。シェルスクリプトと同じディレクトリに置きましょう。
アメダスデータ10分間値の取得方法はこちらを参照してください。
今回は、左から「年、月、日、時、分、秒、日時を結合したやつ、気温、風速、風向、降水、日照」となっているテキストデータ(コンマ区切り)を用意しました。
先ほどの図はこのG列、I列、J列を使用して作図しました。
はい。でました。「日時を結合したやつ」！！！！
ここがとっても大事です。詳しくみていきましょう。
1行目に注目。
2000-03-27T00:00:00
これは2000/3/27 00:00:00を表しています。注意したいのは間をつないでいる「T」です。この形にそろえましょう。
たとえば、Excelを使用して整形する場合…
整形前のifileには「日時を結合したやつ」なんて入っていないので付け足します。
①新しく列を作る(今回はG列を追加しました。)
②G1に「=DATE(A1,B1,C1)+TIME(D1,E1,F1)」を書き込む
→ここで、日時を結合しています。
③G列についてホーム→数値→その他の表示形式→ユーザー定義から「yyyy-mm-ddThh:mm:ss」のように変更する
→これでTを含んだかたちで表示されるようになります。
*****************************************************************************************************
### 日時の結合 ###
starttime=${SYEAR}-${SMONTH}-${SDAY}T${SHOUR}:${SMINUTE}:${SSECOND}
endtime=${SYEAR}-${EMONTH}-${EDAY}T${EHOUR}:${EMINUTE}:${ESECOND}
→ここで、psbasemapの横軸範囲指定を行っています。「作図日時指定」を「T」を含んだ形に変換します。ここでも「T」を含んだ形にしましょう。
*****************************************************************************************************
### GMTのデフォルト変更 ###
#gmtdefaults -D
gmtset ANNOT_FONT_SIZE_PRIMARY 28p
gmtset HEADER_FONT_SIZE 35p
gmtset LABEL_FONT_SIZE 30p
gmtset INPUT_DATE_FORMAT yyyy-mm-dd
gmtset INPUT_CLOCK_FORMAT hh:mm:ss
gmtset PLOT_DATE_FORMAT mm-dd
gmtset PLOT_CLOCK_FORMAT hh:mm
→特に橙色で示したところは必ずシェルスクリプト中に書きましょう。
*****************************************************************************************************
#####1軸＜風速＞ #####
psbasemap -R$starttime/$endtime/$minws/$maxws -Jx7.0i/0.5 -Bs1D::/:: -Ba12h:"JST":/2.5:"Wind Speed [m/s]":WS -X5 -Y5 -K > $ofn
→ベースマップの指定です。大事そうなところだけ。
-R$starttime/$endtime/$minws/$maxwsは図の範囲を指定しています。$starttime/$endtimeで横軸(時間)、$minws/$maxwsで縦軸(風速)です。
-Bs1D::/::のs1Dは横軸の月日を表示させるためのものです。この設定だと1日ごとの表記になります。
-Ba12h:"JST":/2.5:"Wind Speed [m/s]":WSのa12hは横軸の目盛を振る間隔です。12時間ごとの場合はこの表記になります。
*****************************************************************************************************
その他のところは、図の体裁にかかわる部分です。各自好きに変えましょう。
おわりに
時系列グラフを作成する場合、「T」を付けて日時を結合することがポイントとなります。ifile、psbasemap両方。
ifleさえ準備できれば、シェルスクリプト完コピで作図されるはずです。お試しください。
複数ファイルを次々に読み込みGMTで作図したい
みなさんこんにちは。おこめまんです。
今回はGMTでの作図を行う際のTipsを扱います。
はじめに
この記事を読んでいる皆さんは気象データをGMTを用いて図化することが多いかと思います。
このとき「日時のみが異なるifile(描画する元データの入っているファイル)を繰り返して(ループをかけて)シェルスクリプトで指定してGMTで図化したいなあ」と思っていませんか。
描きたい図が少ない場合は気にならないかもしれませんが、大量の図を描く場合は大変です。シェルスクリプトでifileをいちいち指定し直して描画するのでは時間がかかってしまいますよね。ぜひとも自動化したいものです。
今回は、「日時のリストを作成してGMTでつぎつぎに描画する方法」を伝授しましょう。
今回ループで開きたいifile
2018-01-05_00-00-00.txt
2018-05-20_01-00-00.txt
2018-03-31_23-00-00.txt
例として、ファイルの名前が年、月、日、時、分、秒から構成されるifileを扱います。
ただし、各ファイルのデータの構成は同様のものとします。
シェルスクリプト
※適当なシェルスクリプト「example.sh」の中身を以下の通りとします。
----------------------------------------------------------------------------------------------------------
#!/bin/bash
cat << EOF > date.txt
2018 01 05 00 00 00
2018 05 20 01 00 00
2018 03 31 23 00 00
EOF
while read YEAR MONTH DAY HOUR MINUTE SECOND
do
### 作図する日時を表示 ###
echo "draw date ---> " $YEAR $MONTH $DAY $HOUR $MINUTE $SECOND
### 図化スクリプト ###
ifile=$YEAR'-'$MONTH'-'$DAY'_'$HOUR'-'$MINUTE'-'$SECOND'.txt'
psbasemap ~~~~~~~~~
(省略)
ps2raster ~~~~(convert~~~~~)
done < date.txt
----------------------------------------------------------------------------------------------------------
シェルスクリプトの説明
cat << EOF > date.txt
2018 01 05 00 00 00
2018 05 20 01 00 00
2018 03 31 23 00 00
EOF
日時を各項ごとに分けてdate.txtに書き込みます。いちいちExcelでこのリストを作成することは面倒ですので、シェルスクリプトで作成しています。
cat << EOF > date.txt
(ここに1行ずつ日時を記述します。このとき、各項をスペースで区切りましょう。)
EOF
注意！！
2018 03 31 23 00 00
(ここ！空行ダメ！)
EOF
日時とEOFの間に空行を作らないようにしましょう。
-------------------------------------------------------------------------------------------------
while read YEAR MONTH DAY HOUR MINUTE SECOND
do
~~~~~
~~~~~
done < date.txt
この部分で先ほど書き込んだdate.txtを1行ずつ読み込み、処理を行います。
具体的には、date.txtをYEAR MONTH DAY HOUR MINUTE SECONDとして1行ずつ読み込み、doからdoneまでの作図プログラムを実行します。
ifile=$YEAR'-'$MONTH'-'$DAY'_'$HOUR'-'$MINUTE'-'$SECOND'.txt'
で開くファイルを指定します。date.txtを1行ずつ読み込むため、日時の異なるファイルがループで開かれます。
最後に
いかがでしたか。この手法を使えば、大量のファイルでも次々に図化することができるはずです。ぜひお試しください！
Fortranでアメダスデータを処理して平均気温を求める
みなさんこんにちは。おこめまんです。
前回はFortranで大間アメダス10分間値データを読み込むことに挑戦しました。
今回は、Fortranで読み込んだ1か月分の気温データから平均気温を求めてみたいと思います。
今回使用するデータ
今回使用するデータは以前スクレイピングによってゲットした「アメダス大間10分間値データ」です。
アメダス10分間値データをゲットする方法は以下をご覧ください。
＜ファイル名＞
「amd10_31001_201801.csv」
このcsvファイルの中身はこんな感じです。
各列はそれぞれ、
「年」「月」「日」「時」「分」「秒」「降水量」「気温」「平均風速」「風向(平均風速)」「最大瞬間風速」「風向(最大瞬間風速)」「日照時間」
となっています。
よって、最初の1行目には
2018年1月1日0時10分0秒の降水量、気温、平均風速、風向(平均風速)、最大瞬間風速、風向(最大瞬間風速)、日照時間が書かれています。こんな感じ1月31日まで1か月分のデータがあります。
使用した言語
Fortran95
Fortranプログラム
以下を書き込んだプログラム「open_amedas_data.f95」というプログラムを作成します。
program open_amedas_data
implicit none
integer :: year, month, day, hour, minute, second
integer ::counter_temp
real :: prc, temp, ws, wd, gws, gwd, sun
real :: sum_temp, ave_temp
character :: ifile*100
ifile="amd10_31001_201801.csv"
open(30, file=trim(ifile), status='old')
sum_temp = 0
counter_temp = 0
do
read (30, *, end = 99, err=98) year, month, day, hour, minute, second, prc, temp, ws, wd, gws, gwd, sun
if (temp > -900.0) then
sum_temp = sum_temp + temp
counter_temp = counter_temp + 1
end if
!write (6,*) year, month, day, hour, minute, second, prc, temp, ws, wd, gws, gwd, sun
end do
98 write (6, *) "ERROR"
stop
99 write (6, *) "Finish"
close (30)
write(6,*) "counter=", counter_temp
ave_temp = sum_temp/counter_temp
write(6,*) "ave_temp= ", ave_temp
stop
end program open_amedas_data
プログラムの説明
プログラムの説明は以前ご紹介した「Fortranでcsvファイルを開いて1行ずつ読み込む」から追加した内容を扱います。付け足した部分はプログラムをオレンジ色で着色しました。
--------------------------------------------------------
implicit none
integer :: year, month, day, hour, minute, second
integer ::counter_temp
real :: prc, temp, ws, wd, gws, gwd, sun
real :: sum_temp, ave_temp
character :: ifile*100
⇒ 計算処理する際に必要な変数を指定します。counter_tempはデータの個数(10分間値データ1か月分なので、欠測がなければ4464個)を表します。sum_tempは1か月分の気温データを合計したものに相当します。ave_tempは今回求めたい「平均気温」に相当します。
--------------------------------------------------------
sum_temp = 0
counter_temp = 0
⇒ 単純に気温に足せばいいじゃんと思う人もいることでしょう。しかし、はじめのsum_tempとcounterの値を指定していなければうまく計算することができません。
あなたも、「＋４は？」とだけ言われても困ってしまうでしょう。「何」に4を足すのかをはじめに指定しなければなりません。Fortanでは、ループを回す前に最初の値は「0」だと教えてあげましょう。
--------------------------------------------------------
if (temp > -900.0) then
sum_temp = sum_temp + temp
counter_temp = counter_temp + 1
end if
⇒ readで読み込んだ要素のうち、気温に注目して積算していきます。if (temp > -900.0) thenとしているのは、アメダスデータの欠測値を「-999.9」としているので、このときのデータをはじくためです。「もし気温が-900.0℃よりも大きければ、sum_temp = sum_temp + tempとcounter_temp = counter_temp + 1を実行してください」という意味です。正しく観測されていれば、気温は-900.0℃よりは大きい値のはずですよね。counter_tempでは、欠測値を除いたデータ数をカウントしています。欠測がなければ4464になるはずです。
--------------------------------------------------------
write(6,*) "counter=", counter_temp
⇒ さて、doループが終わりました。このとき、sum_tempには1か月分の気温データを積算した値が、counter_tempには欠測値を除いたデータ数が相当していることになります。ここでは、欠測を除いたデータ数を画面に表示させています。
--------------------------------------------------------
ave_temp = sum_temp/counter_temp
write(6,*) "ave_temp= ", ave_temp
⇒ ここで、ついに1か月の平均気温を求めます。平均気温は、積算した気温をデータ数で割れば求まりますよね。計算した平均気温はave_tempに入ります。最後に、write文を書いて、画面に平均気温を表示させましょう。
プログラムを実行して、
counter= 4464
ave_temp= 0.185748115
と画面に出れば正解です。
最後に
いかがでしたか。平均気温は、データ処理のはじめによく扱われます。時間があれば、平均気温のほかに、平均風速も求めてみてください。
それでは、また次回お会いしましょう。おこめまんでした。
Fortranでcsvファイルを開いて1行ずつ読み込む
みなさんこんにちは。おこめまんです。
今回はFortranで処理したいファイルを開いて、読み込む方法をご紹介します。
はじめに
気象屋は多くの場合Fortranを使って気象データの処理を行います。
しかし、Fortran初心者の場合、どうやってファイルの中身を読み込んで処理するんだろうと疑問に思うかと思います。
私もはじめはFortranの仕組みがわからず、ファイルの開き方、データの読み込み方など理解に苦しみました。
そこで、今回は「ダウンロードしたアメダス10分間値のデータをFortranで開いて読み込み、中身を画面に表示する方法」をご紹介したいと思います。
読み込むデータファイルとその中身
今回は、以前スクレイピングでゲットした大間アメダスのデータを書き込んだファイルをFortranで開き、読み込んでみます。
＜ファイル名＞
amd10_31001_201801.csv
＜ファイルの中身＞
はじめの数行を以下にお見せします。
各列はそれぞれ、
「年」「月」「日」「時」「分」「秒」「降水量」「気温」「平均風速」「風向(平均風速)」「最大瞬間風速」「風向(最大瞬間風速)」「日照時間」
となっています。
よって、最初の1行目には
2018年1月1日0時10分0秒の降水量、気温、平均風速、風向(平均風速)、最大瞬間風速、風向(最大瞬間風速)、日照時間が書かれています。
このデータをFortranで読み込んで、画面に表示できるようにしましょう。
使用した言語
Fortran95
Fortranプログラム
以下を書き込んだ「open_amedas_data.f95」というプログラムを作成します。
program open_amedas_data
implicit none
integer :: year, month, day, hour, minute, second
real :: prc, temp, ws, wd, gws, gwd, sun
character :: ifile*100
ifile="amd10_31001_201801.csv"
open(30, file=trim(ifile), status='old')
do
read (30, *, end = 99, err=98) year, month, day, hour, minute, second, prc, temp, ws, wd, gws, gwd, sun
write (6,*) year, month, day, hour, minute, second, prc, temp, ws, wd, gws, gwd, sun
end do
99 write (6, *) "Finish"
98 write (6, *) "ERROR"
close (30)
stop
end program open_amedas_data
プログラムの解説
program open_amedas_data
⇒ Fortranプログラムの名前を書きます。今回は「open_amedas_data.f95」というプログラムを作成したので、「open_amedas_data」となります。
-----------------------------------------------------------------
implicit none
integer :: year, month, day, hour, minute, second
real :: prc, temp, ws, wd, gws, gwd, sun
character :: ifile*100
⇒ A～G列は年、月、日、時、分、秒となっていますので、それぞれyear, month, day, hour, minute, secondとして変数指定を行います。integerとしているのは、A～G列は整数値であるためです。
⇒ H～M列は、気温、平均風速、風向(平均風速)、最大瞬間風速、風向(最大瞬間風速)、日照時間です。これらは実数ですので、realとしてそれぞれ変数指定を行います。
⇒ 開くファイルである「amd10_31001_201801.csv」をifileとしたいためここで変数指定を行います。characterとしているのは、ファイル名を文字列として処理するためです。31001は数字だからintegerなのではと思うかもしれませんが、ファイル名は文字列である必要があります。またifile*100の100という数字ですが、これはifileの文字数が100であるという意味です。今回開くファイルは100文字ではありませんが、多くしていてもこの時点では問題ありません。※open文である処理をする必要があります。
-----------------------------------------------------------------
ifile="amd10_31001_201801.csv"
open(30, file=trim(ifile1), status='old')
⇒ ifile="〇〇"で開くファイルを指定します。
open文でファイルを開きます。30は装置番号で、任意の数字にすることができます。装置番は1～99の間で任意に指定することができますが、5は標準入力(キーボード)、6は
標準出力(画面）となっているため、これらの数字は利用しないようにしましょう。
trimは先ほど、ifile*100として文字列を多くとっていたファイル名の空白をカットしてくれるものです。実際「amd10_31001_201801.csv」は22文字ですから余分な78文字分の空白を切り取ります。trimを書かずそのままfile=ifileとしてしまうとうまく開くことができません。
status='old'は「amd10_31001_201801.csv」はもともと存在するファイルだと宣言するものです。よく利用するのはold、replaceです。詳しい説明は割愛します。
これで、ファイルを開くことができました。
-----------------------------------------------------------------
do
read (30, *, end = 99, err=98) year, month, day, hour, minute, second, prc, temp, ws, wd, gws, gwd, sun
write (6,*) year, month, day, hour, minute, second, prc, temp, ws, wd, gws, gwd, sun
end do
⇒ doループすることで1行ずつファイルの中身を読み込みます。ファイルの読み込みはreadです。装置番号は先ほどopenの際に使用した数字にします。またendで中身を読みこみ、最終行に達した場合に99にジャンプします。errで読み込みの際にエラーがあった場合98にジャンプします。
＜注意＞read文で読み込む際、列の数と変数の数が同じになっているか確かめましょう（year, month, day, hour, minute, second, prc, temp, ws, wd, gws, gwd, sun）。今回開いたcsvファイルは13列であるため、変数も13個です。
writeで先ほど読み込んだ中身を1行ずつ画面に表示させます。装置番号の6は画面表示でしたね。
-----------------------------------------------------------------
99 write (6, *) "Finish"
98 write (6, *) "ERROR"
⇒ 先ほどのend、errのジャンプ先です。endの場合には画面に「Finish」、errの場合には画面に「ERROR」と表示されます。
-----------------------------------------------------------------
close (30)
⇒ 開いていたファイルを閉じます。
-----------------------------------------------------------------
stop
⇒ プログラムをstopさせます。書かなくても今回はOKです。
-----------------------------------------------------------------
end program open_amedas_data
⇒ プログラムの終了を宣言します。
なんとなく理解できたでしょうか。doループのかけ方やデータ型、装置番号など細かい説明は他のサイトを参照してください。
プログラムの実行
さて、先ほどのプログラム「open_amedas_data.f95」を実行してみましょう。
f95 open_amedas_data.f95
./a.out
この2行をコマンドラインで打ってみましょう。
おそらく画面下から上にものすごい勢いで数字が流れているのではないでしょうか。
成功です！
最後に
今回、csvファイルを扱いましたが、txtファイルでも同様に開くことができます（ただし、カンマで区切ったデータである必要があります）。
これで、データ処理の第1歩「ファイルをFortranで開く」ことができました。データ解析の際はファイルをFortranで開き様々な処理を行うことがほとんどです。自由自在にデータを扱えるようにしたいですね。
それではまた次回。おこめまんでした。
Excelで日本語を書いたらLinuxで文字化けする！文字コードを変換しよう
こんにちは、チャリおじさんです。
みなさんはこういう経験ありませんか?
｢viでファイルを開いたら、日本語部分が文字化けしてしまった！｣
もしやあなた、そのファイルはExcelで作ったのでは?
でも大丈夫！文字コードを変換すれば、ちゃんとviで正しく表示されて、編集したり、コマンド操作できるようになりますよ。
なぜ文字化けするのか
Excelというソフトは、【Shift-JIS】という文字コードで保存するようにできています。これは変更できないようです。
対して、viなどLinuxターミナルは、【UTF-8】という文字コードでファイルを読み込むようになっています。
ファイルの文字コードと読み込む側の文字コード設定が一致しないと、文字化けしてしまいます。
viで編集するだけなら
viでは読み込む文字コードを変更することができます。
コマンドviの中で編集する場合は簡単に文字化けを直すことができます。
コマンドモードで、↓を打ってみましょう。
:e ++enc=sjis
どうですか、正しく表示されたと思います。
ターミナル上でコマンド操作するなら
しかし、ターミナルのコマンドライン上で操作する場合はファイルの文字コードを【UTF-8】にしないとエラーになってしまいます。
｢grepコマンドで”気温”という文字列を含むファイルを探す｣とか、
｢sedコマンドで”南”を”北”に置換する｣とか、できないのです。
これを解決するために、nkfコマンドを使って、ファイルの文字モードを変換します。
nkfコマンドをインストールする
おそらく初期状態だと使えないので、nkfをネットからダウンロード後、解答、インストールする必要があります。
やり方は以下のサイトを参考にしてください。
開発ツール/Cygwinにnkfをインストールする方法 - Windowsと暮らす
nkfコマンドの使い方
$ nkf -[オプション] 元ファイル名 > 変更後ファイル名
で、ファイルの文字コードを変更できます。
元ファイルに上書きしたい場合は、
$ nkf -[オプション] --overwrite 元ファイル名
でできます。
オプションには、何の文字コードに変換するかを指定します。
今回は【UTF-8】に変換するので、オプションは -w です。
例えばaaa.csvを変換して上書きする場合は、
です。
ちなみに
以下のオプションも便利そうです。
-g 元ファイルの文字コードを検索する
-s Shift-JISに変換する
-e EUCに変換する
-Lu 改行コードをWindowsからLinuxに変換する (CR+LF→LF)
改行コードについては、また別の記事で紹介しましょう。
最後に
いかがでしたか。
無事に日本語がgrepやsedに引っかかってくれるようになったでしょうか。
それではまた次回。チャリおじさんでした。
バイナリファイルをテキストファイルに変換したい！fortranでやってみよう。
みなさん、こんにちは。チャリおじさんです。
今回は、バイナリファイルをテキストファイルに変換するfortranプログラムをご紹介します。
レーダーエコーとかMSMとかJRAとか、多くの気象データはバイナリファイルで提供されています。
それを統計処理するには、このテキストファイルに変換する操作が必須です。
ここでは例として、｢気象庁が公開している全国合成レーダーGPVのバイナリデータを、fortranでテキストファイルに書き出すプログラム｣をご紹介したいと思います。
バイナリデータについて
レーダーエコー値が記録されている、拡張子が.binになっているファイルのことです。
このバイナリファイル、普通にviとかエディタで開いても、中身を見ることができません。
なので統計解析したいときは、バイナリファイルそのままでは計算できないのです。
よってこのバイナリファイルをfortranを使って、中身が読めるようにテキストファイルに書き出します。
バイナリファイルの構成と変換
バイナリファイルは、010101...のように0と1の羅列がずっと並んでいて、2進法で書かれています。
途中に特に区切りや改行はなく、0と1が永遠に続きます。
これを計算できるようにするためには、10進法に変換する必要があります。
2進法と10進法の関係として、｢2進法の32桁で10進法の1つの値を表す｣というのがあります。(これについて詳しくは他のサイトを参照してください。)
なので、0と1の羅列を32桁ごとに区切って読めば、2進数を10進数に変換することができます。
2進法の1桁は1bitです。なので32桁は32bitです。さらに8bit=1byteです。なので32桁は4byteになります。
よってfortranのプログラムでは、変数定義の際に、｢これは4byteの変数ですよ｣と宣言してあげれば、2進法のデータが10進法に変換されて入力されます。
全国合成レーダーGPVとは
要は雨雲レーダーのことです。
気象庁が提供するデータですが、こちらのサイト(京都大学)で無料で取得することができます。
バイナリファイルの中にはレーダーエコー値(dbZ)が羅列されています。
プログラムの中で、緯度経度を計算して、｢緯度 経度 値｣のように出力できれば、今後の統計解析で役立つでしょう。
プログラム
オレンジ色は入力するデータに合わせて値を入力する箇所です。適宜変更してください。
青色の箇所について、プログラムの下に説明をまとめています。
ーーー
program
implicit none
integer, parameter :: x=(経度グリッド数), y=(緯度グリッド数) !(1)
integer :: i, j
real :: lat, lon
real(kind=4), dimension(x,y) :: figure !(2)
open(unit=11,file='infile.bin',form='unformatted',access='direct',status='old',recl=x*y*4) !(3)
open(unit=31,file='outfile.txt',status='replace') !(4)
read(unit=11,rec=1) figure
do j = 1, y
lat = (南端の緯度) + (緯度間隔)*(j-1) !(5)
do i = 1, x
lon = (西端の経度) + (経度間隔)*(i-1) !(5)
write(31,'(f11.6,f11.6,f10.4)') lat, lon, figure(i,j) !(6)
end do
end do
end program
ーーー
以上です。
ポイント
(1)x,yは定数です。縦横のグリッド数を入力してください。
(2)kindは種別値といって、byte数を指定するものです。
(3)バイナリファイルは｢書式なしファイル｣なので、｢form='unformatted'｣を指定します。
｢access='direct',recl=x*y*4｣によって、レコード長が4byteに固定されます。
(4)出力ファイルは通常通りopenします。別に10進数とか宣言しません。
(5)レーダーエコー値を緯度と経度と並べて出力します。バイナリデータには緯度経度の情報は入っていないので、自分で入力する必要があります。
(6)ほしいレーダーエコー値の桁数に合わせて、編集記述子(f11.6とか)を変更してください。
最後に
いかがでしたか?
これできっと、｢緯度 経度 レーダーエコー値｣という値が羅列したテキストファイルが出力されたでしょう。
それではまた次回、チャリおじさんでした。
