今回用いる白黒画像は多数の1バイトデータからなります．1バイトのデータ（8個の0と1）をそのまま書き込んだファイルをバイナリファイルと呼びます．２の８乗＝２５６個の整数を表現できます．一方で今まで扱ってきた，人間が読めるファイルは「テキストファイル」といって区別されます．シフトJIS（macのデフォルト文字コード）では，半角文字で１バイト，全角文字で２バイトのメモリを使います．２５６は３桁ですから，半角の数字で表現した場合３バイト必要になります．したがって，テキストでファイルに書き出した場合とバイナリで書き出した場合とで，必要な記憶容量が３倍も異なることになります．桁数が多くなればもっと差が大きくなります．
今回あつかうデータはバイナリファイルのため，バイナリファイルの読み書きを行う必要があります．次のプログラムは，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バイトデータを扱う形式が用意されていません．integer(1)は１バイト整数なのですが，符号がついて-128-+127の範囲となります．
0から255までの整数をバイナリデータとしてファイルに書き出し，そのデータを読み出して0から255までの整数を画面に表示するプログラムは以下になります．ファイルに書き出す際は，127よりも大きな数は対応する補数（マイナスの数）に変換します．参考までに．読み出す際には変換前の変数j(-128から+127)も表示されます．128からは，以前に教えた規則通りに二進数から１０進数に変換され，jはマイナスとなります．
implicit none
integer :: i, k
integer(1) :: j
open(1,file='256.dat',access='direct',form='unformatted',recl=1)
do i=0,255 ! 画像データは0-255なので0からループを開始する
if( i > 127 ) then
j = i - 256 ! 128以上の数については補数に変換する
else
j = i
endif
write(1,rec=i+1) j ! rec番号は1から始まるのでiに1を足す
enddo
do i=1,256
read(1,rec=i) j ! 0から127, -128から-1の数がバイナリで格納されている
if( j < 0 ) then
k = j + 256
else
k = j
endif
write(*,*) i, j, k
enddo
end
このデータは，1980年-2000年における全地球の表面気温のデータです．360x180=64800個の1バイトデータ（0-255）が並んだファイルが21個あります．0が-15度，255が25度に対応します．それぞれのデータは地球上の経度(-179-180)，緯度(0-179)における温度です．注意してほしいのは，二次元的にデータが並んでいるわけではなく，一次元的にデータは並んでいることです．データの頭は北極（緯度0）であり，まずは経度方向に360個のデータが並びます．361番目のデータは緯度１，経度-179度（西経179度）となります．360個のデータの組みが180個並んでいる，という構造になってます．
ダウンロードしたら
mv ~/Downloads/data-2/* ./とすることで現在のディレクトリの中にデータファイルが移動します．lsで確認してください．
display -depth 8 -size 360x180 gray:conv1980.raw を実行すると，1980年における温度分布が画像として表示されます．黒が低温，白が高温です．
次に，以下のプログラムを実行し，数値（0-255)が画面に表示されることを確認してください．
implicit none
integer :: i,k
integer(1) :: j
open(1,file='conv1980.raw',access='direct',form='unformatted',recl=1)
do i=1,64800
read(1,rec=i) j
if( j < 0 ) then
k = j + 256
else
k = j
endif
write(*,*) i,k
enddo
end
今回は，1980年から2000年までの21年分のデータを解析し，その温度上昇のトレンドが地域ごとにどう分布しているかを調べます．
そのためには，ファイルを21個開かないといけません．正直にプログラムに記述してもいいのですが効率が悪いので，文字処理のテクニックを用いてファイルを操作します．次のプログラムを実行すると，21個のファイルを順番にオープンしてそのファイル名が画面に表示されます．
implicit none
integer :: i,num,nyear,j, k
character :: name*12,year*4 ! 文字を格納する変数を宣言 12,4は文字数
integer(1) :: m
nyear=21
do i=1,nyear
num=1979+i
write(year,'(I4)') num ! 整数numを文字型に変換しyearに代入
name='conv'//year//'.raw' ! conv, year, .rawからファイル名nameを生成
open(10+i,file=name,access='direct',form='unformatted',recl=1)
! ここでnameには上で作られた文字列が代入される
do j = 1, 64800
read(10+i,rec=j) m
if( m < 0 ) then
k = m + 256
else
k = m
endif
write(*,*) num,j,k
enddo
enddo
end
さてここまでくると，全地球の温度変化を解析できるようになりました．以下の手順を踏んで，地球上のどこで温暖化がより進行しているのか解析しましょう．
implicit none
integer(1) :: m ! バイナリデータを読み込むための変数
integer :: num, nyear ! ファイルを読み込む際の変数
integer :: i, k, j
character :: name*12,year*4
real :: T(21,64800) ! 温度を格納する配列 21年 x 64800点
real :: xsum, ysum, x2sum, xysum ! 最小自乗法で登場する4種の和
real :: A ! 傾きを格納する変数
open(1,file='plus.raw',access='direct',form='unformatted',recl=1) ! 書き込むファイルのオープン
nyear = 21
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個読み込む．読み出し方は
read(10+i,rec=j) m
if( m < 0 ) then
k = m + 256
else
k = m
endif
T(i,j) = k*??? + ??? ! 配列T(i,j)に温度データ（-15度から+25度）が代入される．???は適当な数値で書き換える（kは０から２５５）
enddo
enddo
do j = 1, 64800 ! 場所の番号
do i = 1, nyear ! 年
! この2重ループで，最小二乗法により64800個の傾きを算出する
! 最小二乗法における x が年（1980-2000) y が温度Tとなっている．
enddo
! ここで温度上昇の傾き A から0-255の整数 k に変換する
if( k > 127 ) then
m = k - 256 ! 128以上の整数は補数に変換
else
m = k
endif
write(1,rec=j) m ! 画像ファイルとして書き出し
enddo
画像がだせたら，以下の点を観察してください．