Beam Blockage for SAIL

import pyart
import wradlib as wrl
import cartopy.crs as ccrs
import numpy as np
import matplotlib.pyplot as plt

radar = pyart.io.read('/home/zsherman/Downloads/xprecipradar_guc_volume_20220217-230559.b1.nc')

radar.fields.keys()

display = pyart.graph.RadarMapDisplay(radar)

fig = plt.figure(figsize=(10, 10))
display.plot_ppi_map('DBZ', 0)

radar.altitude['data']

radar.nsweeps

tif_file = '/home/zsherman/Downloads/sail.tif'
#tif_file  = '/home/zsherman/Downloads/dsm_mosaic_min_phase_me.tif'
#tif_file = '/home/zsherman/Downloads/tch_mosaic_min_phase_me.tif'

def beam_block(radar, tif_file, radar_height_offset=10.0,
               beam_width=1.0):
    """
    Beam Block Radar Calculation.
    Parameters
    ----------
    radar : Radar
        Radar object used.
    tif_name : string
        Name of geotiff file to use for the
        calculation.
    radar_height_offset : float
        Add height to the radar altitude for radar towers.
    Other Parameters
    ----------------
    beam_width : float
        Radar's beam width for calculation.
        Default value is 1.0.
    Returns
    -------
    pbb_all : array
        Array of partial beam block fractions for each
        gate in all sweeps.
    cbb_all : array
        Array of cumulative beam block fractions for
        each gate in all sweeps.
    References
    ----------
    Bech, J., B. Codina, J. Lorente, and D. Bebbington,
    2003: The sensitivity of single polarization weather
    radar beam blockage correction to variability in the
    vertical refractivity gradient. J. Atmos. Oceanic
    Technol., 20, 845–855
    Heistermann, M., Jacobi, S., and Pfaff, T., 2013:
    Technical Note: An open source library for processing
    weather radar data (wradlib), Hydrol. Earth Syst.
    Sci., 17, 863-871, doi:10.5194/hess-17-863-2013
    Helmus, J.J. & Collis, S.M., (2016). The Python ARM
    Radar Toolkit (Py-ART), a Library for Working with
    Weather Radar Data in the Python Programming Language.
    Journal of Open Research Software. 4(1), p.e25.
    DOI: http://doi.org/10.5334/jors.119
    """
    # Opening the tif file and getting the values ready to be
    # converted into polar values.
    rasterfile = tif_file
    data_raster = wrl.io.open_raster(rasterfile)
    rastervalues, rastercoords, proj = wrl.georef.extract_raster_dataset(
        data_raster, nodata=None)
    #rastervalues_, rastercoords_, proj = wrl.georef.extract_raster_dataset(data_raster, nodata=-32768.)
    sitecoords = (np.float(radar.longitude['data']),
                  np.float(radar.latitude['data']),
                  np.float(radar.altitude['data'] + radar_height_offset))
    pbb_arrays = []
    cbb_arrays = []
    _range = radar.range['data']
    beamradius = wrl.util.half_power_radius(_range, beam_width)
    # Cycling through all sweeps in the radar object.
    print('Calculating beam blockage.')
    del data_raster
    for i in range(len(radar.sweep_start_ray_index['data'])):
        index_start = radar.sweep_start_ray_index['data'][i]
        index_end = radar.sweep_end_ray_index['data'][i] + 1
        elevs = radar.elevation['data'][index_start:index_end]
        azimuths = radar.azimuth['data'][index_start:index_end]
        rg, azg = np.meshgrid(_range, azimuths)
        rg, eleg = np.meshgrid(_range, elevs)
        nrays = azimuths.shape[0]              # number of rays
        nbins = radar.ngates                   # number of range bins
        bw = beam_width                        # half power beam width (deg)
        range_res = 60.                       # range resolution (meters)
        el = radar.fixed_angle['data'][i]
        coord = wrl.georef.sweep_centroids(nrays, range_res, nbins, el)
        coords = wrl.georef.spherical_to_proj(rg, azg, eleg,
                                              sitecoords, proj=proj)
        lon = coords[..., 0]
        lat = coords[..., 1]
        alt = coords[..., 2]
        polcoords = coords[..., :2]
        rlimits = (lon.min(), lat.min(), lon.max(), lat.max())

        #Clip the region inside our bounding box
        ind = wrl.util.find_bbox_indices(rastercoords, rlimits)
        rastercoords = rastercoords[ind[0]:ind[3], ind[0]:ind[2], ...]
        rastervalues = rastervalues[ind[0]:ind[3], ind[0]:ind[2]]
        polarvalues = wrl.ipol.cart_to_irregular_spline(
            rastercoords, rastervalues, polcoords, order=3,
            prefilter=False)
        # Calculate partial beam blockage using wradlib.
        pbb = wrl.qual.beam_block_frac(polarvalues, alt, beamradius)
        pbb = np.ma.masked_invalid(pbb)
        pbb[pbb < 0] = 0.0
        pbb[pbb > 1] = 1.0
        pbb_arrays.append(pbb)
        # Calculate cumulative beam blockage using wradlib.
        cbb = wrl.qual.cum_beam_block_frac(pbb)
        cbb_arrays.append(cbb)
    pbb_all = np.ma.concatenate(pbb_arrays)
    cbb_all = np.ma.concatenate(cbb_arrays)
    return pbb_all, cbb_all

def pbb_to_dict(pbb_all):
    """ Function that takes the pbb_all array and turns
    it into a dictionary to be used and added to the
    pyart radar object. """
    pbb_dict = {}
    pbb_dict['coordinates'] = 'elevation azimuth range'
    pbb_dict['units'] = '1'
    pbb_dict['data'] = pbb_all
    pbb_dict['long_name'] = 'Partial Beam Block Fraction'
    pbb_dict['comment'] = 'Partial beam block fraction due to terrain.'
    return pbb_dict


def cbb_to_dict(cbb_all):
    """ Function that takes the cbb_all array and turns
    it into a dictionary to be used and added to the
    pyart radar object. """
    cbb_dict = {}
    cbb_dict['coordinates'] = 'elevation azimuth range'
    cbb_dict['units'] = '1'
    cbb_dict['data'] = cbb_all
    cbb_dict['long_name'] = 'Cumulative Beam Block Fraction'
    cbb_dict['comment'] = 'Cumulative beam block fraction due to terrain.'
    return cbb_dict


def cbb_flags_to_dict(cbb_flags):
    """ Function that takes the cbb_flag array and turns
    it into a dictionary to be used and added to the
    pyart radar object. """
    cbb_dict = {}
    cbb_dict['coordinates'] = 'elevation azimuth range'
    cbb_dict['units'] = '1'
    cbb_dict['data'] = cbb_flags
    cbb_dict['long_name'] = 'Cumulative Beam Block Fraction Flag'
    cbb_dict['comment'] = 'Cumulative beam block flag due to terrain.'
    return cbb_dict

radar.range

radar.range['data'] = radar.range['data'] + (59.94095*2)

radar.init_gate_altitude()
radar.init_gate_longitude_latitude()

radar.range['data'][0:4]

pbb, cbb = beam_block(radar, tif_file, radar_height_offset=10.0, beam_width=1.0)

pbb_dict = pbb_to_dict(pbb)
cbb_dict = cbb_to_dict(cbb)
radar.add_field('partial_beam_blockage', pbb_dict, replace_existing=True)
radar.add_field('cumulative_beam_blockage', cbb_dict, replace_existing=True)

radar.fields['partial_beam_blockage']['data'].min()

fig = plt.figure(figsize=(10, 8))
display.plot_ppi('partial_beam_blockage', 0, cmap='pyart_HomeyerRainbow')
plt.plot(0, 0, marker="o", markersize=20, markeredgecolor="red", markerfacecolor="green")
plt.xlim(-5, 5)
plt.ylim(-5, 5)
plt.savefig('/home/zsherman/pbb_sail_zoomed.png')

pbb[:, 0:5]

fig = plt.figure(figsize=(10, 8))
display.plot_ppi('cumulative_beam_blockage', 2, cmap='pyart_HomeyerRainbow')
plt.savefig('/home/zsherman/cbb.png')

def beam_block_flags(cbb_all, complete_block_thresh=0.95):
    """
    Takes CBB arrays created from the beam_block function and creates beam block flags.

    Parameters
    ----------
    cbb_all : array
        Array of cumulative beam blockage values created from the beam_block
        function.

    Other Parameters
    ----------------
    complete_block_thresh : float
        Value for the cutoff for complete blockage flag value of 1. Anything
        above the complete_block_thres is considered blocked. Default value
        is 0.95.

    Returns
    -------
    cbb_flags : array
        Array of integer values depicting no, partial, and complete beam
        blockage based on the cumulative beam blockage data. This array can
        then be used to create a cumulative beam block flags field.
    """
    cbb_flags = np.ones(cbb_all.shape)
    cbb_flags[cbb_all < complete_block_thresh] = 0
    return cbb_flags

cbb_flags = beam_block_flags(cbb, complete_block_thresh=0.80)

cbb_flag_dict = cbb_flags_to_dict(cbb_flags)

radar.add_field('cbb_flag', cbb_flag_dict, replace_existing=True)

fig = plt.figure(figsize=(10, 8))
display.plot_ppi('cbb_flag', 2, cmap='pyart_HomeyerRainbow')
plt.savefig('/home/zsherman/cbb_flag.png')

fields = []
for field in radar.fields.keys():
    fields.append(field)

fields

for field in fields:
    if field == 'cbb_flag':
        continue
    elif field == 'cumulative_beam_blockage':
        continue
    else:
        radar.fields.pop(field)

radar.fields.keys()

pyart.io.write_cfradial('/home/zsherman/flag_radar.nc', radar)