#!/usr/bin/env ruby

require "getopts"
require "numru/netcdf"
module NumRu
  class NetCDFVar
    def scaled_get(hash=nil)
      sf = att('scale_factor')
      ao = att('add_offset')
      if ( sf == nil && ao == nil ) then
	# no scaling --> just call get
	simple_get(hash)
      else
	if (sf != nil) 
	  csf = sf.get
	       if csf.is_a?(NArray) then  # --> should be a numeric
		 csf = csf[0]
	       elsif csf.is_a?(String)
		 raise NetcdfError, "scale_factor is not a numeric"
	       end
	  if(csf == 0) then; raise NetcdfError, "zero scale_factor"; end
	else
	  csf = 1.0      # assume 1 if not defined
	end
	if (ao != nil) 
	  cao = ao.get
	  if cao.is_a?(NArray) then  # --> should be a numeric
	    cao = cao[0]
	  elsif csf.is_a?(String)
	    raise NetcdfError, "add_offset is not a numeric"
	  end
	    else
	  cao = 0.0      # assume 0 if not defined
	end
	var = simple_get(hash)
	var = var.to_type( NArray::FLOAT )
	csf*var+cao
      end
    end  
    alias put scaled_put
    alias get scaled_get
  end
end

# $LOAD_PATH に /work31/daktu32/ecmwf/bin を追加. これで draw.rb を require 可能.
# ローカルの PC で使う場合は, 任意のディレクトリを指定すること.


$: << '/GFD_Dennou_Club/dc-arch/dcchart/daktu32/ERA40/bin'
$: << '/work31/daktu32/ERA40/bin'
$: << '~/work/ecmwf/eva/bin'

require "libgphys-e" # 内部で numru/ggraph を呼び出している
require "libdraw"
include Draw


# オプション解析 ----------------------------------------------------
unless getopts("hHfnc", "help", "file", "native", "mean:", "itr:", "title:", "range:", "nocont", "line", "map", "exchange")
  print "#{$0}:illegal options. please exec this program with -h/--help. \n"
  exit 1
end

file = Array.new

# 引数中の nc ファイルを
ARGV.each do |i| 
  if i =~ /.nc$/ then
    p i
    file << i
  end
end


# 変数 : デフォルトはヘッダ の一番最後の変数を自動的に選択
if ARGV[ARGV.size-1] =~ /.nc$/ then
  var = NetCDF.open(file[0], "r").var_names[-1]
else # 陽に指定も可能. その場合, 一番最後に与える.
  var = ARGV[ARGV.size-1]
end


if ($OPT_h) || ($OPT_H) || ($OPT_help) || ARGV.size == 0then
  print <<EOF

====================================================
           #{File.basename($0)} ver.0.1
====================================================

Describe:

  ECMWF の 再解析データ netCDF ファイルから絵を描くスクリプト.
  指定するデータは, 帯状平均, 時間, 月平均がかかったものを前提とする.

  (ex. mean-V_2001_01.nc, mean_T_2001_01.nc ... 頭に mean がくっついているファイル) 

need:

  同ディレクトリの draw.rb から関数を呼び出している. draw.rb は Draw モジュールを定義している.

USAGE: #{File.basename($0)} <opt> [ncfile [ncfile ..] ] [var]
  
  * 引数解説

    [ncfile] : 対象とする netCDF ファイル名を指定. (ex. T_2001_01.nc)
               複数指定可能. その場合, 変数値を各格子毎に平均したものを出力する.
    [var]    : 対象とする変数名を指定. デフォルトの値は 適当な変数を出力.

  * オプション

    -h, -H, --help : このメッセージを表示.
    -f, --file     : 出力結果は dcl.ps へ. (デフォルトは標準出力.) 

HISTORY:

  2003/12/10  created  by daktu32@ep.sci.hokudai.ac.jp

====================================================

EOF
exit 1
end

# オプション解析 終了----------------------------------------------------

p var

if file.size == 1 then
  gphys = GPhys::NetCDF_IO.open(file[0].to_s, var)
 p  $title = File::basename(file[0].to_s, ".nc")

else 
  if ($OPT_c)
    gphys = GPhys::NetCDF_IO.open(file, var)
  else
    gphys = mean_gphys(file, var)
  end
end


# range
if ($OPT_range)
  range = ($OPT_range).to_s.split(/\s*,\s*/)
  range.each{|rg|
    ran = rg.to_s.split(/\s*=\s*/)
    if ran[1] =~ /\s*\.\.\s*/
      rg1=ran[1].split(/\s*\.\.\s*/)[0].to_f
      rg2=ran[1].split(/\s*\.\.\s*/)[1].to_f
      gphys = gphys.cut(ran[0] => rg1..rg2)
    else
      gphys = gphys.cut(ran[0] => ran[1].to_f)
    end
  }
end


# title
$title = ($OPT_title) if ($OPT_title)

# mean
if ($OPT_mean)
  mean = ($OPT_mean).to_s.split(/\s*,\s*/)
  mean.each{|mn|
    gphys = gphys.mean(gphys.axnames.index(mn))
  }
end


if  var == "t" then
  Draw::settei_T
elsif var == "u" || var == "uwnd" then 
  Draw::settei_U(gphys)
#  gphys = gphys.cut('latitude'=>0, 'time'=>1200, 'levelist'=>1000)
elsif var == "v" then 
  Draw::settei_V(gphys)
elsif var == "w" then 
  Draw::settei_W(gphys)
elsif var == "cp" then
  Draw::settei_CP
  Draw::mkwin
  Draw::draw(gphys, $levels, $patterns, true, true, false)
  exit
elsif var == "tsrc" then
  Draw::settei_tsrc
  Draw::mkwin
elsif var == "lsp" then  
  Draw::settei_LSP
  Draw::mkwin
  Draw::draw(gphys, $levels, $patterns, true, true, false)
  exit
elsif var == "tp" then  
  Draw::settei_TP
#  Draw::settei_CP
  Draw::mkwin
  if (gphys.axis(0).pos.name == "longitude") && (gphys.axis(1).pos.name == "latitude")
    Draw::draw(gphys, $levels, $patterns, true, true, false)
  else
    Draw::draw(gphys, $levels, $patterns, true, false, false) if !($OPT_line)
    Draw::line(gphys) if ($OPT_line)
  end
  exit
elsif ($OPT_native) || ($OPT_n) then
  Draw::settei_native
elsif var == "psi" then
  Draw::settei_strm(gphys)
  Draw::mkwin
  Draw::draw(gphys, $levels, $patterns, true, false, true)
  exit
elsif var == "theta" then
  Draw::settei_Theta
  Draw::mkwin
  Draw::draw(gphys, $levels, $patterns, false, false, true)
  exit
elsif var == "pbar" then
  Draw::settei_Available
elsif var == "pws" then
  Draw::settei_PWS
else
  Draw::settei_native
end



Draw::mkwin
#Draw::line(gphys)
Draw::draw(gphys, $levels, $patterns, true, ($OPT_map), !($OPT_nocont), true, ($OPT_line))

#############################################################
 #    一時的に描きたくなった絵用の描画メソッド呼び出し部   #
#############################################################

#p gphys = gphys.cut("time"=>1200, "longitude"=>200, "levelist"=>1000)
#p gphys.data.val.max
#exit
#  Draw::draw(gphys)
#Draw::line(gphys.cut('levelist'=>10..30).cut("latitude"=>0))
#Draw::draw_with_all(gphys)
#, "levelist"=>1000)
#Draw::draw(gphys, $levels, $patterns, true)
#Draw::draw(gphys.cut("levelist"=>10..1000), $levels, $patterns)
#Draw::draw_cont(gphys)