From ae14194fed3888e14b5b3bfa777a8d15e39527ff Mon Sep 17 00:00:00 2001
From: HenryTux <h.schoeller@fu-berlin.de>
Date: Fri, 12 Jan 2024 17:19:57 +0100
Subject: [PATCH] updates

---
 calc_E.py | 158 ++++++++++++++++++++++++++++--------------------------
 calc_E.sh |   6 ++-
 main.py   |  12 +++--
 3 files changed, 94 insertions(+), 82 deletions(-)

diff --git a/calc_E.py b/calc_E.py
index 6846bdd..58acdf9 100644
--- a/calc_E.py
+++ b/calc_E.py
@@ -29,104 +29,110 @@ if argv[0]=='':
 else:
     date_0 = datetime(int(argv[2]),int(argv[3]),int(argv[4]),int(argv[5]))
     trapath = argv[1]
+    print(date_0)
 
 fname = "traj" + date_0.strftime('_%Y%m%d_%H') + ".npy"
 trajs = np.load(trapath + fname)
 
+for a in [0,1]:
+    if a == 0:
 
-t_sel = np.arange(0,int(trajs.shape[1]/2)+1)
-
-z_coord_type = "p"
-
-n = 1
+        t_sel = np.arange(0,int(trajs.shape[1]/2)+1)
+    else:
+        t_sel = np.arange(int(trajs.shape[1]/2),trajs.shape[1])
 
-path = "./" + date_0.strftime("dat_%Y%m%d_%H") + "/" + str(n) + "/"
-if not exists(path):
-    makedirs(path)
 
-print(n)
+    z_coord_type = "p"
 
-Lon = trajs['lon'][::n]
-Lat = trajs['lat'][::n]
-P = trajs['p'][::n]
-U = trajs['U'][::n]
-V = trajs['V'][::n]
-Omg = trajs['OMEGA'][::n]
-T = trajs['T'][::n]
+    n = 1
 
-# %% Initials
+    path = "./" + date_0.strftime("dat_%Y%m%d_%H") + "/" + str(n) + "/"
+    if not exists(path):
+        makedirs(path)
 
-H = cc.ptoh(P, 1013.25, 8.435) # standard atmosphere
+    print(n)
 
-# 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))
+    Lon = trajs['lon'][::n]
+    Lat = trajs['lat'][::n]
+    P = trajs['p'][::n]
+    U = trajs['U'][::n]
+    V = trajs['V'][::n]
+    Omg = trajs['OMEGA'][::n]
+    T = trajs['T'][::n]
 
-D = list((list(), list(), list()))
+    # %% Initials
 
-epsilon_opt = np.array([20000, 50000, 100000, 200000, 500000])
-bound_meth = "$\\alpha$"
-alpha = 1e-3
+    H = cc.ptoh(P, 1013.25, 8.435) # standard atmosphere
 
-for epsilon in epsilon_opt:
+    # 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))
 
-    print(epsilon)
+    D = list((list(), list(), list()))
 
-    if z_coord_type=="p":
-        x, y, z = cc.coord_trans(Lon, Lat, P, "None", 1,
-                                 proj="stereo")
+    epsilon_opt = np.array([20000, 50000, 100000, 200000, 500000])
+    bound_meth = "$\\alpha$"
+    alpha = 1e-3
 
-        if (bound_meth == "Concave"):
-            boundscs, hullscs = dd.io_bounds(path + "stereop",# + "norm" ,
-                                         "opt. $\\alpha$", x, y, z, None)
-            E = dd.io_eigen((path + "spnopt" + str(t_sel[0]) +
-                         str(t_sel[-1])),
-                        Lon.shape[0], t_sel[0], t_sel[-1], D[0], D[1], D[2],
-                        epsilon, boundscs, 20, Lon, Lat, P, True, k_p, "p")
+    for epsilon in epsilon_opt:
 
-        elif (bound_meth == "Convex"):
-            boundscs, hullscs = dd.io_bounds(path + "stereop",#  + "norm" ,
-                                             "Convex", x, y, z, None)
-            E = dd.io_eigen((path + "spncvx" + str(t_sel[0]) +
-                         str(t_sel[-1])),
-                        Lon.shape[0], t_sel[0], t_sel[-1], D[0], D[1], D[2],
-                        epsilon, boundscs, 20, Lon, Lat, P, True, k_p, "p")
-
-        elif (bound_meth == "$\\alpha$"):
-            # for alpha in np.array([1e-4, 2e-4, 5e-4, 0.001, 0.002, 0.005, 0.01]):
-            boundscs, hullscs = dd.io_bounds(path + "stereop", "$\\alpha$",
-                                             x, y, z, alpha)
-            E = dd.io_eigen((path + "spnalp" + str(alpha) + str(t_sel[0]) +
-                     str(t_sel[-1])),
-                    Lon.shape[0], t_sel[0], t_sel[-1], D[0], D[1], D[2],
-                    epsilon, boundscs, 20, Lon, Lat, P, True, k_p, "p",
-                    force_calc=True)
+        print(epsilon)
 
-    else:
-        x, y, z = cc.coord_trans(Lon, Lat, H, "None", 1,
-                                 proj="stereo")
+        if z_coord_type=="p":
+            x, y, z = cc.coord_trans(Lon, Lat, P, "None", 1,
+                                     proj="stereo")
 
-        if (bound_meth == "Concave"):
-            boundscs, hullscs = dd.io_bounds(path + "stereoh",# + "norm" ,
-                                         "opt. $\\alpha$", x, y, z, None)
-            E = dd.io_eigen((path + "shnopt" + str(t_sel[0]) +
-                         str(t_sel[-1])),
-                        Lon.shape[0], t_sel[0], t_sel[-1], D[0], D[1], D[2],
-                        epsilon, boundscs, 20, Lon, Lat, H, True, k_h, "H")
+            if (bound_meth == "Concave"):
+                boundscs, hullscs = dd.io_bounds(path + "stereop",# + "norm" ,
+                                             "opt. $\\alpha$", x, y, z, None)
+                E = dd.io_eigen((path + "spnopt" + str(t_sel[0]) +
+                             str(t_sel[-1])),
+                            Lon.shape[0], t_sel[0], t_sel[-1], D[0], D[1], D[2],
+                            epsilon, boundscs, 20, Lon, Lat, P, True, k_p, "p")
 
-        elif (bound_meth == "Convex"):
-            boundscs, hullscs = dd.io_bounds(path + "stereoh",#  + "norm" ,
-                                             "Convex", x, y, z, None)
-            E = dd.io_eigen((path + "shncvx" + str(t_sel[0]) +
-                         str(t_sel[-1])),
-                        Lon.shape[0], t_sel[0], t_sel[-1], D[0], D[1], D[2],
-                        epsilon, boundscs, 20, Lon, Lat, H, True, k_h, "H")
+            elif (bound_meth == "Convex"):
+                boundscs, hullscs = dd.io_bounds(path + "stereop",#  + "norm" ,
+                                                 "Convex", x, y, z, None)
+                E = dd.io_eigen((path + "spncvx" + str(t_sel[0]) +
+                             str(t_sel[-1])),
+                            Lon.shape[0], t_sel[0], t_sel[-1], D[0], D[1], D[2],
+                            epsilon, boundscs, 20, Lon, Lat, P, True, k_p, "p")
 
-        elif (bound_meth == "$\\alpha$"):
-            for alpha in np.array([1e-4, 2e-4, 5e-4, 0.001, 0.002, 0.005, 0.01]):
-                boundscs, hullscs = dd.io_bounds(path + "stereoh", "$\\alpha$",
+            elif (bound_meth == "$\\alpha$"):
+                # for alpha in np.array([1e-4, 2e-4, 5e-4, 0.001, 0.002, 0.005, 0.01]):
+                boundscs, hullscs = dd.io_bounds(path + "stereop", "$\\alpha$",
                                                  x, y, z, alpha)
-                E = dd.io_eigen((path + "shnalp" + str(alpha) + str(t_sel[0]) +
+                E = dd.io_eigen((path + "spnalp" + str(alpha) + str(t_sel[0]) +
                          str(t_sel[-1])),
                         Lon.shape[0], t_sel[0], t_sel[-1], D[0], D[1], D[2],
-                        epsilon, boundscs, 20, Lon, Lat, H, True, k_h, "H")
\ No newline at end of file
+                        epsilon, boundscs, 20, Lon, Lat, P, True, k_p, "p",
+                        force_calc=True)
+
+        else:
+            x, y, z = cc.coord_trans(Lon, Lat, H, "None", 1,
+                                     proj="stereo")
+
+            if (bound_meth == "Concave"):
+                boundscs, hullscs = dd.io_bounds(path + "stereoh",# + "norm" ,
+                                             "opt. $\\alpha$", x, y, z, None)
+                E = dd.io_eigen((path + "shnopt" + str(t_sel[0]) +
+                             str(t_sel[-1])),
+                            Lon.shape[0], t_sel[0], t_sel[-1], D[0], D[1], D[2],
+                            epsilon, boundscs, 20, Lon, Lat, H, True, k_h, "H")
+
+            elif (bound_meth == "Convex"):
+                boundscs, hullscs = dd.io_bounds(path + "stereoh",#  + "norm" ,
+                                                 "Convex", x, y, z, None)
+                E = dd.io_eigen((path + "shncvx" + str(t_sel[0]) +
+                             str(t_sel[-1])),
+                            Lon.shape[0], t_sel[0], t_sel[-1], D[0], D[1], D[2],
+                            epsilon, boundscs, 20, Lon, Lat, H, True, k_h, "H")
+
+            elif (bound_meth == "$\\alpha$"):
+                for alpha in np.array([1e-4, 2e-4, 5e-4, 0.001, 0.002, 0.005, 0.01]):
+                    boundscs, hullscs = dd.io_bounds(path + "stereoh", "$\\alpha$",
+                                                     x, y, z, alpha)
+                    E = dd.io_eigen((path + "shnalp" + str(alpha) + str(t_sel[0]) +
+                             str(t_sel[-1])),
+                            Lon.shape[0], t_sel[0], t_sel[-1], D[0], D[1], D[2],
+                            epsilon, boundscs, 20, Lon, Lat, H, True, k_h, "H")
\ No newline at end of file
diff --git a/calc_E.sh b/calc_E.sh
index eb4436c..cff4aa9 100755
--- a/calc_E.sh
+++ b/calc_E.sh
@@ -1,8 +1,8 @@
 #!/bin/bash
 
-#SBATCH --array=1-24                                                                                   
+#SBATCH --array=9-23                                                                                   
 
-year=2016
+year=2017
 TRAPATH=/net/scratch/schoelleh96/WP2/WP2.1/LAGRANTO/wp21/era5/traj
 CSPATH=/net/scratch/schoelleh96/WP2/WP2.1/LAGRANTO/csstandalone
 
@@ -16,4 +16,6 @@ M1=${filename:9:2}
 D1=${filename:11:2}
 H1=${filename:14:2}
 
+source /home/schoelleh96/Applications/envs/wp21env/bin/activate
+
 python ${CSPATH}/calc_E.py ${TRAPATH}/${year}/ ${Y1} ${M1} ${D1} ${H1}
diff --git a/main.py b/main.py
index 1d40f41..37c3d2a 100644
--- a/main.py
+++ b/main.py
@@ -32,12 +32,12 @@ n = 1
 
 force_calc = False
 
-date_0 = datetime(2017,1,25,0)
+date_0 = datetime(2016,5,2,0)
 fname = "../wp21/era5/traj/" + date_0.strftime('%Y/traj_%Y%m%d_%H') + ".npy"
 
 B = dd.get_block_dat(date_0)
 if date_0.strftime("%Y") == "2016":
-    coasts = dd.get_coasts(10, 90, -60, -240)
+    coasts = dd.get_coasts(10, 90, 0, -360)
 elif date_0.strftime("%Y") == "2017":
     coasts = dd.get_coasts(10, 90, 120, -60)
 
@@ -60,6 +60,7 @@ if direc=="whole":
 elif direc=="forward":
     t_sel = np.arange(int(trajs.shape[1]/2),trajs.shape[1])
     plot_0 = 0
+    B = list(B[12:,:,:], B[1], B[2])
 
 elif direc=="backward":
     t_sel = np.arange(0,int(trajs.shape[1]/2)+1)
@@ -488,16 +489,19 @@ if True:
 
     keys, inv = np.unique(kclust.labels_, return_inverse=True)
     vals = np.array([colors[key] for key in keys])
+    # to remove Null-Level-Set cluster (assuming its cluster 5):
+    # dat = dat[inv!=5, :]
+    # inv = inv[inv !=5]
     pp.plot_spaghetti(ax_spag, dat, alpha_spag, colors=colors, inv=inv)
 
-    v = pp.draw_vert_t(ax_spag, 0, lon.shape[1], t_slider2.val)
+    v = pp.draw_vert_t(ax_spag, t_sel[0]-72, t_sel[-1] -72, t_slider2.val -72)
     plt.show()
 
     def var_changed(val):
         dat = trajs[var_check.value_selected][::n, t_sel]
         ax_spag.clear()
         pp.plot_spaghetti(ax_spag, dat, alpha_spag, colors=colors, inv=inv)
-        pp.draw_vert_t(ax_spag, 0, lon.shape[1], t_slider2.val)
+        pp.draw_vert_t(ax_spag, t_sel[0]-72, t_sel[-1] -72, t_slider2.val -72)
         fig_spag.canvas.draw_idle()
 
     var_check.on_clicked(var_changed)
-- 
GitLab