updates

calc_E.py

@@ -16,12 +16,9 @@ chdir(dirname(abspath(__file__)))
 from sys import argv, path
 import numpy as np
 path.append("../wp21/pyscripts/LIB")
-import plot as pp
 import calc as cc
 import data as dd
-import matplotlib.pyplot as plt
 from datetime import datetime
-import itertools
 
 # %% Get the data
 
@@ -37,7 +34,7 @@ fname = "traj" + date_0.strftime('_%Y%m%d_%H') + ".npy"
 trajs = np.load(trapath + fname)
 
 
-t_sel = np.arange(0,trajs.shape[1])
+t_sel = np.arange(0,int(trajs.shape[1]/2)+1)
 
 z_coord_type = "p"
 
@@ -61,13 +58,13 @@ T = trajs['T'][::n]
 
 H = cc.ptoh(P, 1013.25, 8.435) # standard atmosphere
 
-k_p = cc.calc_k(U/1000, V/1000, Omg/100)
+# k_p = cc.calc_k(U/1000, V/1000, Omg/100)
 k_p=15
 # k_h = cc.calc_k(U, V, cc.omg2w(Omg/100, T, P))
 
 D = list((list(), list(), list()))
 
-epsilon_opt = np.array([25000, 50000, 100000, 250000, 500000])
+epsilon_opt = np.array([20000, 50000, 100000, 200000, 500000])
 bound_meth = "$\\alpha$"
 alpha = 1e-3
 

calc_bound.py

@@ -16,10 +16,8 @@ chdir(dirname(abspath(__file__)))
 from sys import argv, path, stdout
 import numpy as np
 path.append("../wp21/pyscripts/LIB")
-import plot as pp
 import calc as cc
 import data as dd
-import matplotlib.pyplot as plt
 from datetime import datetime
 import itertools
 
@@ -81,8 +79,8 @@ elif coord == "x, y, p (stereo)":
     x, y, z = cc.coord_trans(Lon, Lat, P, "None", 1, proj="stereo")
     d_path = path + "stereop"
     if (not norm):
-        z = z * cc.calc_k(trajs['U']/1000, trajs['V']/1000,
-                          trajs['OMEGA']/100)
+        z = z * 15 #cc.calc_k(trajs['U']/1000, trajs['V']/1000,
+                          # trajs['OMEGA']/100)
 elif coord == "x, y, h (stereo)":
     x, y, z = cc.coord_trans(Lon, Lat, P, "std_atm", 1, proj="stereo")
     d_path = path + "stereoh"

calc_dist.py

@@ -16,12 +16,9 @@ chdir(dirname(abspath(__file__)))
 from sys import argv, path
 import numpy as np
 path.append("../wp21/pyscripts/LIB")
-import plot as pp
 import calc as cc
 import data as dd
-import matplotlib.pyplot as plt
 from datetime import datetime
-import itertools
 
 # %% Get the data
 

main.py

@@ -32,11 +32,14 @@ n = 1
 
 force_calc = False
 
-date_0 = datetime(2016,4,30,12)
+date_0 = datetime(2017,1,25,0)
 fname = "../wp21/era5/traj/" + date_0.strftime('%Y/traj_%Y%m%d_%H') + ".npy"
 
 B = dd.get_block_dat(date_0)
-coasts = dd.get_coasts(0, 90, -180, 180)
+if date_0.strftime("%Y") == "2016":
+    coasts = dd.get_coasts(10, 90, -60, -240)
+elif date_0.strftime("%Y") == "2017":
+    coasts = dd.get_coasts(10, 90, 120, -60)
 
 #fname = "traj" + date_0.strftime('_%Y%m%d_%H') + ".npy"
 path = "./" + date_0.strftime("dat_%Y%m%d_%H") + "/" + str(n) + "/"
@@ -46,7 +49,7 @@ if not exists(path):
 trajs = np.load(fname)
 
 if argv[0]=='':
-     direc = 'whole'
+     direc = 'backward'
 else:
     direc = argv[1]
 
@@ -59,8 +62,8 @@ elif direc=="forward":
     plot_0 = 0
 
 elif direc=="backward":
-    t_sel = np.arange(0,int(trajs.shape[1]/2))
-    plot_0 = -1
+    t_sel = np.arange(0,int(trajs.shape[1]/2)+1)
+    plot_0 = int(trajs.shape[1]/2)
 
 # t_sel = np.arange(73, 74)
 
@@ -86,14 +89,14 @@ elev = 45
 dist = 7
 t_i = 0
 alpha = 1e-3
-e = 100000
-epsilon_opt = np.array([25000, 50000, 100000, 250000, 500000])
+e = 20000
+epsilon_opt = np.array([20000, 50000, 100000, 200000, 500000])
 alph_opt = np.array([1e-4, 2e-4, 5e-4, 0.001, 0.002, 0.005, 0.01, 0.05,
                      0.1, 0.5, 1, 5, 10])
-N_k = 3
+N_k = 6
 # u and v is in m/s, omg in P/s; result is in km/hPa
-k_p = cc.calc_k(trajs['U']/1000, trajs['V']/1000, trajs['OMEGA']/100)
-
+# k_p = cc.calc_k(trajs['U']/1000, trajs['V']/1000, trajs['OMEGA']/100)
+k_p = 15
 # Coordinates for visualization
 X = lon
 Y = lat
@@ -103,6 +106,7 @@ k_h = cc.calc_k(trajs['U'], trajs['V'], cc.omg2w(trajs['OMEGA']/100,
                                                  trajs['T'], trajs['p']))
 
 # %% Boundary
+
 # if False:
 if argv[0] == '' or argv[1] == "Boundary":
 
@@ -307,11 +311,11 @@ if True:
     kclust = cc.kcluster_idx(E, N_k)
 
     # upper case letters for visualization
-    Xcs, Ycs, Zcs = cc.coord_trans(lon, lat, p, "None", k_p, proj="stereo")
+    Xcs, Ycs, Zcs = cc.coord_trans(lon, lat, p, "None", 1, proj="stereo")
     points = pp.plot_clustering(ax3dcs, Xcs, Ycs, Zcs, azim, elev, dist, t_i,
                                 c=kclust.labels_, plot_0=False,
-                                coord="stereo", mean_traj=True, block=B,
-                                coast=coasts)
+                                coord="stereo", mean_traj=True, block=None,
+                                coast=None)
     ax3dcs.invert_zaxis()
 
     # Add widget items
@@ -338,7 +342,7 @@ if True:
     null_ax = fig_cs.add_axes([0.02, 0.43, 0.1, 0.05])
     null_check = CheckButtons(null_ax, ["Plot start?"], [False])
     geo_ax = fig_cs.add_axes([0.1, 0.43, 0.1, 0.05])
-    geo_check = CheckButtons(geo_ax, ["Plot features?"], [True])
+    geo_check = CheckButtons(geo_ax, ["Plot features?"], [False])
     run_ax = fig_cs.add_axes([0.02, 0.35, 0.1, 0.05])
     run_check2 = Button(run_ax, "Calculate")
 
@@ -474,7 +478,7 @@ if True:
     ax_spag = fig_spag.add_subplot(111)
     fig_spag.subplots_adjust(left=0.2)
     var_ax = fig_spag.add_axes([0.025, 0.1, 0.1, 0.8])
-    var_check = RadioButtons(var_ax, trajs.dtype.names, 3)
+    var_check = RadioButtons(var_ax, trajs.dtype.names, 6)
     var_ax.set_title("Variable")
     dat = trajs[var_check.value_selected][::n, t_sel]