improved plots

Author: sdat2

fill-in.ipynb

@@ -1982,7 +1982,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 94,
+   "execution_count": 96,
    "id": "a10d439a",
    "metadata": {},
    "outputs": [
@@ -1992,25 +1992,25 @@
      "text": [
       "Interpolating specific humidity 'q' as well.\n",
       "t -79.06939697265625 30.296119589309058\n",
-      "'calculate_pi'  0.03452 s \n",
+      "'calculate_pi'  0.03005 s \n",
       "\n",
       "t -79.06939697265625 30.296119589309058\n",
-      "'calculate_pi'  0.03012 s \n",
+      "'calculate_pi'  0.02750 s \n",
       "\n",
       "t -79.06939697265625 29.628082275390625\n",
-      "'calculate_pi'  0.03337 s \n",
+      "'calculate_pi'  0.03104 s \n",
       "\n",
       "t -79.06939697265625 25.721405029296875\n",
-      "'calculate_pi'  0.02879 s \n",
+      "'calculate_pi'  0.02786 s \n",
       "\n",
       "t -79.06939697265625 17.55535888671875\n",
-      "'calculate_pi'  0.02905 s \n",
+      "'calculate_pi'  0.02669 s \n",
       "\n",
       "t -79.06939697265625 30.296119589309065\n",
-      "'calculate_pi'  0.02621 s \n",
+      "'calculate_pi'  0.02583 s \n",
       "\n",
       "t -79.06939697265625 30.296119589309058\n",
-      "'calculate_pi'  0.03303 s \n",
+      "'calculate_pi'  0.03278 s \n",
       "\n"
      ]
     },
@@ -2031,7 +2031,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 93,
+   "execution_count": 95,
    "id": "8483fb42",
    "metadata": {},
    "outputs": [
@@ -2041,25 +2041,25 @@
      "text": [
       "Interpolating specific humidity 'q' as well.\n",
       "t -69.27500915527344 31.399871826171875\n",
-      "'calculate_pi'  0.07401 s \n",
+      "'calculate_pi'  0.04714 s \n",
       "\n",
       "t -69.27500915527344 31.399871826171875\n",
-      "'calculate_pi'  0.05028 s \n",
+      "'calculate_pi'  0.03215 s \n",
       "\n",
       "t -69.27500915527344 26.84356689453125\n",
-      "'calculate_pi'  0.04249 s \n",
+      "'calculate_pi'  0.03131 s \n",
       "\n",
       "t -69.27500915527344 23.313018798828125\n",
-      "'calculate_pi'  0.04869 s \n",
+      "'calculate_pi'  0.02618 s \n",
       "\n",
       "t -69.27500915527344 15.280517578125\n",
-      "'calculate_pi'  0.03119 s \n",
+      "'calculate_pi'  0.02575 s \n",
       "\n",
       "t -69.27500915527344 31.399871826171875\n",
-      "'calculate_pi'  0.03782 s \n",
+      "'calculate_pi'  0.03177 s \n",
       "\n",
       "t -69.27500915527344 31.399871826171875\n",
-      "'calculate_pi'  0.03270 s \n",
+      "'calculate_pi'  0.03051 s \n",
       "\n"
      ]
     },

tcpips/pi.py

@@ -124,6 +124,9 @@ def fix_profile(
     2.  'lapse_rate': Fills NaNs at the bottom of the profile by extrapolating
         downwards from the lowest valid data point using a dry adiabatic lapse rate.
         This is ideal for fixing missing surface-level pressure data.
+    3.  'well_mixed': Fills NaNs at the bottom of the profile by assuming a well-mixed
+        boundary layer. This method fills both temperature and specific humidity
+        profiles consistently and is physically principled for tropical atmospheres.
 
     Args:
         ds (xr.Dataset): Xarray dataset containing temperature 't' with a

w22/img/new_orleans_timeseries.pdf

@@ -402,14 +402,43 @@ def place_to_position(place: str) -> tuple:
 
 def plot_spatials(
     axs: np.ndarray,
-    places: str = "new_orleans",
-    vars: Tuple[str, str] = ("vmax", "r0"),
+    place: str = "new_orleans",
+    vars: Tuple[str, str] = ("vmax_3", "r0_3", "rmax_3", "rmax_1"),
+    labels: tuple = (
+        "Potential intensity, $V_{p}$ [m s$^{-1}$]",
+        "Corresponding potential outer size, $r_{a3}$ [km]",
+        "Corresponding potential inner size, $r_{3}$ [km]",
+        "Cat1 potential inner size, $r_{1}$ [km]",
+    ),
     pi_version: int = 4,
     trial=1,
     pressure_assumption="isothermal",
 ) -> None:
     assert len(axs) == len(vars)
 
+    name = f"august_cmip6_pi{pi_version}_{pressure_assumption}_trial{trial}.nc"
+
+    if place in {"new_orleans", "galverston", "miami"}:
+        name = "gom_" + name
+    if place in {"shanghai", "hong_kong", "hanoi"}:
+        name = "scs_" + name
+    ds = xr.open_dataset(os.path.join(DATA_PATH, name))
+    print("spatial ds", ds)
+
+    ds["lon"].attrs = {"units": "$^{\circ}E$", "long_name": "Longitude"}
+    ds["lat"].attrs = {"units": "$^{\circ}N$", "long_name": "Latitude"}
+    for i, var in enumerate(vars):
+        if var in ds:
+            if var in {"r0", "rmax", "rmax_1", "rmax_3", "r0_1", "r0_3"}:
+                (ds[var] / 1000).plot(ax=axs[i], cbar_kwargs={"label": ""})
+            else:
+                ds[var].plot(ax=axs[i], cbar_kwargs={"label": ""})
+            axs[i].set_title(labels[i])
+        point = place_to_position(place)
+        axs[i].scatter(*point, color="black", s=75, marker="x")
+        if i != len(vars) - 1:
+            axs[i].set_xlabel("")
+
 
 def plot_two_spatial(
     axs: np.ndarray,
@@ -553,7 +582,13 @@ def get_cmip6_timeseries(
 
 def plot_timeserii(
     axs: np.ndarray,
-    serii: tuple = ("vmax_3", "r0_3", "rmax_3", "rmax_1"),
+    vars: tuple = ("vmax_3", "r0_3", "rmax_3", "rmax_1"),
+    labels: tuple = (
+        "Potential intensity, $V_{p}$ [m s$^{-1}$]",
+        "PI potential outer size, $r_{a3}$ [km]",
+        "PI potential inner size, $r_{3}$ [km]",
+        "Cat1 potential inner size, $r_{1}$ [km]",
+    ),
     color: str = "black",
     member: int = 4,
     model: str = "CESM2",
@@ -563,24 +598,34 @@ def plot_timeserii(
     year_min: int = 2014,
     year_max: int = 2100,
 ) -> None:
-    assert len(axs) == len(serii)
+    assert len(axs) == len(vars)
+    assert len(labels) == len(vars)
     ds = get_cmip6_timeseries(
         place=place,
         pressure_assumption=pressure_assumption,
         member=member,
         model=model,
         pi_version=pi_version,
     )
-    for i in range(len(serii)):
-        axs[i].set_title(serii[i])
+    for i in range(len(vars)):
+        axs[i].set_title(vars[i])
+        var_np = ds[vars[i]].values
+        if vars[i] in {"r0", "rmax", "rmax_1", "rmax_3", "r0_1", "r0_3"}:
+            var_np /= 1000
+
         axs[i].plot(
-            np.array([t.astype("datetime64[Y]").astype(int) + 1970 for t in ds.values]),
-            ds[serii[i]].values,
+            # np.array(
+            #     [t.astype("datetime64[Y]").astype(int) + 1970 for t in ds.time.values]
+            # ),
+            ds["time"].values,
+            var_np,
             color=color,
             linewidth=0.5,
+            alpha=0.5,
         )
+        axs[i].set_title(labels[i])
         axs[i].set_xlabel("")
-        if i == len(serii) - 1:
+        if i == len(vars) - 1:
             axs[i].set_xlabel("Year")
         axs[i].set_xlim([1850, 2100])
         # vertical black line at year_min
@@ -839,6 +884,73 @@ def figure_two(place: str = "new_orleans") -> None:
     plt.close()
 
 
+def multipanel(
+    place: str = "new_orleans",
+    vars: Tuple[str] = ("vmax_3", "r0_3", "rmax_3", "rmax_1"),
+):
+    """Plot the multipanel figure for the CMIP6 timeseries.
+
+    Args:
+        vars (tuple): The variables to plot (default is ("vmax_3", "r0_3", "rmax_3", "rmax_1")).
+    """
+    plot_defaults()
+    _, axs = plt.subplots(
+        len(vars),
+        2,
+        figsize=get_dim(ratio=1.5),
+    )
+    members = [4, 10, 11]
+    colors = ["orange", "green", "grey"]
+    plot_spatials(axs[:, 0], place="new_orleans", vars=vars)
+    for i, member in enumerate(members):
+        plot_timeserii(
+            axs[:, 1],
+            vars=vars,
+            color="orange",
+            member=member,
+            model="CESM2",
+            pressure_assumption="isothermal",
+            place="new_orleans",
+            pi_version=4,
+            year_min=2014,
+            year_max=2100,
+        )
+
+    members = ["r1i1p1f3", "r2i1p1f3", "r3i1p1f3"]
+    for i, member in enumerate(members):
+        plot_timeserii(
+            axs[:, 1],
+            vars=vars,
+            color="green",
+            member=member,
+            model="HADGEM3-GC31-MM",
+            pressure_assumption="isothermal",
+            place="new_orleans",
+            pi_version=4,
+            year_min=2014,
+            year_max=2100,
+        )
+    members = ["r1i1p1f1", "r2i1p1f1", "r3i1p1f1"]
+    for i, member in enumerate(members):
+        plot_timeserii(
+            axs[:, 1],
+            vars=vars,
+            color="grey",
+            member=member,
+            model="MIROC6",
+            pressure_assumption="isothermal",
+            place="new_orleans",
+            pi_version=4,
+            year_min=2014,
+            year_max=2100,
+        )
+
+    label_subplots(axs)
+    plt.savefig(os.path.join(FIGURE_PATH, f"{place}_multipanel.pdf"))
+    plt.clf()
+    plt.close()
+
+
 def temporal_relationship_data(place: str = "new_orleans", pi_version: int = 4) -> None:
     """Get the temporal relationships data for the given place and potential intensity version.
 
@@ -1131,8 +1243,8 @@ def _generate_header_map(columns: list[str]) -> dict[str, str]:
     # python -m w22.plot
     # plot_panels()
     #
-    figure_two(place="new_orleans")
-    figure_two(place="hong_kong")
+    # figure_two(place="new_orleans")
+    # figure_two(place="hong_kong")
     # temporal_relationship_data(place="new_orleans", pi_version=4)
     # plot_seasonal_profiles()
     # years = [2015, 2100]
@@ -1149,3 +1261,5 @@ def _generate_header_map(columns: list[str]) -> dict[str, str]:
     #     years=years,
     #     members=[4, 10],
     # )
+    multipanel(place="new_orleans")
+    multipanel(place="hong_kong")