水面のトレーサー分布を表示するスクリプト

import numpy as np
import glob
import matplotlib.pyplot as plt
from numpy import linspace, meshgrid
from scipy.interpolate import griddata
 
def grid(x, y, z, resX=100, resY=100):
    xi = linspace(min(x), max(x), resX)
    yi = linspace(min(y), max(y), resY)
    X, Y = meshgrid(xi, yi)
    Z = griddata((x, y), z, (X, Y), method='linear') 
    return X, Y, Z
 
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":
      key_names =items[1:]
    elif items[0] == "column":
      dic = {"variables":{}}
      for n, name in enumerate(key_names):
        dic[name] = items[n+1]
      columns.append(dic)
    else:
      columns[int(items[0])]["variables"][items[1]]=[float(x) for x in items[2:]]
  return time, dzs, columns
 
##### specific setting
file_names = glob.glob("./abashiri*.dat")
target_tracer = "temperature"
reso_x = 200
reso_y = 200
##### 

fig  = plt.figure(figsize=(10,10))

for num,fname in enumerate(file_names):
  data = []
  x    = []
  y    = []
 
  time, dzs, columns = get_data(fname)
 
  for column in columns:
    k_level = int(column["sfc_k"])
    if column["dry_wet"]=="wet":
      data.append(float(column["variables"][target_tracer][k_level]))
    else:
      data.append(np.nan)
    x.append(float(column["center_x"]))
    y.append(float(column["center_y"]))
 
  X, Y, Z = grid(x, y, data, reso_x, reso_y)
  cs = plt.contourf(X, Y, Z)
  plt.colorbar(cs)
  plt.show()