ATM 509 Final Project Figures - Evan Belkin¶
Imports¶
import pandas as pd
import xarray as xr
import numpy as np
import matplotlib.pyplot as plt
from metpy.plots import USCOUNTIES
import cartopy.crs as ccrs
from scipy.interpolate import griddata
import cartopy.feature as cfeature
from matplotlib.colors import ListedColormap, BoundaryNorm
import pickle
import cmweather
import matplotlib.dates as mdates
from matplotlib.dates import DateFormatter, AutoDateLocator,HourLocator,DayLocator
from matplotlib.ticker import FuncFormatter
from herbie import Herbie
from matplotlib.colors import Normalize
from datetime import datetime, timedeltaLoading in Data:¶
monthly_percentiles = pd.read_pickle("../data/509_Cases.pkl")
Figure #1: Extract 90th percentile column and create a 4x3 plot.¶
def plot_lat_lon_12panel(dataframe):
# Define month names
month_names = ["January", "February", "March", "April", "May", "June",
"July", "August", "September", "October", "November", "December"]
# Ensure month column is integer for sorting
dataframe['Month'] = dataframe['Month'].astype(int)
months = sorted(dataframe['Month'].unique())
# Set up the figure with 12 subplots (4 rows x 3 columns)
fig, axes = plt.subplots(nrows=4, ncols=3, figsize=(20, 14),
subplot_kw={'projection': ccrs.PlateCarree()},
dpi=150)
axes = axes.flatten()
# Define fixed color bins
bounds = np.arange(5, 50, 5)
cmap = plt.get_cmap('NWSRef', len(bounds)-1)
norm = BoundaryNorm(bounds, cmap.N)
for i, month in enumerate(months):
ax = axes[i]
month_data = dataframe[dataframe['Month'] == month]
latitudes = month_data['center_lat'].values
longitudes = month_data['center_lon'].values
variable = month_data['90'].values ## Extract 90th percentile!
ax.set_extent([-81, -67, 38.8, 48], crs=ccrs.PlateCarree())
# Map features
ax.add_feature(cfeature.STATES, linestyle='--', linewidth=3)
ax.add_feature(cfeature.BORDERS, linestyle=':', linewidth=0.5)
ax.add_feature(cfeature.LAND, edgecolor='gray')
ax.add_feature(USCOUNTIES.with_scale('20m'), edgecolor='gray', linewidth=0.4)
# Scatter plot with fixed color bins
sc = ax.scatter(longitudes, latitudes, c=variable, s=2,
cmap=cmap, norm=norm, transform=ccrs.PlateCarree())
# Month label
ax.text(-76.5, 47.0, month_names[month-1],
fontsize=24, fontweight='bold',
transform=ccrs.PlateCarree(),
bbox=dict(facecolor='white', edgecolor='none', alpha=0.7))
# Add a shared colorbar
cbar_ax = fig.add_axes([0.92, 0.1, 0.015, 0.8])
cbar = fig.colorbar(sc, cax=cbar_ax, extend='both', ticks=bounds)
cbar.set_label('90th Percentile (mm)', fontsize=24)
cbar.ax.tick_params(labelsize=24)
# Layout
plt.suptitle("Monthly 90th Percentile of Measurable Precipitation", fontsize=36)
plt.tight_layout(rect=[0, 0, 0.9, 1], pad=1.5, w_pad=0.3, h_pad=0.3)
#plt.savefig('evan.png')
plt.show()
plot_lat_lon_12panel(monthly_percentiles)
Figure #2: JFK Airport Example Event¶
time_series = pd.read_pickle("../data/509_JFK.pkl")
def custom_hour_formatter(x, pos):
dt = mdates.num2date(x)
if dt.hour == 0:
return dt.strftime('%H %d %b')
else:
return dt.strftime('%H')
rolling_24h = time_series['24h_rolling_sum']
time = time_series['Datetime']
precipitation = time_series['Precipitation']
event_rolling24h = time_series['event_rolling_24h']
fig, ax = plt.subplots(figsize=(14, 6))
ax.axhline(y=25.2380, c='red', alpha=0.7, linewidth=2, label='June 90th percentile')
ax.plot(time, rolling_24h, marker='o', linestyle='-', label='24 h rolling precip', color='b')
ax.plot(time, precipitation, marker='o', linestyle='-', label='Hourly precip', color='g')
ax.set_xlabel('Time (UTC)', fontsize=15)
ax.set_ylabel('Precip (mm)', fontsize=15)
ax.xaxis.set_major_locator(mdates.HourLocator(byhour=[0, 6, 12, 18]))
ax.xaxis.set_minor_locator(mdates.HourLocator(interval=1))
ax.xaxis.set_major_formatter(FuncFormatter(custom_hour_formatter))
specific_time = pd.Timestamp('2018-04-16 00:00:00')
ax.axvline(x=specific_time, color='purple', linewidth=5.0)
specific_time = pd.Timestamp('2018-04-16 05:00:00')
ax.axvline(x=specific_time, color='orange', linewidth=3.5, zorder=1)
specific_time = pd.Timestamp('2018-04-16 17:00:00')
ax.axvline(x=specific_time, color='purple', linewidth=5.0)
plt.ylim(-1,80)
ax.tick_params(axis='x', pad=8, size=10)
plt.xticks(fontsize=15)
plt.yticks(fontsize=15)
plt.grid()
Figure #3A: Topography. Here is the code! Elev data is data from Herbie! :)¶
H = Herbie(
date= datetime(2020, 12, 24, 6, 0, 0), # initialization date and time: year, month, day, hour, miunte, second
model="hrrr",
product="sfc",
fxx=6,
)
elev = H.xarray(':HGT:surface').orogfig = plt.figure(figsize=(10, 10), layout='constrained')
ax = fig.add_subplot(111, projection=ccrs.PlateCarree())
ax.set_title('Topography of the Northeast United States', fontsize=16)
plot_bounds = [-81, -67, 38.8, 48]
ax.set_extent(plot_bounds)
# Define levels using full 0–1200 scale
levels = np.arange(100, 1201, 50)
# Mask values < 60 so they are not plotted
masked_elev = elev.where(elev >= 100)
elev_cb = ax.contourf(
elev.longitude, elev.latitude, masked_elev,
levels=levels, cmap='gist_earth', extend='max'
)
cbar = plt.colorbar(
elev_cb, shrink=0.6,
ticks=np.arange(100, 1201, 300),
extend='max'
)
cbar.ax.tick_params(labelsize=12)
cbar.set_label('Elevation (m)', fontsize=12)
ax.add_feature(cfeature.STATES, linestyle='--', linewidth=3)
# Add county boundaries
# Add red stars and labels
ax.plot(-73.7562, 42.6526, marker='*', color='red', markersize=15,
transform=ccrs.PlateCarree(), label='Albany, NY')
ax.plot(-75.85, 43.65, marker='*', color='red', markersize=15,
transform=ccrs.PlateCarree(), label='Tug Hill Plateau')
ax.text(-73.5562, 42.4526, 'Albany, NY', transform=ccrs.PlateCarree(),
fontsize=14, color='red', weight='bold')
ax.text(-75.65, 43.45, 'Tug Hill Plateau', transform=ccrs.PlateCarree(),
fontsize=14, color='red', weight='bold')
plt.show()
Figure #3B: Total # of Cases.¶
latitudes = monthly_percentiles['center_lat'].values
longitudes = monthly_percentiles['center_lon'].values
cases = monthly_percentiles['Cases'].values
# Define bounds for discrete colorbar
bounds = np.arange(250, 525, 25) # 250–500 in steps of 25
cmap = plt.get_cmap('NWSRef', len(bounds) - 1)
norm = plt.matplotlib.colors.BoundaryNorm(bounds, cmap.N)
# Plot
fig, ax = plt.subplots(figsize=(10, 8), subplot_kw={'projection': ccrs.PlateCarree()})
ax.set_extent([-81, -67, 38.8, 48], crs=ccrs.PlateCarree())
ax.add_feature(cfeature.STATES, linestyle='--')
ax.add_feature(cfeature.BORDERS, linestyle=':')
ax.add_feature(cfeature.LAND, edgecolor='gray')
# Counties (optional, if USCOUNTIES is available)
ax.add_feature(USCOUNTIES.with_scale('20m'), edgecolor='gray')
# Scatter plot
sc = ax.scatter(longitudes, latitudes, c=cases, s=1,
cmap=cmap, norm=norm, transform=ccrs.PlateCarree())
# Colorbar with fixed bounds
cbar = plt.colorbar(sc, ax=ax, boundaries=bounds, ticks=bounds,
shrink=0.65, extend='max')
cbar.set_label("Cases")
plt.title("Total Number of Cases", fontsize=14)
plt.show()
Figure #4: Event Duration by Month¶
duration = pd.read_pickle("../data/509_duration.pkl")
durationdef plot_lat_lon_12panel(dataframe):
# Define month names
month_names = ["January", "February", "March", "April", "May", "June",
"July", "August", "September", "October", "November", "December"]
# Ensure month column is integer for sorting
dataframe['month'] = dataframe['month'].astype(int)
months = sorted(dataframe['month'].unique())
# Set up the figure with 12 subplots (4 rows x 3 columns)
fig, axes = plt.subplots(nrows=4, ncols=3, figsize=(20, 14),
subplot_kw={'projection': ccrs.PlateCarree()},
dpi=150)
axes = axes.flatten()
# Define fixed color bins
bounds = np.arange(8, 40, 4)
cmap = plt.get_cmap('NWSRef', len(bounds)-1)
norm = BoundaryNorm(bounds, cmap.N)
for i, month in enumerate(months):
ax = axes[i]
month_data = dataframe[dataframe['month'] == month]
latitudes = month_data['center_lat'].values
longitudes = month_data['center_lon'].values
variable = month_data['duration'].values ## Extract 90th percentile!
ax.set_extent([-81, -67, 38.8, 48], crs=ccrs.PlateCarree())
# Map features
ax.add_feature(cfeature.STATES, linestyle='--', linewidth=3)
ax.add_feature(cfeature.BORDERS, linestyle=':', linewidth=0.5)
ax.add_feature(cfeature.LAND, edgecolor='gray')
ax.add_feature(USCOUNTIES.with_scale('20m'), edgecolor='gray', linewidth=0.4)
# Scatter plot with fixed color bins
sc = ax.scatter(longitudes, latitudes, c=variable, s=2,
cmap=cmap, norm=norm, transform=ccrs.PlateCarree())
# Month label
ax.text(-76.5, 47.0, month_names[month-1],
fontsize=24, fontweight='bold',
transform=ccrs.PlateCarree(),
bbox=dict(facecolor='white', edgecolor='none', alpha=0.7))
# Add a shared colorbar
cbar_ax = fig.add_axes([0.92, 0.1, 0.015, 0.8])
cbar = fig.colorbar(sc, cax=cbar_ax, extend='max', ticks=bounds)
cbar.set_label('Monthly Mean Event Duration (h)', fontsize=24)
cbar.ax.tick_params(labelsize=24)
# Layout
plt.suptitle("Monthly Mean Event Duration", fontsize=36)
plt.tight_layout(rect=[0, 0, 0.9, 1], pad=1.5, w_pad=0.3, h_pad=0.3)
#plt.savefig('evan.png')
plt.show()plot_lat_lon_12panel(duration)
Figure 5: Event Duration Plots¶
fall = pd.read_pickle("../data/509_fall.pkl")
winter = pd.read_pickle("../data/509_winter.pkl")
spring = pd.read_pickle("../data/509_spring.pkl")
summer = pd.read_pickle("../data/509_summer.pkl")
fig, axes = plt.subplots(2, 2, figsize=(14, 10))
# Font sizes
title_fontsize = 18
label_fontsize = 16
tick_label_size = 14
supertitle_fontsize = 20
# Define a helper function to plot each season
def plot_season(ax, data, title):
# Plot non-AR events in default color
non_ar = data[data['AR_association'] == False]
ax.scatter(non_ar['duration'], non_ar['avgrate_1h'], label='Non-AR', alpha=0.7)
# Plot AR-associated events in red
ar = data[data['AR_association'] == True]
ax.scatter(ar['duration'], ar['avgrate_1h'], color='red', label='AR-associated', alpha=0.7)
ax.set_xlabel('Event Duration', fontsize=label_fontsize)
ax.set_ylabel('Event Average 1 h Precip Rate', fontsize=label_fontsize)
ax.set_title(title, fontsize=title_fontsize)
ax.set_ylim(0, 20)
ax.set_xlim(0, 70)
ax.grid(True)
ax.legend(fontsize=12)
ax.tick_params(axis='both', labelsize=tick_label_size)
# Plot each season
plot_season(axes[0, 0], fall, 'FALL')
plot_season(axes[0, 1], winter, 'WINTER')
plot_season(axes[1, 0], spring, 'SPRING')
plot_season(axes[1, 1], summer, 'SUMMER')
# Set the main title
fig.suptitle('Event Duration vs. Average 1 h Precipitation Rate by Season – Albany, NY',
fontsize=supertitle_fontsize)
# Adjust layout
plt.tight_layout(rect=[0, 0, 1, 0.95])
# Show plot
plt.show()
Figure 6: Seasonal Breakdown of Events:¶
seasonal_stats = pd.read_pickle("../data/509_seasonal_stats.pkl")
# Calculate non-AR event counts
non_ar_count = seasonal_stats['Total_Count'] - seasonal_stats['True_Count']
# Plot
plt.figure(figsize=(10, 6))
plt.bar(seasonal_stats.index, seasonal_stats['True_Count'], label='AR-Associated Events', color='orange')
plt.bar(seasonal_stats.index, non_ar_count, bottom=seasonal_stats['True_Count'], label='Non-AR Events', color='steelblue')
# Labels and formatting
plt.xlabel('Season', fontsize=14)
plt.ylabel('Number of Events', fontsize=14)
plt.title(f'Seasonal Breakdown of Events – Albany, NY', fontsize=16)
plt.xticks(fontsize=12)
plt.yticks(fontsize=12)
plt.legend(fontsize=12)
plt.grid(axis='y')
plt.ylim(0,150)
plt.tight_layout()
plt.show()
Figure 7: AR Event Footprint¶
percent_AR_seasonal_df = pd.read_pickle("../data/509_percents.pkl")point_lon, point_lat = -73.819, 42.735 ###ALBANY
seasons = ['DJF', 'MAM', 'JJA', 'SON']
norm = Normalize(vmin=10, vmax=80)
fig, axs = plt.subplots(2, 2, figsize=(16, 12), subplot_kw={'projection': ccrs.PlateCarree()})
axs = axs.flatten()
# Loop over each season and plot
for i, season in enumerate(seasons):
ax = axs[i]
season_df = percent_AR_seasonal_df[percent_AR_seasonal_df['season'] == season]
scatter = ax.scatter(
season_df['lon'],
season_df['lat'],
c=season_df['percent_AR'],
cmap='viridis',
norm=norm,
s=35,
transform=ccrs.PlateCarree()
)
ax.plot(
point_lon, point_lat,
marker='*',
color='red',
alpha=0.5,
markersize=12,
markeredgecolor='black',
transform=ccrs.PlateCarree(),
zorder=5
)
#ax.set_title(season, fontsize=18)
ax.text(-75.5, 46.2, season,
fontsize=24, fontweight='bold',
transform=ccrs.PlateCarree(),
bbox=dict(facecolor='white', edgecolor='none', alpha=0.7))
ax.add_feature(cfeature.STATES, linestyle='--', edgecolor='black', linewidth=0.8)
ax.add_feature(cfeature.BORDERS, linestyle=':', edgecolor='black', linewidth=0.8)
ax.add_feature(cfeature.LAND, edgecolor='gray')
ax.add_feature(USCOUNTIES.with_scale('20m'), edgecolor='gray', linewidth=0.5)
cbar_ax = fig.add_axes([0.92, 0.15, 0.03, 0.8]) # [left, bottom, width, height]
cbar = fig.colorbar(scatter, cax=cbar_ax, extend='both', label='% AR Events')
cbar.set_label('% AR Events', fontsize=18)
cbar.ax.tick_params(labelsize=18) # Set tick labels' font size
#fig.suptitle('Seasonal Percent of Events Associated with Atmospheric Rivers', fontsize=18)
#fig.suptitle('Seasonal AR Frequency During Albany Extreme Precip Events', fontsize=18)
fig.subplots_adjust(wspace=0.01, hspace=0.00001, right=0.9, top=0.92)
plt.show()