2010-7-23:全観測点の年平均降水量とその標準偏差を出力する - BioEoS::環境解析学研究室::

全観測点の年平均降水量とその標準偏差を出力する

 #(1)観測点番号　平均値　標準偏差を求める(MonthlyAnnual/)
 #(2)観測点番号　経度　緯度を引き出す(STATION.DAT)
 #(3)観測点番号　経度　緯度　平均値　標準偏差を含むファイルを作る
 ################################################################################
 =begin
 まず基本的な考え方ですが、
 「観測点番号」をkeyにして、2つ(またはそれ以上)のデータ
 （今の場合は、STATION.DATと降水量データ)を結び付ける。
 これをするには、(2)より観測点番号をkeyにして、経度・緯度・高度・観測地点名
 をvalueとするような配列を作成し、そこにkeyとvalueを記憶させる。
 こうした配列は、普通の配列(Array)を使ってもできるが、
 連想配列(Hash)を使うと、もっと便利で計算速度も速い。
 こうして作ったmethodが、下の def station_datです。
 ここで、Hash(連想配列)名は、stationです。
 なお、緯度・経度の計算が正しくできるように、以前のものを変更し
 もう一つ別のddmmというmethodを作りました。
 =end
 ################################################################################
 #(2)観測点番号　経度　緯度を引き出す(STATION.DAT)
 def station_dat
   station = {}                         # station (Hash) の初期化
   station_file = "STATIONS_new.dat"        # STATIONS.DATを  station_new.dat
   open(station_file).each do |line|    # オープンし、各行読み取り
     array = line.chop.split
     if array.size == 8
       station_number = array[0]        # 観測点番号
       longitude = ddmm(array[2])       # 経度
       latitude = ddmm(array[1])        # 緯度
       altitude = array[3]              # 高度
       name = array[7]                  # 観測点名
       # station_numberをキーにして、観測点名、経度、緯度、高度の値を
       # 配列にしたものを、連想配列stationに格納
       station[station_number] = [name,longitude,latitude,altitude]
     end
   end
   station                              # 連想配列stationを戻す
 end
 # ddmmの表現の緯度経度の値をdegにする
 def ddmm(dm)
   x = dm.to_i                          # 元の値を整数化
   y = x / 100                          # 元の値を100で割る→度の値
   z = (x % 100).to_f/60.0   # 元の値を100で割ったときの余り(分に相当)を60で割る
   y + z                     # どの値(整数)に小数点以下を加えて戻りにする
 end
 ################################################################################
 class Array                #Arrayという配列を作る
   # 配列の全要素の平均を算出する
   def average　　　　　　　　　#averageは変数、初期化　　　
     inject(0.0){|r,i| r += i}/size.to_f　　　　　　#sizeは要素の数
   end
   # 配列の全要素の分散を算出する
   def variance
     a = average
     inject(0.0){|r,i| r += (i - a)**2 }/size.to_f
   end
   # 配列の全要素の標準偏差を算出する
   def standard_deviation
     Math.sqrt(variance)
   end
 end
 ################################################################################
 #(1)観測点番号　平均値　標準偏差を求める(MonthlyAnnual/)
 # ある一つのファイルに対して、47年分の年降水量を読み取り、
 # 配列に格納した後、平均年降水量とその標準偏差を返す
  def mean_stdv(dir_file)
   n = 0                                # ヘッダを読み飛ばすための準備
   array = []                           # 各年の年降水量を格納する配列の初期化
   open(dir_file).each do |line|        # dir_fileをオープンし、各行を読み取る
     if n > 0                           # 2行目以降(ヘッダでない)なら
       v = line.chop.split              # lineから空白で分割して配列vを作成
       array << v[13].to_f              # 年降水量を配列arrayに格納
     end
     n += 1                             # 各行の読み取りと作業が済めば、行数加算
   end
 # 平均年降水量とその標準偏差を返す
   [array.average,array.standard_deviation]　　#arrayは配列名、averageは上に定義したmethod
 end
 #############################################################################
 #(2)観測点番号　経度　緯度を引き出す(STATION.DAT)
 station = station_dat                  # 戻り値station: 連想配列
 ##########################################################################
 #(3)観測点番号　経度　緯度　平均値　標準偏差を含むファイルを作る
 # 全観測点の年平均降水量とその標準偏差を出力するファイル
 #out_file =  "c:/GIS/china/buchashuju" + "/" + "mean_annual_prec" + ".yr"
 out_file = "new_mean_annual_prec"
 out = open(out_file,"w")
 # ヘッダの出力
 #out.puts "st_num  name           longitude  latitude   altitude   mean   stdv" 
 dir = "MonthlyAnnual/"                 # 元ファイルのあるディレクトリ
 Dir.foreach(dir) do |file|
   if /mon$/ =~file
   # この観測点の47年分の年降水量から、平均年降水量とその標準偏差を計算
     mean,stdv = mean_stdv(dir + file)
     station_number = file[0..-5]       # 観測点番号
   # この観測点の結果の出力
     if station[station_number]         # 降水量観測点が全観測点記録中にあれば
     # 連想配列stationから、station_numberの経度緯度等の値を取り出す
       name,longitude,latitude,altitude = station[station_number]
       out.printf "%5s|%11.5f|%10.5f|%7.1f|%7.1f|%7.1f|%  s\n",station_number,longitude,latitude,altitude,mean,stdv,name
     end
   end
 end
 out.close                              # 出力ファイルのクローズ

結果は以下の通りです

QGISで改めでstationのベクタレイヤを作ります。

cat new_mean_annual_prec | v.in.ascii cat=1 out=new_mean_annual_prec x=2 y=3 col="cat integer, longitude double, latitude double, altitude varchar(10), mean double, stdv double, name varchar(30)"

今回はここまでになります。
お疲れ様でした。