Front Lifecycle Analysis¶

# -----------------------------------
# LIBRARY IMPORTS
# -----------------------------------
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from descartes import PolygonPatch

import cartopy
from matplotlib.colors import ListedColormap

from fronttracker import FrontTracker

# -----------------------------------
# LOAD DATA
# -----------------------------------
gm = np.transpose(np.load(".../gm.npy"), (2, 1, 0))
latitude = np.load(".../latitude.npy")
longitude = np.load(".../longitude.npy")
time = np.load(".../time.npy", allow_pickle=True)
time = np.arange(np.datetime64(time[0].strftime("%Y-%m-%d")), np.datetime64(time[-1].strftime("%Y-%m-%d")) + np.timedelta64(1, 'D'),np.timedelta64(1, 'D'))

# -----------------------------------
# FRONT THRESHOLD
# -----------------------------------
treshold = 0.08

# -----------------------------------
# FRONT IDENTIFICATION AND TRACKING
# -----------------------------------

# Initialize the front tracker
front_tracker = FrontTracker()

# Identify fronts in the GM field using the threshold defined earlier
front_tracker.identify_fronts(time, longitude, latitude, gm, treshold=treshold)

# Retrieve the tabular data of all detected front pixels
front_tracker.data

# Link fronts between consecutive time steps (temporal tracking)
front_tracker.track()

# Retrieve a DataFrame that maps each front label to its next label in the tracking path
front_tracker.data_labels

# Build the frontal path by combining automated detection with manual corrections
frontal_path = [None]
frontal_path += front_tracker.get_frontal_path(7)
frontal_path += front_tracker.get_frontal_path(548)[:2]
frontal_path += front_tracker.get_frontal_path(559)
frontal_path += front_tracker.get_frontal_path(579)[:3]
frontal_path += front_tracker.get_frontal_path(598)
frontal_path += front_tracker.get_frontal_path(700)
frontal_path += front_tracker.get_frontal_path(728)
frontal_path += front_tracker.get_frontal_path(789)
frontal_path.append(None)

print(f"Frontal path length: {len(frontal_path)}")
print(frontal_path)

# -----------------------------------
# VISUALIZATION SETUP
# -----------------------------------

# Create a large figure with a grid of 11 rows × 14 columns
fig = plt.figure(figsize=(4.8*4, 6.2*4))
nrows, ncols = 11, 14

for k in range(len(time)):
    ax = fig.add_subplot(nrows, ncols, k+1, projection=cartopy.crs.PlateCarree())
    
    land = cartopy.feature.NaturalEarthFeature('physical', 'land', scale="10m", edgecolor='k', facecolor=cartopy.feature.COLORS['land'])
    ax.add_feature(land, facecolor='lightgray', linewidth=0.25)
    gl = ax.gridlines(draw_labels=True, linewidth=0.25, color='gray', alpha=0.5, linestyle='--')
    gl.right_labels = False
    if k % ncols != 0:
        gl.left_labels = False 
    
    # Plot gradient magnitude (GM) field
    pc = ax.pcolor(longitude, latitude, gm[k], vmin=0.0, vmax=0.3, cmap="jet")
    
    # If a frontal path exists for this timestep
    if frontal_path[k] is not None:
        label = frontal_path[k]
        
        # Get all detected fronts at this time
        _df = front_tracker.data[front_tracker.data.time == time[k]]
        uniq = _df.labels.unique()
    
        for _label in uniq:
            front = front_tracker.fronts[_label]
            
            # Try to compute contours
            try:
                front.contouring()
            except Exception:
                continue
                
            # Highlight the tracked front in pink, others in orange
            color = "deeppink" if _label == label else "orange"

            # Draw each polygon of the front contours
            for poly in front.contours:
                x, y = poly.exterior.xy
                ax.plot(x, y, c=color, linewidth=2)
                
    ax.text(102, 4.8, k, color="red", va="center", ha="center", fontsize=16)

plt.show()

# -----------------------------------
# PLOT FRONT METRICS
# -----------------------------------

# Create a large figure with a grid of 2 rows × 1 column
fig, axs = plt.subplots(2, 1, figsize=(6.2*2, 4.8*0.75))

# --- Panel 1: Frontal area ---
axs[0].plot(np.arange(0, len(time)), metrics["area"])
axs[0].set_ylim(0,65000)
axs[0].set_yticks(np.arange(0, 60001, 15000))
axs[0].set_ylabel(r"Area (km$^2$)")
axs[0].set_xlim(0,len(time)-1)
axs[0].set_xticks(np.arange(0, 150, 5))

# --- Panel 2: Gradient magnitude statistics ---
axs[1].plot(time[1:-1], metrics["gm_mean"][1:-1])
axs[1].fill_between(time[1:-1], metrics["gm_min"][1:-1], metrics["gm_max"][1:-1], alpha=0.3)
axs[1].set_ylim(0.05,0.3)
axs[1].set_yticks(np.arange(0.05, 0.31, 0.05))
axs[1].set_ylabel(r"GM (°C (10 km)$^{-1}$)")
axs[1].set_xlim(time[0],time[-1])
axs[1].set_xticks(time[::5])
axs[1].tick_params(axis='x', rotation=90)

plt.show()