Road Icing Index Python Code

Appendix B: Road Icing Index Python Code.

##################################################################################

# Hybrid Precipitation Type Calculation Function (Bourgouin, Ramer, Baldwin, etc.)

##################################################################################

def calculate_hybrid_precipitation_type(temp_3d, pressure_3d, t2_da, psfc_da, qrain, qsnow, qgraup, wet_bulb_temp):

"""

Determine precipitation type using a hybrid method combining key aspects of the

Bourgouin, Ramer, and Baldwin algorithms, plus an ensemble approach for final classification.

This function examines vertical temperature profiles up to ~700 mb, layer thicknesses,

and surface mixing ratios (rain, snow, graupel) to classify precipitation types

such as Rain, Freezing Rain, Snow, Sleet/Graupel, or None.

Parameters

----------

temp_3d : xarray.DataArray

3D temperature in Kelvin (perturbation potential temperature already converted to real T).

pressure_3d : xarray.DataArray

3D pressure in Pa (perturbation + base state).

t2_da : xarray.DataArray

2 m temperature (K).

psfc_da : xarray.DataArray

Surface pressure (Pa).

qrain : xarray.DataArray

3D rain mixing ratio (kg/kg).

qsnow : xarray.DataArray

3D snow mixing ratio (kg/kg).

qgraup : xarray.DataArray

3D graupel mixing ratio (kg/kg).

wet_bulb_temp : ndarray

2D wet-bulb temperature at the surface (K).

Returns

-------

ptype : ndarray (2D)

Precipitation type classification array (strings):

['Rain', 'Freezing Rain', 'Snow', 'Sleet/Graupel', 'None']

Notes

-----

- Truncates the temperature profile near ~700 mb instead of 850 mb, capturing deeper warm layers

that could influence freezing rain/sleet formation.

- Uses ensemble-like voting from multiple methods (Bourgouin, Ramer, Baldwin).

- Resolves ambiguous classification with wet-bulb temperature checks and mixing ratios at surface.

"""

# Convert near-surface temperature to Celsius

t2_c = t2_da.values - 273.15

wb_c = wet_bulb_temp - 273.15 # 2D array of wet-bulb in °C

# Truncate temperature profile to ~700 mb

ref_pressure_profile = pressure_3d.isel(south_north=0, west_east=0).values

idx_700 = np.argmin(np.abs(ref_pressure_profile - 700.0)) # find 700 mb

temp_profile_700 = temp_3d.isel(bottom_top=slice(0, idx_700 + 1))

# Extract surface-level mixing ratios

qrain_surf = qrain.isel(bottom_top=0).values

qsnow_surf = qsnow.isel(bottom_top=0).values

qgraup_surf = qgraup.isel(bottom_top=0).values

# Convert truncated profile to Celsius

temp_c_profile = temp_profile_700.values - 273.15 # shape: [levels, y, x]

# Identify any layer above freezing vs. entire below freezing (Bourgouin approach)

above_freezing_any = np.any(temp_c_profile > 0, axis=0)

below_freezing_all = np.all(temp_c_profile <= 0, axis=0)

# Initialize precipitation type array

ptype_hybrid = np.full(t2_da.shape, "None", dtype=object)

# ----------------------------------------------------------

# 1) Preliminary "Raw" Classifications from Microphysics

# ----------------------------------------------------------

# If there's a significant mixing ratio of rain, snow, or graupel at surface level

is_rain_surf = (qrain_surf > 1e-6)

is_snow_surf = (qsnow_surf > 1e-6)

is_graup_surf = (qgraup_surf > 1e-6)

# ----------------------------------------------------------

# 2) Ramer/Baldwin Criteria (simplified)

# ----------------------------------------------------------

# Ramer-like: if temperature above 0°C in upper slices but <=0°C at surface => freezing rain

# Baldwin-like logic can also consider thickness/delta T, but here simplified

level_temp_high = temp_profile_700.isel(bottom_top=-1).values - 273.15

surface_temp_c = t2_c

ramer_freezing_rain = (level_temp_high > 0) & (surface_temp_c <= 0)

ramer_snow = (surface_temp_c <= 0)

ramer_rain = (surface_temp_c > 0)

# ----------------------------------------------------------

# 3) Bourgouin Warm/Cold Layer Detection (Simplified)

# ----------------------------------------------------------

# Warming layer above, subfreezing at surface => freezing rain or sleet

warm_layer_bourgouin = above_freezing_any

cold_surface_bourgouin = (surface_temp_c <= 0)

bourgouin_freezing_rain = warm_layer_bourgouin & cold_surface_bourgouin

# ----------------------------------------------------------

# 4) Ensemble Voting

# ----------------------------------------------------------

votes_rain = np.zeros_like(t2_da.values, dtype=int)

votes_freezing_rain = np.zeros_like(t2_da.values, dtype=int)

votes_snow = np.zeros_like(t2_da.values, dtype=int)

votes_sleet = np.zeros_like(t2_da.values, dtype=int)

# Based on microphysics presence

votes_snow[is_snow_surf] += 1

votes_sleet[is_graup_surf] += 1

votes_rain[is_rain_surf & (surface_temp_c > 0)] += 1

votes_freezing_rain[is_rain_surf & (surface_temp_c <= 0)] += 1

# Based on Ramer approach

votes_freezing_rain[ramer_freezing_rain] += 1

votes_snow[ramer_snow] += 1

votes_rain[ramer_rain] += 1

# Based on Bourgouin

votes_freezing_rain[bourgouin_freezing_rain] += 1

# Decide final classification by majority ptype_hybrid[votes_freezing_rain >= 2] = "Freezing Rain"

ptype_hybrid[votes_snow >= 2] = "Snow"

ptype_hybrid[votes_sleet >= 2] = "Sleet/Graupel"

ptype_hybrid[votes_rain >= 2] = "Rain"

# ----------------------------------------------------------

# 5) Resolve Ambiguous Grid Cells with Wet-Bulb Temperature

# ----------------------------------------------------------

ambiguous_mask = (ptype_hybrid == "None")

ptype_hybrid[ambiguous_mask & (wb_c < 0)] = "Freezing Rain"

ptype_hybrid[ambiguous_mask & (wb_c >= 0) & (wb_c <= 2)] = "Rain"

ptype_hybrid[ambiguous_mask & (wb_c > 2)] = "Rain"

# If no surface microphysics present, set to "None"

no_precip_mask = (~is_rain_surf & ~is_snow_surf & ~is_graup_surf)

ptype_hybrid[no_precip_mask] = "None"

return ptype_hybrid

###################################################################

# Vectorized Function to Calculate Road Icing Index (RII)

###################################################################

def calculate_icing_index_vectorized(t2k, wbk, ptype, precip_rate, rh, ws, top_soil_temp_k):

"""

Compute Road Icing Index (RII) in a vectorized manner.

Parameters

----------

t2k : ndarray

Near-surface air temperature in Kelvin.

wbk : ndarray

Wet-bulb temperature in Kelvin.

ptype : ndarray

Precipitation type classification (e.g., 'Rain', 'Snow', 'Freezing Rain', 'Sleet/Graupel', 'None').

precip_rate : ndarray

Precipitation rate in mm/h.

rh : ndarray

Relative humidity in %.

ws : ndarray

Wind speed in m/s.

top_soil_temp_k : ndarray

Top-layer soil temperature in Kelvin.

Returns

-------

icing_index : ndarray

Calculated Road Icing Index for each grid cell.

# Convert temperatures from Kelvin to Celsius

t2_c = t2k - 273.15

wb_c = wbk - 273.15

soil_c = top_soil_temp_k - 273.15

# Temperature Score

temp_score = np.zeros_like(t2_c)

temp_score[t2_c < -5] = 5 # Change from 3 to 4 to further increase the probability of freezing taking into account probable underestimation of the model

temp_score[(t2_c >= -5) & (t2_c <= 0)] = 3 # Change from 2 to 3 to further increase the probability of freezing taking into account probable underestimation of the model

temp_score[(t2_c > 0) & (t2_c <= 3)] = 2

temp_score[(t2_c > 3 ) & (t2_c <= 5)] = -1 # These lines are added to penalize and reduce the probability of freezing when there are high temperatures

temp_score[(t2_c > 5) & (t2_c <= 10)] = -3

temp_score[t2_c > 10] = -5

# Wet-Bulb Temperature Score

wet_bulb_score = np.zeros_like(wb_c)

wet_bulb_score[wb_c < 0] = 2

wet_bulb_score[(wb_c >= 0) & (wb_c <= 2)] = 1

# Precipitation Type Score

ptype_score_map = {

"Rain": 1,

"Freezing Rain": 3,

"Snow": 2,

"Sleet/Graupel": 2,

"None": 0

}

ptype_score = np.zeros_like(t2_c)

for pt, val in ptype_score_map.items():

ptype_score[ptype == pt] = val

# Precipitation Rate Score

precip_rate_score = np.zeros_like(precip_rate)

# Masks for precipitation type

freezing_rain_mask = (ptype == "Freezing Rain")

snow_mask = (ptype == "Snow") | (ptype == "Sleet/Graupel")

rain_mask = (ptype == "Rain")

# Freezing Rain thresholds

precip_rate_score[freezing_rain_mask & (precip_rate <= 2.5)] = 1

precip_rate_score[freezing_rain_mask & (precip_rate > 2.5) & (precip_rate <= 7.5)] = 2

precip_rate_score[freezing_rain_mask & (precip_rate > 7.5)] = 3

# Snow thresholds (liquid equivalent)

precip_rate_score[snow_mask & (precip_rate <= 1.0)] = 1

precip_rate_score[snow_mask & (precip_rate > 1.0) & (precip_rate <= 3.0)] = 2

precip_rate_score[snow_mask & (precip_rate > 3.0)] = 3

# Rain thresholds

precip_rate_score[rain_mask & (precip_rate <= 5.0)] = 1

precip_rate_score[rain_mask & (precip_rate > 5.0) & (precip_rate <= 15.0)] = 2

precip_rate_score[rain_mask & (precip_rate > 15.0)] = 3

# Relative Humidity Score

rh_score = np.zeros_like(rh)

rh_score[rh >= 90] = 1

# Wind Speed Score

wind_score = np.zeros_like(ws)

wind_score[ws < 10] = 1

wind_score[(ws >= 10) & (ws <= 20)] = 0

wind_score[ws > 20] = -1

# Soil Temperature Score

soil_score = np.zeros_like(soil_c)

soil_score[soil_c < 0] = 1

# Summation for Road Icing Index

# Each factor is added, reflecting the relative contribution to icing risk.

icing_index = (

temp_score

+ wet_bulb_score

+ ptype_score

+ precip_rate_score

+ rh_score

+ wind_score

+ soil_score

)

return icing_index,temp_score,wet_bulb_score,ptype_score,precip_rate_score,rh_score,wind_score,soil_score

This python code demonstrates how the RII calculated and implemented as well as how the hybrid ensemble voting equations.