Black reformat

CPyS/B.py

@@ -66,7 +66,7 @@ def area_weights(field):
     return w
 
 
-def B(th, geopt, SH=False, names=["snap_z900", "snap_z600"]): # TODO: Useless?
+def B(th, geopt, SH=False, names=["snap_z900", "snap_z600"]):  # TODO: Useless?
     """
     Computes the B parameter for a point, with the corresponding snapshot of geopt at 600hPa and 900hPa
 

CPyS/CPS.py

@@ -4,12 +4,13 @@
 import pandas as pd
 import numpy as np
 
+
 def compute_CPS_parameters(
     tracks,
     geopt,
     geopt_name="snap_zg",
     plev_name="level",
-    verbose = True,
+    verbose=True,
 ):
     """
     Computes the three (+ theta) Hart parameters for all the points in tracks.
@@ -29,30 +30,41 @@ def compute_CPS_parameters(
 
     # Curate input
     ## Snapshots
-    geopt = geopt.rename({plev_name: "plev"}) # Change vertical coordinate name
+    geopt = geopt.rename({plev_name: "plev"})  # Change vertical coordinate name
     geopt = geopt.where(np.abs(geopt.snap_zg) < 1e10)
-	
+
     # tracks
-    if "time" not in tracks.columns :
-        tracks["time"] = pd.to_datetime(tracks.year.astype(str) + '-' + tracks.month.astype(str) + '-' + tracks.day.astype(str) + ' ' + tracks.hour.astype(str) + ':00:00')
-
+    if "time" not in tracks.columns:
+        tracks["time"] = pd.to_datetime(
+            tracks.year.astype(str)
+            + "-"
+            + tracks.month.astype(str)
+            + "-"
+            + tracks.day.astype(str)
+            + " "
+            + tracks.hour.astype(str)
+            + ":00:00"
+        )
+
     # 1/ B computation
-    if verbose: print("Computing B...")
+    if verbose:
+        print("Computing B...")
     ## Select 900 & 600 hPa levels
     z900, z600 = (
         geopt[geopt_name].sel(plev=900e2, method="nearest"),
         geopt[geopt_name].sel(plev=600e2, method="nearest"),
     )
-    if verbose: print(
-        "Level "
-        + str(z900.plev.values)
-        + " is taken for 900hPa"
-        + "\n"
-        + "Level "
-        + str(z600.plev.values)
-        + " is taken for 600hPa"
-        + "\n"
-    )
+    if verbose:
+        print(
+            "Level "
+            + str(z900.plev.values)
+            + " is taken for 900hPa"
+            + "\n"
+            + "Level "
+            + str(z600.plev.values)
+            + " is taken for 600hPa"
+            + "\n"
+        )
 
     ## theta computation
     if "theta" not in tracks.columns:
@@ -64,7 +76,8 @@ def compute_CPS_parameters(
     )
 
     # 2/ VTL & VTU computation
-    if verbose: print("Computing VTL & VTU...")
+    if verbose:
+        print("Computing VTL & VTU...")
     geopt = geopt.sortby("plev", ascending=False)
     VTL, VTU = VT(geopt, name=geopt_name)
 

CPyS/VT.py

@@ -1,5 +1,6 @@
 import numpy as np
 
+
 def VT(geopt, name="snap_zg"):
     """
     Parameters

CPyS/plot.py

@@ -1,27 +1,28 @@
 import matplotlib.pyplot as plt
 
-def plot_CPS(tracks, title = ""):
-    fig, axs = plt.subplots(1,2, figsize = [10,5])
+
+def plot_CPS(tracks, title=""):
+    fig, axs = plt.subplots(1, 2, figsize=[10, 5])
     fig.suptitle(title)
 
     # Left plot (B vs. VTL)
     ## Data
-    axs[0].plot(tracks.VTL, tracks.B, marker = "o", color = "k")
+    axs[0].plot(tracks.VTL, tracks.B, marker="o", color="k")
     ## x-axis
-    axs[0].axvline(x=0, color = "k", alpha = 0.5, linestyle = "--", linewidth = 1)
+    axs[0].axvline(x=0, color="k", alpha=0.5, linestyle="--", linewidth=1)
     axs[0].set_xlabel("$-V_T^L$ / m")
     ## y-axis
-    axs[0].axhline(y=10, color = "k", alpha = 0.5, linestyle = "--", linewidth = 1)
+    axs[0].axhline(y=10, color="k", alpha=0.5, linestyle="--", linewidth=1)
     axs[0].set_ylabel("B / m")
 
-    #Right plot (VTU vs. VTL)
+    # Right plot (VTU vs. VTL)
     ## Data
-    axs[1].plot(tracks.VTL, tracks.VTU, marker = "o", color = "k")
+    axs[1].plot(tracks.VTL, tracks.VTU, marker="o", color="k")
     ## x-axis
-    axs[1].axvline(x=0, color = "k", alpha = 0.5, linestyle = "--", linewidth = 1)
+    axs[1].axvline(x=0, color="k", alpha=0.5, linestyle="--", linewidth=1)
     axs[1].set_xlabel("$-V_T^L$ / m")
     ## y-axis
-    axs[1].axhline(y=0, color = "k", alpha = 0.5, linestyle = "--", linewidth = 1)
+    axs[1].axhline(y=0, color="k", alpha=0.5, linestyle="--", linewidth=1)
     axs[1].set_ylabel("$-V_T^U$ / m")
 
     plt.tight_layout()

CPyS/theta.py

@@ -1,6 +1,7 @@
 import numpy as np
 import pyproj
 
+
 def theta(x0=120, x1=130, y0=12, y1=10):
     """
     Computes the angular direction between to points.
@@ -18,8 +19,8 @@ def theta(x0=120, x1=130, y0=12, y1=10):
     th (float): The directional angle between the current and the next point, in degrees.
         0 for eastward, 90 for northward.
     """
-    geodesic = pyproj.Geod(ellps='WGS84')
-    fwd_azimuth,back_azimuth,distance = geodesic.inv(x0, y0, x1, y1)
+    geodesic = pyproj.Geod(ellps="WGS84")
+    fwd_azimuth, back_azimuth, distance = geodesic.inv(x0, y0, x1, y1)
     return -1 * (fwd_azimuth - 90) % 360
 
 

build/lib/CPyS/B.py

@@ -66,7 +66,7 @@ def area_weights(field):
     return w
 
 
-def B(th, geopt, SH=False, names=["snap_z900", "snap_z600"]): # TODO: Useless?
+def B(th, geopt, SH=False, names=["snap_z900", "snap_z600"]):  # TODO: Useless?
     """
     Computes the B parameter for a point, with the corresponding snapshot of geopt at 600hPa and 900hPa
 

build/lib/CPyS/CPS.py

@@ -1,13 +1,16 @@
 from .theta import theta_multitrack
 from .B import B_vector
 from .VT import VT
+import pandas as pd
+import numpy as np
 
 
 def compute_CPS_parameters(
     tracks,
     geopt,
     geopt_name="snap_zg",
     plev_name="level",
+    verbose=True,
 ):
     """
     Computes the three (+ theta) Hart parameters for all the points in tracks.
@@ -25,25 +28,43 @@ def compute_CPS_parameters(
     tracks (pd.DataFrame): The set of TC points with four new columns corresponding to the parameters
     """
 
-    geopt = geopt.rename({plev_name: "plev"}) # Change vertical coordinate name
+    # Curate input
+    ## Snapshots
+    geopt = geopt.rename({plev_name: "plev"})  # Change vertical coordinate name
+    geopt = geopt.where(np.abs(geopt.snap_zg) < 1e10)
+
+    # tracks
+    if "time" not in tracks.columns:
+        tracks["time"] = pd.to_datetime(
+            tracks.year.astype(str)
+            + "-"
+            + tracks.month.astype(str)
+            + "-"
+            + tracks.day.astype(str)
+            + " "
+            + tracks.hour.astype(str)
+            + ":00:00"
+        )
 
     # 1/ B computation
-    print("Computing B...")
+    if verbose:
+        print("Computing B...")
     ## Select 900 & 600 hPa levels
     z900, z600 = (
         geopt[geopt_name].sel(plev=900e2, method="nearest"),
         geopt[geopt_name].sel(plev=600e2, method="nearest"),
     )
-    print(
-        "Level "
-        + str(z900.plev.values)
-        + " is taken for 900hPa"
-        + "\n"
-        + "Level "
-        + str(z600.plev.values)
-        + " is taken for 600hPa"
-        + "\n"
-    )
+    if verbose:
+        print(
+            "Level "
+            + str(z900.plev.values)
+            + " is taken for 900hPa"
+            + "\n"
+            + "Level "
+            + str(z600.plev.values)
+            + " is taken for 600hPa"
+            + "\n"
+        )
 
     ## theta computation
     if "theta" not in tracks.columns:
@@ -55,7 +76,8 @@ def compute_CPS_parameters(
     )
 
     # 2/ VTL & VTU computation
-    print("Computing VTL & VTU...")
+    if verbose:
+        print("Computing VTL & VTU...")
     geopt = geopt.sortby("plev", ascending=False)
     VTL, VTU = VT(geopt, name=geopt_name)
 

build/lib/CPyS/VT.py

@@ -1,5 +1,6 @@
 import numpy as np
 
+
 def VT(geopt, name="snap_zg"):
     """
     Parameters

build/lib/CPyS/plot.py

@@ -1,27 +1,28 @@
 import matplotlib.pyplot as plt
 
-def plot_CPS(tracks, title = ""):
-    fig, axs = plt.subplots(1,2, figsize = [10,5])
+
+def plot_CPS(tracks, title=""):
+    fig, axs = plt.subplots(1, 2, figsize=[10, 5])
     fig.suptitle(title)
 
     # Left plot (B vs. VTL)
     ## Data
-    axs[0].plot(tracks.VTL, tracks.B, marker = "o", color = "k")
+    axs[0].plot(tracks.VTL, tracks.B, marker="o", color="k")
     ## x-axis
-    axs[0].axvline(x=0, color = "k", alpha = 0.5, linestyle = "--", linewidth = 1)
+    axs[0].axvline(x=0, color="k", alpha=0.5, linestyle="--", linewidth=1)
     axs[0].set_xlabel("$-V_T^L$ / m")
     ## y-axis
-    axs[0].axhline(y=10, color = "k", alpha = 0.5, linestyle = "--", linewidth = 1)
+    axs[0].axhline(y=10, color="k", alpha=0.5, linestyle="--", linewidth=1)
     axs[0].set_ylabel("B / m")
 
-    #Right plot (VTU vs. VTL)
+    # Right plot (VTU vs. VTL)
     ## Data
-    axs[1].plot(tracks.VTL, tracks.VTU, marker = "o", color = "k")
+    axs[1].plot(tracks.VTL, tracks.VTU, marker="o", color="k")
     ## x-axis
-    axs[1].axvline(x=0, color = "k", alpha = 0.5, linestyle = "--", linewidth = 1)
+    axs[1].axvline(x=0, color="k", alpha=0.5, linestyle="--", linewidth=1)
     axs[1].set_xlabel("$-V_T^L$ / m")
     ## y-axis
-    axs[1].axhline(y=0, color = "k", alpha = 0.5, linestyle = "--", linewidth = 1)
+    axs[1].axhline(y=0, color="k", alpha=0.5, linestyle="--", linewidth=1)
     axs[1].set_ylabel("$-V_T^U$ / m")
 
     plt.tight_layout()

build/lib/CPyS/theta.py

@@ -1,4 +1,5 @@
 import numpy as np
+import pyproj
 
 
 def theta(x0=120, x1=130, y0=12, y1=10):
@@ -16,20 +17,11 @@ def theta(x0=120, x1=130, y0=12, y1=10):
     Returns
     -------
     th (float): The directional angle between the current and the next point, in degrees.
+        0 for eastward, 90 for northward.
     """
-    v = [x1 - x0, y1 - y0]  # Vector corresponding to the track movement
-    if v != [0, 0]:  # If not stagnating
-        cos = (x1 - x0) / (np.linalg.norm(v))  # Cosinus with eastward vector [1,0]
-        if cos == -1:
-            th = 180
-        else:
-            th = np.sign(y1 - y0) * np.arccos(cos) * 180 / np.pi
-    else:  # If stagnating: Set to NaN
-        th = np.nan
-    if th < 0:
-        th += 360
-        # Make all angles between 0 and 360
-    return th
+    geodesic = pyproj.Geod(ellps="WGS84")
+    fwd_azimuth, back_azimuth, distance = geodesic.inv(x0, y0, x1, y1)
+    return -1 * (fwd_azimuth - 90) % 360
 
 
 def theta_track(lon, lat):