12/13/2019

1. time: 計算結果の時間を表しています．


2. dz: 計算セル鉛直方向の厚さを，下（底）から上（水面）方向に列記しています．
   
   
3. limit_depth: 冠水・干出の限界水深を表しています．この水深よりも浅いときは陸地(dry)になります．
   
   
4. items: fantomの出力データは鉛直セルコラム構造が基本となっています．itemsはこのcolumnの基本項目を示しています．
   
   
   • number : 全コラムに付けた通し番号です（構造格子ではないので，番号順に注意が必要です．）．
   • k_offset: 現在使っておりませんので，無視して問題ありません．
   • center_x, center_y: コラム中心のx, y座標です（この情報を利用して場所を特定することが多いです）．
   • delta_x, delta_y: コラム（セルと同じ）の水平方向サイズです．
   • btm_k: 水底が存在する鉛直インデックスです（最下層セルが0になります）．
   • sfc_k: 水面が存在する鉛直インデックスです（最上層セルはdzの数-1になります）．
   • btm_z: コラム内の水底の鉛直座標です．
   • sfc_z: コラム内水面の鉛直座標です．
   • container_i, container_j: コラムが属しているコンテナのデカルト座標インデックスです．
   • local_i, local_j: コンテナ内におけるコラムのデカルト座標インデックスです．
   • dry_wet: コラムの冠水・干出を示します(dry, wetのどちらか)．
   
   
5. columnで始まる行は，各コラムのitemsで示した情報を記載しています．この行はコラムの数だけあります．
   
   
6. 数字から始まる行は，数字に対応する通し番号のコラムのデータを示しています．数字の後に続く名前が変数名で，その後の数値の列は，最下層のセルから最上層のセルまでのデータです（水底より下，水面より上のセルも含みます）．
   
   最後に，コラム中心の水平座標から対象コラムを選び出し，そのコラム（一点）における水面の時間変動をプロットするPythonスクリプトの例を示します．
   
   

   	import numpy as np
   	import matplotlib.pyplot as plt
   	import glob
   	# extract data from a file
   	def get_data(file_name):
   	columns = []
   	for line in open(file_name).readlines():
   	items = line.rstrip().split()
   	if items[0] == "time":
   	time = float(items[1])
   	elif items[0] == "dz":
   	dzs = [float(x) for x in items[1:]]
   	elif items[0] == "limit_depth":
   	pass
   	elif items[0] == "items":
   	pass
   	elif items[0] == "column":
   	columns.append({"number" :int(items[1]), "k_offset" :int(items[2]),
   	"center_x" :float(items[3]), "center_y" :float(items[4]),
   	"delta_x" : float(items[5]), "delta_y" :float(items[6]),
   	"btm_k" :int(items[7]), "sfc_k" :int(items[8]),
   	"btm_z" :float(items[9]), "sfc_z" :float(items[10]),
   	"container_i":int(items[11]), "container_j":int(items[12]),
   	"container_i":int(items[13]), "container_j":int(items[14]),
   	"status" :items[15], "variables" :{}})
   	else:
   	columns[int(items[0])]["variables"][items[1]]=[float(x) for x in items[2:]]
   	return time, dzs, columns
   	if __name__ == '__main__':
   	# gathering dat files from folder
   	files = sorted(glob.glob("kasumi*.dat"))
   	# specifying the location for plotting surface level
   	target_x = 11375
   	target_y = 10625
   	# temporary read first file to obtain columns structure
   	t, dzs, columns = get_data(files[0])
   	# finding the target column from location
   	for column in columns:
   	if column["center_x"]==target_x and column["center_y"]==target_y:
   	target_number = column["number"]
   	break;
   	# collecting surface level from all the files
   	time = []
   	surface_level = []
   	for file in files[:]:
   	print ("processing... ", file)
   	t, dzs, columns = get_data(file)
   	time.append(t/3600.0) # sec -> hour
   	surface_level.append(columns[target_number]["sfc_z"])
   	# plot time variation of surface level
   	plt.title("Time variation of surface level")
   	plt.plot(time, surface_level, 'k-x')
   	plt.xlabel("Time (hours)")
   	plt.ylabel("Level (m)")
   	plt.savefig("surface_variation.png")
   	plt.show()

   view raw surface_displacement.py hosted with ❤ by GitHub