今回用いる白黒画像は多数の1バイトデータからなります.1バイトのデータ(8個の0と1)をそのまま書き込んだファイルをバイナリファイルと呼びます.2の8乗=256個の整数を表現できます.一方で今まで扱ってきた,人間が読めるファイルは「テキストファイル」といって区別されます.シフトJIS(macのデフォルト文字コード)では,半角文字で1バイト,全角文字で2バイトのメモリを使います.256は3桁ですから,半角の数字で表現した場合3バイト必要になります.したがって,テキストでファイルに書き出した場合とバイナリで書き出した場合とで,必要な記憶容量が3倍も異なることになります.桁数が多くなればもっと差が大きくなります.
今回あつかうデータはバイナリファイルのため,バイナリファイルの読み書きを行う必要があります.次のプログラムは,0.1から1.0までの10個の単精度実数を生成し,それをバイナリ形式でファイル10.datに書き込むプログラムです.
implicit none
integer i
real x
open(1,file='10.dat',access='direct',form='unformatted',recl=4)
do i=1,10
x=0.1*i
write(1,rec=i) x
enddo
end
ここで,form='unformatted'でバイナリ形式を指定し,recl=4で一つのデータ(レコード)の大きさ(今の場合単精度実数=4バイト=32ビット)を指定します.access='direct'は,直接参照ファイルであることを指定します.ファイル中の場所を指定することで,特定の数値に直接アクセスします.今までは'sequential'となっていて,ファイルの先頭から順番にアクセスすることが暗に指定されていました.またwrite文では,rec=iでレコード番号を指定します.このプログラムを実行し,生成された10.datのサイズが40バイトになっていることを確認しましょう.
先ほど生成したデータを読み込んで画面に表示するプログラムは以下になります.
implicit none
integer i
real x
open(1,file='10.dat',access='direct',form='unformatted',recl=4)
do i=1,10
read(1,rec=i) x
write(*,*) x
enddo
end
実行して結果を確認しましょう.
今回扱う画像は白黒256階調なので,各ピクセルは1バイトのデータで表現されます.FORTRANでは符号なし1バイトデータを扱う適当な形式が用意されていないので,character文で宣言される文字型変数に数値を格納します.
それぞれの文字には文字コードと呼ばれる数字が対応しています.そこで,0から255までの整数を文字に一度変換し,その文字をファイルに書き込むと,0から255までの整数(本来は文字コード)がファイルに書き出されます.
0から255までの整数をバイナリデータとしてファイル256.datに書き込み,そのデータを読み出して0から255までの整数を画面に表示するプログラムは以下になります.
implicit none
character n
integer i
open(1,file='256.dat',access='direct',form='unformatted',recl=1)
do i=0,255
n=char(i)
write(1,rec=i+1) n
enddo
do i=1,256
read(1,rec=i) n
write(*,*) i,ichar(n)
enddo
end
ポイントは,256階調の整数はchar()で文字型変数 n に変換してファイルに書き込む.読み出す際は文字型変数 n にまずは代入し,それをichar()関数で256の整数に変換する です.
このデータは,1980年-2000年における全地球の表面気温のデータです.360x180=64800個の1バイトデータ(0-255)が並んだファイルが21個あります.0が-15度,255が25度に対応します.それぞれのデータは地球上のある経度(0-359),緯度(0-179)における温度です.注意してほしいのは,二次元的にデータが並んでいるわけではなく,一次元的にデータは並んでいることです.データの頭は北極(緯度0)であり,まずは経度方向に360個のデータが並びます.361番目のデータは緯度1,経度0となります.360個のデータの組みが180個並んでいる,という構造になってます.
ダウンロードしたら
mv ~/Downloads/data/* ./とすることで現在のディレクトリの中にデータファイルが移動します.lsで確認してください.
display -depth 8 -size 360x180 gray:conv1980.raw を実行すると,1980年における温度分布が画像として表示されます.黒が低温,白が高温です.displayコマンドが使えない場合はこちら
次に,以下のプログラムを実行し,数値(0-255)が画面に表示されることを確認してください.
implicit none
integer i
character n
open(1,file='conv1980.raw',access='direct',form='unformatted',recl=1)
do i=1,64800
read(1,rec=i) n
write(*,*) i,ichar(n)
enddo
end
今回は,1980年から2000年までの21年分のデータを解析し,その温度上昇のトレンドが地域ごとにどう分布しているかを調べます.
そのためには,ファイルを21個開かないといけません.正直にプログラムに記述してもいいのですが効率が悪いので,文字処理のテクニックを用いてファイルを操作します.次のプログラムを実行すると,21個のファイルを順番にオープンしてそのファイル名が画面に表示されます.
implicit none
integer i,num,nyear,j
character name*12,year*4,n
nyear=21
do i=1,nyear
num=1979+i
write(year,'(I4)') num ! 整数numを文字型に変換しyearに代入
name='conv'//year//'.raw' ! conv, year, .rawからファイル名を生成
open(10+i,file=name,access='direct',form='unformatted',recl=1)
! ここでnameには上で作られた文字列が代入される
do j = 1, 64800
read(10+i,rec=j) n
write(*,*) num,j,ichar(n)
enddo
enddo
end
write(year,'(I4)') num で,4桁の整数numが文字となって文字型変数yearに格納されます.
さてここまでくると,全地球の温度変化を解析できるようになりました.以下の手順を踏んで,地球上のどこで温暖化がより進行しているのか解析しましょう.
implicit none
! 宣言文など
open(1,file='plus.raw',access='direct',form='unformatted',recl=1) ! 書き込むファイルのオープン
do i=1,nyear
num=1979+i
write(year,'(I4)') num
name='conv'//year//'.raw'
write(*,*) name
open(10+i,file=name,access='direct',form='unformatted',recl=1) !各年のファイルのオープン
do j = 1, 64800
! 10+i番のファイルから数値(0-255)を64800個読み込み,温度に変換し配列T(i,j)に格納する.読み出し方は
read(10+i,rec=j) n
T(i,j) = ichar(n) ! 配列T(i,j)に温度データ(-30度-30度)が代入される
enddo
enddo
do j = 1, 64800
do i = 1, nyear
! この2重ループで,最小二乗法により64800個の傾きを算出する
! xが年(1980-2000) yが温度T(年i,場所の番号j)となっている.
enddo
enddo
do j = 1, 64800 ! 64800個の傾きの書き出し
! ここで温度上昇の傾きから0-255の整数に変換する
! プラスの場合とマイナスの場合と二通り実行する必要あり
n = char(変換後の整数がここに入る)
write(1,rec=i) n ! ファイルへの書き込み
enddo
画像がだせたら,以下の点を観察してください.