init commit from UWLCM plotting

Energy_spectrum/spectrum.py

@@ -0,0 +1,75 @@
+# (C) Maciej Waruszewski
+
+import h5py
+import numpy as np
+from sys import argv
+import matplotlib.pyplot as plt
+
+#velocities = ["u", "v", "w"]
+#velocities = ["u", "v", "w", "cloud_rw_mom3", "rv", "th", "RH", "aerosol_rw_mom3"]
+#velocities = ["cloud_rw_mom3"]
+velocities = ["w"]
+
+time_start = int(argv[1])
+time_end = int(argv[2])
+outfreq = int(argv[3])
+from_lvl = int(argv[4])
+to_lvl = int(argv[5])
+
+directories = argv[6:len(argv):2]
+labels = argv[7:len(argv):2]
+print directories, labels
+
+# read in nx, ny, nz
+for directory, lab in zip(directories, labels):
+  w3d = h5py.File(directory + "/timestep" + str(time_start).zfill(10) + ".h5", "r")["u"][:,:,:]
+  nx, ny, nz = w3d.shape
+  Exy_avg = {}
+  for vel in velocities:
+    Exy_avg[vel] = np.zeros(((nx+1)/2))
+
+  for t in range(time_start, time_end+1, outfreq):
+    filename = directory + "/timestep" + str(t).zfill(10) + ".h5"
+    print filename
+
+    for vel in velocities:
+
+      w3d = h5py.File(filename, "r")[vel][:,:,:] # * 4. / 3. * 3.1416 * 1e3 
+
+      for lvl in range(from_lvl, to_lvl+1):
+        w2d = w3d[:, :, lvl]
+
+        wkx = 1.0 / np.sqrt(nx - 1) * np.fft.rfft(w2d, axis = 0)
+        wky = 1.0 / np.sqrt(ny - 1) * np.fft.rfft(w2d, axis = 1)
+
+        #Ex = 0.5 * (np.abs(wkx) ** 2)
+        Ex = (np.abs(wkx) ** 2)
+        Ex = np.mean(Ex, axis = 1)
+        #Ey = 0.5 * (np.abs(wky) ** 2)
+        Ey = (np.abs(wky) ** 2)
+        Ey = np.mean(Ey, axis = 0)
+
+        Exy = 0.5 * (Ex + Ey)
+        Exy_avg[vel] += Exy
+
+      K = np.fft.rfftfreq(nx - 1)
+  #    plt.loglog(K, Exy)
+      lmbd = 50. / K # assume dx=50m
+
+    if (t == time_start and lab==labels[0]):
+      plt.loglog(lmbd, 2e-1* K**(-5./3.) )
+
+  for vel in velocities:
+    Exy_avg[vel] /= (time_end - time_start) / outfreq + 1
+    Exy_avg[vel] /= to_lvl+1 - from_lvl
+  #crudely scale
+  #  Exy_avg[vel] /= Exy_avg[vel][len(Exy_avg[vel])-1]
+    plt.loglog(lmbd, Exy_avg[vel] , linewidth=2, label=lab+"_"+vel)
+
+plt.xlim(10**4,10**2)
+plt.xlabel("l[m]")
+plt.ylabel("PSD")
+plt.legend()
+plt.grid(True, which='both', linestyle='--')
+plt.title("Mean PSD of w 322m<z<642m @3h")
+plt.show()

Energy_spectrum/spectrum_shells.py

@@ -0,0 +1,46 @@
+# (C) Maciej Waruszewski
+
+import h5py
+import numpy as np
+from sys import argv
+import matplotlib.pyplot as plt
+
+N = 65
+
+h = h5py.File(argv[1], "r")
+
+NN = N / 2 + 1
+
+u = h["u"][0:N,0:N,0:N]
+v = h["v"][0:N,0:N,0:N]
+w = h["w"][0:N,0:N,0:N]
+
+
+o = N * np.fft.fftfreq(N)
+
+dk = 2.0
+dkh = dk / 2
+
+K = N * np.fft.rfftfreq(N)
+print K
+E = np.zeros(len(K))
+
+
+uk = np.fft.fftn(u) / N ** 3
+vk = np.fft.fftn(v) / N ** 3
+wk = np.fft.fftn(w) / N ** 3
+
+for i in xrange(N):
+    for j in xrange(N):
+        for k in xrange(N):
+            q = np.sqrt(o[i] ** 2 + o[j] ** 2 + o[k] ** 2)
+            for ix in xrange(len(K)):
+                if abs(K[ix] - q) < 0.5:
+                    e = abs(uk[i,j,k]) ** 2 + abs(vk[i,j,k]) ** 2 + abs(wk[i,j,k]) ** 2
+                    E[ix] += 0.5 * e
+
+K = K[2:]
+E = E[2:]
+
+plt.loglog(K, E, 'k-', lw = 2)
+plt.savefig('spectrum_shells.svg')

Matplotlib_common/latex_labels.py

@@ -0,0 +1,45 @@
+var_labels = {
+  "lwp" : 'LWP [g m$^{-2}$]',
+  "rwp" : 'RWP [g m$^{-2}$]',
+  "surf_precip" : 'Surface precip. [mm/day]',
+  "acc_precip" : 'Accumulated precip. [mm]',
+  "cl_nc" : '$N_c$ [cm$^{-3}$] (cloudy cells)',
+  "cl_nr" : '$N_{r>25\mu m}$ [cm$^{-3}$]',
+  "cl_gccn_conc" : '$N_{GCCN}$ [cm$^{-3}$] (cloudy cells)',
+  "thl" : r'$\theta_l$ [K]',
+  "00rtot" : '$q_{t}$ [g/kg]',
+  "rliq" : '$q_{l}$ [g/kg]',
+  "clfrac" : 'Cloud fraction',
+  "prflux" : 'Precip. flux [W m$^{-2}$]',
+  "wvar" : r'Var$\left(w\right)$ [m$^2$ s$^{-2}$]',
+  "w3rd" : 'Third mom. of $w$ [m$^3$ s$^{-3}$]',
+  "sat_RH" : 'Supersaturation [\%]',
+  "rad_flx" : 'radiative flux [W m$^{-2}$]',
+  "cl_nc_zoom" : '$N_c$ [cm$^{-3}$]',
+  "base_prflux_vs_clhght" : 'pr. flux at column base [W m$^{-2}$]',
+  "base_prflux_vs_clhght number of occurances" : 'number of occurances',
+  "er" : 'Entrainment rate [cm s$^{-1}$]',
+  "wvarmax" : 'Max. $w$ variance [m$^{2}$ s$^{-2}$]',
+  "cloud_base" : 'Cloud base height [m]',
+  "gccn_rw_cl" : '$<r>$ of GCCN droplets (cloudy cells) [um]',
+  "non_gccn_rw_cl" : '$<r>$ of CCN droplets (cloudy cells) [um]',
+  "clb_bigrain_mean_rd" : '$<r_d>$ of (r$>$40um) @ clbase [m]',
+  "clb_bigrain_mean_kappa" : '$\kappa$ of (r$>$40um) @ clbase',
+  "clb_bigrain_mean_conc" : 'conc. of (r$>$40um) @ clbase [1/cc]',
+  "clb_bigrain_mean_inclt" : 'time since act. of (r$>$40um) @ clbase [s]',
+  "clb_bigrain_mean_gccn_fraction" : 'frac. of (r$>$40um) on gccn @ clbase'
+}
+
+
+labeldict = {
+ 0 : "(a)",
+ 1 : "(b)",
+ 2 : "(c)",
+ 3 : "(d)",
+ 4 : "(e)",
+ 5 : "(f)",
+ 6 : "(g)",
+ 7 : "(h)",
+ 8 : "(i)",
+ 9 : "(j)"
+}

Matplotlib_common/plot_profs.py

@@ -0,0 +1,44 @@
+import numpy as np
+from sys import argv
+
+from latex_labels import var_labels
+from read_UWLCM_arrays import *
+
+def plot_profiles(var_list, plot_iter, nplotx, nploty, axarr, xscaledict, yscaledict, xlimdict, ylimdict, suffix='', ylabel='', file_names=[], file_labels=[]):
+  # if file names are not defined, read them and labels from command line
+  if len(file_names)==0:
+    file_no = np.arange(1, len(argv)-1 , 2)
+    print file_no
+    for no in file_no:
+      file_names.append(argv[no]+suffix)
+      file_labels.append(argv[no+1])
+
+  for var in var_list:
+    label_counter = 0
+    for file_name in file_names:
+      print file_name
+    #  try:
+      profiles_file = open(file_name, "r")
+      my_pos = read_my_var(profiles_file, "position")
+      my_res = read_my_var(profiles_file, var)
+      profiles_file.close()
+
+      # remove artificial values of cloudy-cell variables near the surface due to incorrect cloudiness mask (e.g. RICO mask used in DYCOMS)
+      if(var == "gccn_rw_cl" or var == "non_gccn_rw_cl" or var == "cl_nc"):
+        my_res[0:10] = 0
+
+      linestyles = ['--', '-.', ':']
+      dashList = [(3,1),(1,1),(4,1,1,1),(4,2)] 
+
+      if(var == "base_prflux_vs_clhght" or var == "base_prflux_vs_clhght number of occurances"):
+        plot_my_array(axarr, plot_iter, my_res, my_pos, nploty, xlabel=var_labels[var], ylabel="cloudy column height [m]", varlabel=file_labels[label_counter], dashes = dashList[label_counter % len(dashList)], xscale=xscaledict[var], yscale=yscaledict[var], xlim=xlimdict[var], ylim=ylimdict[var])
+      else:
+        plot_my_array(axarr, plot_iter, my_res, my_pos, nploty, xlabel=var_labels[var], ylabel=ylabel, varlabel=file_labels[label_counter], dashes = dashList[label_counter % len(dashList)], xscale=xscaledict[var], yscale=yscaledict[var], xlim=xlimdict[var], ylim=ylimdict[var])
+   #   except:
+   #     print 'error opening file: ', file_name
+   #     my_pos = 0
+   #     my_res = 0
+      label_counter = label_counter+1
+    plot_iter += 1
+  return plot_iter
+

Matplotlib_common/plot_series.py

@@ -0,0 +1,33 @@
+import numpy as np
+from sys import argv
+
+from latex_labels import var_labels
+from read_UWLCM_arrays import *
+
+def plot_series(var_list, plot_iter, nplotx, nploty, axarr, xscaledict, yscaledict, xlimdict, ylimdict, show_bin=False, suffix='', xlabel='', ylabeldict=var_labels, file_names=[], file_labels=[]):
+  # if file names are not defined, read them and labels from command line
+  if len(file_names)==0:
+    file_no = np.arange(1, len(argv)-1 , 2)
+    for no in file_no:
+      file_names.append(argv[no] + suffix)
+      file_labels.append(argv[no+1])
+
+  for var in var_list:
+    label_counter=0
+    for file_name in file_names:
+      print file_name, var
+      series_file = open(file_name, "r")
+      my_times = read_my_var(series_file, "position")
+      my_res = read_my_var(series_file, var)
+
+      # rescale time to hours
+      my_times = my_times / 3600.
+
+      series_file.close()
+
+      linestyles = ['--', '-.', ':']
+      dashList = [(3,1),(1,1),(4,1,1,1),(4,2)] 
+      plot_my_array(axarr, plot_iter, my_times, my_res, nploty, xlabel=xlabel, ylabel=ylabeldict[var], varlabel=file_labels[label_counter], dashes = dashList[label_counter % len(dashList)], xscale=xscaledict[var], yscale=yscaledict[var], xlim=xlimdict[var], ylim=ylimdict[var])
+      label_counter+=1
+    plot_iter = plot_iter + 1
+  return plot_iter

Matplotlib_common/read_UWLCM_arrays.py

@@ -0,0 +1,54 @@
+import numpy as np
+
+def read_my_array(file_obj):
+  arr_name = file_obj.readline()
+  file_obj.readline() # discarded line with size of the array
+  line = file_obj.readline()
+  line = line.split(" ")
+  del line[0]
+  del line[len(line)-1]
+  arr = map(float,line)
+  return np.array(arr), arr_name
+
+def read_my_var(file_obj, var_name):
+  file_obj.seek(0)
+  while True:
+    arr, name = read_my_array(file_obj)
+    if(str(name).strip() == str(var_name).strip()):
+      break
+  return arr
+
+#def plot_my_array(axarr, plot_iter, time, val, nploty, xlabel=None, ylabel=None, varlabel=None , linestyle='--', dashes=(5,2), xlim=None, ylim=None, xscale="linear", yscale="linear"):
+#  x = int(plot_iter / nploty)
+#  y = plot_iter % nploty
+#  if varlabel != None:
+#    axarr[x, y].plot(time, val, label=varlabel, linestyle=linestyle, linewidth=1, dashes=dashes)
+#  else:
+#    axarr[x, y].plot(time, val, linestyle=linestyle, linewidth=1, dashes=dashes)
+#  if xlabel:
+#    axarr[x, y].set_xlabel(xlabel)
+#  if ylabel:
+#    axarr[x, y].set_ylabel(ylabel)
+#  if xlim:
+#    axarr[x, y].set_xlim(xlim)
+#  if ylim:
+#    axarr[x, y].set_ylim(ylim)
+#  axarr[x, y].set_xscale(xscale)
+#  axarr[x, y].set_yscale(yscale)
+
+def plot_my_array(axarr, plot_iter, time, val, nploty, xlabel=None, ylabel=None, varlabel=None , linestyle='--', dashes=(5,2), xlim=None, ylim=None, xscale="linear", yscale="linear"):
+  axflat = axarr.flatten() # flattening makes it work with 1d axarr
+  if varlabel != None:
+    axflat[plot_iter].plot(time, val, label=varlabel, linestyle=linestyle, linewidth=1, dashes=dashes)
+  else:
+    axflat[plot_iter].plot(time, val, linestyle=linestyle, linewidth=1, dashes=dashes)
+  if xlabel:
+    axflat[plot_iter].set_xlabel(xlabel)
+  if ylabel:
+    axflat[plot_iter].set_ylabel(ylabel)
+  if xlim:
+    axflat[plot_iter].set_xlim(xlim)
+  if ylim:
+    axflat[plot_iter].set_ylim(ylim)
+  axflat[plot_iter].set_xscale(xscale)
+  axflat[plot_iter].set_yscale(yscale)