内部波の可視化ツール

境界面変動の時系列図

import glob
import matplotlib.pyplot as plt
import numpy as np
from scipy import interpolate

def get_z_centers(dzs):
  z_pos  = 0.5*dzs[0]
  z_centers = [z_pos]
  for k in range(len(dzs)-1):
    z_pos += 0.5*(dzs[k]+dzs[k+1])
    z_centers.append(z_pos)
  return z_centers

def get_data(file):
  print ("processing... ", file)
  columns = []
  for line in open(file).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]),
                      "bottom_k":int(items[7]),    "surface_k":int(items[8]),
                      "bottom_z":float(items[9]),  "surface_z":float(items[10]),
                      "variables":{}})
    else:
      columns[int(items[0])]["variables"][items[1]]=[float(x) for x in items[2:]]
  return time, dzs, columns

def find_nearest_column(columns, x, y):
  min_dist = np.finfo(np.float64).max
  n = 0
  for num, column in enumerate(columns):
    dist = np.sqrt((x-column["center_x"])**2+(y-column["center_y"])**2)
    if dist < min_dist:
      min_dist = dist
      n        = num
  return n

def get_depth(var, z_centers, value):
  f = interpolate.interp1d(var, z_centers)
  return f(value)


def main():

  files = glob.glob("intwave-0*.dat")
  files.sort()

  # wave gage setting
  x_list=[0.5,1.0,1.5,2.0,2.5,3.0,3.5,4.0,4.5,5.0,5.5,6.0]# m
  data = get_data(files[0])
  target_column_numbers = [find_nearest_column(data[2], target_x,0.05) for target_x in x_list]
  z_centers             = get_z_centers(data[1])

  upper_s  = 10  # salinity
  lower_s  = 20  # salinity

  upper_z   = [[] for i in range(len(x_list))]
  lower_z   = [[] for i in range(len(x_list))]

  t = []
  for file in files:
    data    = get_data(file)
    time    = data[0]
    columns = data[2]
    for n, target_num in enumerate(target_column_numbers):
      var = columns[target_num]["variables"]["salinity"]
      upper_z[n].append(get_depth(var, z_centers, upper_s))
      lower_z[n].append(get_depth(var, z_centers, lower_s))

    t.append(time)


  for n in range(len(x_list)):
    plt.plot(t, upper_z[n], label="salinity = "+str(upper_s))
    plt.plot(t, lower_z[n], label="salinity = "+str(lower_s))
    plt.gca().set_ylim(0,0.3)
    plt.title("Temporal variation of upper and lower interfaces at "+str(x_list[n])+" m")
    plt.ylabel("z (m)")
    plt.xlabel("time (s)")
    plt.legend()
    plt.savefig("time_variation_"+str(x_list[n])+"m.png")
    plt.clf()

if __name__ == '__main__':
  main()

import glob
import matplotlib.pyplot as plt
import numpy as np
from scipy import interpolate

def get_z_centers(dzs):
  z_pos  = 0.5*dzs[0]
  z_centers = [z_pos]
  for k in range(len(dzs)-1):
	z_pos += 0.5*(dzs[k]+dzs[k+1])
	z_centers.append(z_pos)
  return z_centers

def get_xz(columns, dzs):
  x_centers  = []
  for column in columns:
	  x_centers.append(column["center_x"])
  z_centers = get_z_centers(dzs)
  x_matrix = np.array([x_centers]*len(dzs)).transpose()
  z_matrix = np.array([z_centers]*len(x_centers))
  return x_matrix, z_matrix

def get_var(columns, name):
  return np.array([column["variables"][name] for column in columns])

def get_data(file):
  print ("processing... ", file)
  columns = []
  for line in open(file).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]),
					  "bottom_k":int(items[7]),    "surface_k":int(items[8]),
					  "bottom_z":float(items[9]),  "surface_z":float(items[10]),
					  "variables":{}})
	else:
	  columns[int(items[0])]["variables"][items[1]]=[float(x) for x in items[2:]]
  return time, dzs, columns

def get_depth(var, z_centers, value):
  f = interpolate.interp1d(var, z_centers)
  return f(value)

def main():
  files = glob.glob("intwave-0*.dat")
  files.sort()
  data    = get_data(files[0])
  x, z    = get_xz(data[2], data[1])

  for file in files:
	data   = get_data(file)
	s_data = get_var(data[2], "salinity")
	fig    = plt.figure(figsize=(8,4))
	plt.title("Second mode internal wave")
	plt.xlabel("$x$(m)")
	plt.ylabel("$z$(m)")
	plt.contourf(x,z,s_data,np.arange(0,30,1),cmap="jet")
	plt.gca().set_xlim(0,6)
	plt.gca().set_ylim(0,0.3)
	plt.fill_between([3, 6], [0.0, 0.3], color = "0.5") # sloping boundary
	plt.savefig("fig_"+file[-9:-4]+".png")
	plt.clf()

if __name__ == '__main__':
  main()