|
| 1 | +""" pyplots.ai |
| 2 | +skewt-logp-atmospheric: Skew-T Log-P Atmospheric Diagram |
| 3 | +Library: letsplot 4.8.2 | Python 3.13.11 |
| 4 | +Quality: 90/100 | Created: 2026-01-17 |
| 5 | +""" |
| 6 | + |
| 7 | +import numpy as np |
| 8 | +import pandas as pd |
| 9 | +from lets_plot import ( |
| 10 | + LetsPlot, |
| 11 | + aes, |
| 12 | + element_blank, |
| 13 | + element_text, |
| 14 | + geom_path, |
| 15 | + geom_point, |
| 16 | + geom_segment, |
| 17 | + geom_text, |
| 18 | + ggplot, |
| 19 | + ggsave, |
| 20 | + ggsize, |
| 21 | + labs, |
| 22 | + scale_y_log10, |
| 23 | + scale_y_reverse, |
| 24 | + theme, |
| 25 | +) |
| 26 | + |
| 27 | + |
| 28 | +LetsPlot.setup_html() |
| 29 | + |
| 30 | +# Atmospheric sounding data (synthetic radiosonde profile) |
| 31 | +np.random.seed(42) |
| 32 | + |
| 33 | +# Pressure levels from surface to upper troposphere (hPa) |
| 34 | +pressure = np.array([1000, 950, 900, 850, 800, 750, 700, 650, 600, 550, 500, 450, 400, 350, 300, 250, 200, 150, 100]) |
| 35 | + |
| 36 | +# Temperature profile (°C) - typical mid-latitude sounding |
| 37 | +# Decreasing with height, with a tropopause around 200 hPa |
| 38 | +temperature = np.array([25, 22, 18, 14, 10, 6, 2, -3, -8, -14, -21, -29, -38, -47, -55, -58, -56, -55, -56]) |
| 39 | + |
| 40 | +# Dewpoint profile (°C) - typically lower than temperature |
| 41 | +# Shows moisture decreasing with height |
| 42 | +dewpoint = np.array([18, 16, 12, 8, 4, 0, -5, -12, -20, -28, -35, -42, -50, -58, -65, -70, -72, -75, -78]) |
| 43 | + |
| 44 | +# For Skew-T: we need to skew the temperature by adding an offset based on pressure |
| 45 | +# The skew factor converts the diagram coordinates |
| 46 | +# In a Skew-T, x_plot = T + skew_factor * log(p0/p) |
| 47 | +p0 = 1000 # Reference pressure |
| 48 | +skew_factor = 40 # Degrees of skew per pressure decade |
| 49 | + |
| 50 | +# Calculate skewed coordinates for plotting |
| 51 | +log_pressure = np.log10(p0 / pressure) |
| 52 | +temp_skewed = temperature + skew_factor * log_pressure |
| 53 | +dewpoint_skewed = dewpoint + skew_factor * log_pressure |
| 54 | + |
| 55 | +# Create DataFrame for main profiles with type for legend |
| 56 | +df_temp = pd.DataFrame({"pressure": pressure, "value": temp_skewed, "type": "Temperature"}) |
| 57 | + |
| 58 | +df_dewpoint = pd.DataFrame({"pressure": pressure, "value": dewpoint_skewed, "type": "Dewpoint"}) |
| 59 | + |
| 60 | +# Generate isotherms (temperature reference lines, skewed 45 degrees) |
| 61 | +isotherm_temps = np.arange(-80, 50, 10) # °C |
| 62 | +isotherm_data = [] |
| 63 | +p_range = np.array([1000, 100]) |
| 64 | +for t in isotherm_temps: |
| 65 | + log_p_vals = np.log10(p0 / p_range) |
| 66 | + t_skewed = t + skew_factor * log_p_vals |
| 67 | + isotherm_data.append( |
| 68 | + {"x_start": t_skewed[0], "x_end": t_skewed[1], "y_start": p_range[0], "y_end": p_range[1], "temp": t} |
| 69 | + ) |
| 70 | +df_isotherms = pd.DataFrame(isotherm_data) |
| 71 | + |
| 72 | +# Generate dry adiabats (lines of constant potential temperature) |
| 73 | +# theta = T * (p0/p)^(R/cp) where R/cp ≈ 0.286 |
| 74 | +dry_adiabat_thetas = np.arange(250, 450, 20) # K |
| 75 | +dry_adiabat_data = [] |
| 76 | +p_levels = np.linspace(1000, 100, 50) |
| 77 | +for theta in dry_adiabat_thetas: |
| 78 | + temps = (theta * (p_levels / p0) ** 0.286) - 273.15 # Convert to Celsius |
| 79 | + log_p_vals = np.log10(p0 / p_levels) |
| 80 | + temps_skewed = temps + skew_factor * log_p_vals |
| 81 | + for i in range(len(p_levels)): |
| 82 | + dry_adiabat_data.append( |
| 83 | + {"pressure": p_levels[i], "temp_skewed": temps_skewed[i], "theta": theta, "line_type": "Dry Adiabat"} |
| 84 | + ) |
| 85 | +df_dry_adiabats = pd.DataFrame(dry_adiabat_data) |
| 86 | + |
| 87 | +# Generate moist adiabats (equivalent potential temperature lines) |
| 88 | +# Moist adiabats follow a different curve - steeper at low levels, flatter aloft |
| 89 | +# theta_e is conserved along moist adiabats |
| 90 | +moist_adiabat_theta_e = np.arange(280, 360, 10) # K |
| 91 | +moist_adiabat_data = [] |
| 92 | +for theta_e in moist_adiabat_theta_e: |
| 93 | + # Approximate moist adiabat using simplified formula |
| 94 | + # At surface, T ≈ theta_e - 273.15, then follows moist lapse rate (~6.5°C/km) |
| 95 | + t_surface = theta_e - 273.15 |
| 96 | + for p in p_levels: |
| 97 | + # Moist adiabatic lapse rate is variable, roughly 6.5 C/km |
| 98 | + # Use pseudo-adiabatic approximation |
| 99 | + height_factor = np.log(p0 / p) * 2.5 # Scale factor for height |
| 100 | + t_moist = t_surface - 6.5 * height_factor * 0.8 # Moist lapse rate effect |
| 101 | + log_p_val = np.log10(p0 / p) |
| 102 | + t_skewed = t_moist + skew_factor * log_p_val |
| 103 | + moist_adiabat_data.append( |
| 104 | + {"pressure": p, "temp_skewed": t_skewed, "theta_e": theta_e, "line_type": "Moist Adiabat"} |
| 105 | + ) |
| 106 | +df_moist_adiabats = pd.DataFrame(moist_adiabat_data) |
| 107 | + |
| 108 | +# Generate mixing ratio lines (lines of constant water vapor mixing ratio) |
| 109 | +# Simplified: ws = 622 * es(T) / (p - es(T)) where es is saturation vapor pressure |
| 110 | +mixing_ratios = [1, 2, 4, 7, 10, 15, 20] # g/kg |
| 111 | +mixing_data = [] |
| 112 | +for ws in mixing_ratios: |
| 113 | + for p in p_levels[::5]: # Sample every 5th pressure level for smooth lines |
| 114 | + # Approximate temperature for given mixing ratio and pressure |
| 115 | + # Using simplified formula: T ≈ 35 * log(ws * p / 622) - 20 |
| 116 | + t = 35 * np.log10(ws * p / 622) - 20 |
| 117 | + log_p_val = np.log10(p0 / p) |
| 118 | + t_skewed = t + skew_factor * log_p_val |
| 119 | + mixing_data.append({"pressure": p, "temp_skewed": t_skewed, "ws": ws, "line_type": "Mixing Ratio"}) |
| 120 | +df_mixing = pd.DataFrame(mixing_data) |
| 121 | + |
| 122 | +# Build the plot using lets-plot |
| 123 | +plot = ( |
| 124 | + ggplot() |
| 125 | + # Isotherms (skewed temperature lines) - light gray |
| 126 | + + geom_segment( |
| 127 | + aes(x="x_start", xend="x_end", y="y_start", yend="y_end"), |
| 128 | + data=df_isotherms, |
| 129 | + color="#999999", |
| 130 | + size=0.6, |
| 131 | + alpha=0.8, |
| 132 | + ) |
| 133 | + # Dry adiabats - orange dashed (increased visibility) |
| 134 | + + geom_path( |
| 135 | + aes(x="temp_skewed", y="pressure", group="theta"), |
| 136 | + data=df_dry_adiabats, |
| 137 | + color="#D97706", |
| 138 | + size=0.8, |
| 139 | + alpha=0.7, |
| 140 | + linetype="dashed", |
| 141 | + ) |
| 142 | + # Moist adiabats - purple dash-dot |
| 143 | + + geom_path( |
| 144 | + aes(x="temp_skewed", y="pressure", group="theta_e"), |
| 145 | + data=df_moist_adiabats, |
| 146 | + color="#7C3AED", |
| 147 | + size=0.8, |
| 148 | + alpha=0.7, |
| 149 | + linetype="dotdash", |
| 150 | + ) |
| 151 | + # Mixing ratio lines - green dotted (increased visibility) |
| 152 | + + geom_path( |
| 153 | + aes(x="temp_skewed", y="pressure", group="ws"), |
| 154 | + data=df_mixing, |
| 155 | + color="#059669", |
| 156 | + size=0.8, |
| 157 | + alpha=0.7, |
| 158 | + linetype="dotted", |
| 159 | + ) |
| 160 | + # Temperature profile - solid red line |
| 161 | + + geom_path(aes(x="value", y="pressure"), data=df_temp, color="#DC2626", size=2.5) |
| 162 | + # Dewpoint profile - dashed blue line |
| 163 | + + geom_path(aes(x="value", y="pressure"), data=df_dewpoint, color="#2563EB", size=2.5, linetype="dashed") |
| 164 | + # Points on profiles for clarity |
| 165 | + + geom_point(aes(x="value", y="pressure"), data=df_temp, color="#DC2626", size=4) |
| 166 | + + geom_point(aes(x="value", y="pressure"), data=df_dewpoint, color="#2563EB", size=4) |
| 167 | + # Reverse y-axis (pressure decreases upward) and log scale |
| 168 | + + scale_y_log10() |
| 169 | + + scale_y_reverse(limits=[1000, 100]) |
| 170 | + # Labels and title |
| 171 | + + labs(x="Temperature (°C) - Skewed", y="Pressure (hPa)", title="skewt-logp-atmospheric · letsplot · pyplots.ai") |
| 172 | + + ggsize(1600, 900) |
| 173 | + + theme( |
| 174 | + axis_title=element_text(size=20), |
| 175 | + axis_text=element_text(size=16), |
| 176 | + plot_title=element_text(size=24), |
| 177 | + panel_grid_major=element_blank(), |
| 178 | + panel_grid_minor=element_blank(), |
| 179 | + panel_background=element_blank(), |
| 180 | + axis_line=element_blank(), |
| 181 | + ) |
| 182 | +) |
| 183 | + |
| 184 | +# Create legend using actual line segments with matching styles |
| 185 | +# Position legend in upper right area |
| 186 | +legend_x_start = 95 |
| 187 | +legend_x_end = 115 |
| 188 | +legend_text_x = 118 |
| 189 | + |
| 190 | +# Legend items with their properties matching the actual plot lines |
| 191 | +legend_entries = [ |
| 192 | + {"label": "Temp", "color": "#DC2626", "linetype": "solid", "size": 2.5, "y": 125}, |
| 193 | + {"label": "Dewpt", "color": "#2563EB", "linetype": "dashed", "size": 2.5, "y": 140}, |
| 194 | + {"label": "Dry Ad.", "color": "#D97706", "linetype": "dashed", "size": 0.8, "y": 158}, |
| 195 | + {"label": "Moist Ad.", "color": "#7C3AED", "linetype": "dotdash", "size": 0.8, "y": 178}, |
| 196 | + {"label": "Mix. Ratio", "color": "#059669", "linetype": "dotted", "size": 0.8, "y": 200}, |
| 197 | +] |
| 198 | + |
| 199 | +# Add legend line segments with actual linetypes |
| 200 | +for entry in legend_entries: |
| 201 | + plot = plot + geom_segment( |
| 202 | + aes(x="x_start", xend="x_end", y="y", yend="y"), |
| 203 | + data=pd.DataFrame([{"x_start": legend_x_start, "x_end": legend_x_end, "y": entry["y"]}]), |
| 204 | + color=entry["color"], |
| 205 | + size=entry["size"], |
| 206 | + linetype=entry["linetype"], |
| 207 | + ) |
| 208 | + plot = plot + geom_text(x=legend_text_x, y=entry["y"], label=entry["label"], color="#333333", size=12, hjust=0) |
| 209 | + |
| 210 | +# Save the plot |
| 211 | +ggsave(plot, "plot.png", path=".", scale=3) |
| 212 | + |
| 213 | +# Save HTML for interactive version |
| 214 | +ggsave(plot, "plot.html", path=".") |
0 commit comments