import numpy as np
import xarray as xr
import pyart
def column_profile(radar: pyart.core.Radar,
latitude: float = 33.8,
longitude: float = -88.3,
azimuth_spread: float = 3,
spatial_spread: float = 3,
v_res: float = 100,
min_alt: float = None,
max_alt: float = None) -> xr.Dataset:
"""
Computes a vertical column profile from radar data and interpolates it
onto a uniform height grid.
Parameters
----------
radar : pyart.core.Radar
Py-ART radar object containing the volume scan data.
latitude : float, optional
Latitude of the point of interest, by default 33.8.
longitude : float, optional
Longitude of the point of interest, by default -88.3.
azimuth_spread : float, optional
Azimuthal spread in degrees to consider around the point of interest, by default 3.
spatial_spread : float, optional
Horizontal spatial spread in kilometers for averaging, by default 3.
v_res : float, optional
Vertical resolution in meters, by default 100.
min_alt : float
Minimum altitude (height) in meters for the vertical profile.
max_alt : float
Maximum altitude (height) in meters for the vertical profile.
Returns
-------
xr.Dataset
Interpolated columnar vertical profile on a uniform height grid.
"""
if min_alt is None or max_alt is None:
raise ValueError("Both min_alt and max_alt must be specified.")
# Compute the column profile using Py-ART
col_prof = pyart.util.column_vertical_profile(
radar, latitude=latitude, longitude=longitude,
azimuth_spread=azimuth_spread, spatial_spread=spatial_spread
)
# Define new height levels for interpolation
new_heights = xr.Dataset(
coords={'height': ('height', np.arange(min_alt, max_alt + v_res, v_res))}
)
# Interpolate the column profile onto the uniform height grid
col_prof_interp = col_prof.interp_like(new_heights)
return col_prof_interp
def rain_rate_mp(dbz, a=200.0, b=1.6):
"""Calculates the inverse of input decibel values
Convert to rain rate using Marshall Palmer Relation
"""
Z = 10.0 ** (dbz / 10.0)
return (Z / a) ** (1.0 / b)
def rain_rate_categorized(dbz, conv_threshold=40.0, a_conv=300.0, b_conv=1.4,
a_strat=200.0, b_strat=1.6):
"""
Calculates rain rate by categorizing reflectivity into convective and stratiform regions,
then applying different coefficients for each category using the Power-law relation.
Parameters:
dbz (xarray.DataArray): Input reflectivity values in dBZ.
conv_threshold (float): Threshold dBZ value for convective classification. Default is 40.0 dBZ.
a_conv (float): Coefficient 'a' for convective rain. Default is 300.0.
b_conv (float): Exponent 'b' for convective rain. Default is 1.4.
a_strat (float): Coefficient 'a' for stratiform rain. Default is 200.0.
b_strat (float): Exponent 'b' for stratiform rain. Default is 1.6.
Returns:
xarray.DataArray: Rain rate in mm/h, with the same dimensions and coordinates as the input dbz.
"""
# Convert dBZ to reflectivity factor Z
Z = 10.0 ** (dbz / 10.0)
# Classify into convective and stratiform regions
is_conv = dbz >= conv_threshold
# Initialize rain rate array
rr = xr.zeros_like(Z)
# Calculate convective rain rate
rr = xr.where(is_conv, (Z / a_conv) ** (1.0 / b_conv), rr)
# Calculate stratiform rain rate
rr = xr.where(~is_conv, (Z / a_strat) ** (1.0 / b_strat), rr)
# Preserve metadata and coordinates
rr.attrs['standard_name'] = "Rain rate"
rr.name = "rain_rate"
rr.attrs["units"] = "mm/h"
rr.attrs["description"] = "Rain rate calculated using categorized reflectivity"
return rr
import glob
files = sorted(glob.glob("path/to/your/data*.nc"))
print(len(files))
rain_rate_list = []
for file in files[-5:]:
radar = pyart.io.read(file)
col_prof_interp = column_profile(radar,
latitude=33.5, longitude=-88.3,
azimuth_spread=3, spatial_spread=3,
v_res=100, min_alt=200, max_alt=8e3)
rain_rate = rain_rate_categorized(col_prof_interp['REF'])
rain_rate['time'] = col_prof_interp['base_time']
rain_rate_list.append(rain_rate)
da_rain_rate = xr.concat(rain_rate_list, dim='time')
da_rain_rate.plot.contourf(x='time', levels=range(0, 55, 5))
import os
os.makedirs("~/Downloads/output_qpe_folder", exist_ok=True)
da_rain_rate.to_netcdf("~/Downloads/output_qpe_folder/rain_rates_vert.nc")
