Help/guidance on cumulative beam blockage

[zssherman]

Hello Wradlib folks!

wradlib=1.16 python=3.9

I was doing some beam blockage calculations with my code here:

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)
        print(pbb.min())
        print(pbb.max())
        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

The partial beam blockage looks as expected, images attached.

However when I do cumulative beam blockage, its all blocked. The radar has no blockage to the north until certain points to the northwest and northeast so I would expect CBB rays to the south and eventually starting to the north.

Was wondering, if there is something I'm missing? Its been awhile since i've worked with the beam blockage code, but checking to see if I have cumulative calculated incorrectly.

[kmuehlbauer]

@zssherman I'm assuming the green dot with red border in the image center is just a indication of the radar location and is not within the data, right? Maybe you can just leave it out to see the data behind it.

I suspect that your pbb data has it's maximum in the very first bin. This would lead to complete blocked output. If that's the case we need to look at your source data of ranges. If the first range bin is 0 (zero), the pbb will be 1. That's because of a division by zero (half-power beam radius would be 0) which in turn leads to an inf value which is treated as full blocked value in the pbb calculation.

[zssherman]

Hi @kmuehlbauer ! The green dot was just to indicate where the radar was, to make sure it wasn't being placed in the mountain, due to incorrect coordinate info in the radar file. I'll check the radar range now, and get back to you here shortly, thanks again for your guidance on this!

[zssherman]

@kmuehlbauer Possible problem, the first radar range is -112 for this radar: 'meters_to_center_of_first_gate': -112.830795, 'meters_between_gates': 59.94095

I'm wondering if this is the cause of everything then as you mentioned the radar range causing inf value

[zssherman]

That fixed it, fixing the offset of range!

[zssherman]

Thank you @kmuehlbauer for the help with this! I'll close this for now!

[kmuehlbauer]

That fixed it, fixing the offset of range!

Great you figured it out and my gut feeling was correct. But, what a strange place for a radar. Could you comment on that special location?

[zssherman]

@kmuehlbauer Its a part of ARM's SAIL campaign found here: https://www.arm.gov/research/campaigns/amf2021sail

Gives a background on the region chosen, but for that specific area for the radar, i'm not 100% sure to be honest. I'm going to ask and see if I can give a better answer

[kmuehlbauer]

Thanks Zach, I'll have a closer look later. But that already looks very much like mountainous river basin.

[zssherman]

@kmuehlbauer More info from Adam from ARM: This view gave us the most coverage over the main ARM site, but also to the west as well to support the NOAA SPLASH radar and possible multi-doppler efforts. In the end, this study is really just looking at that water basin and not the broader region per se. Chandra did the beam blockage analysis for a number of sites on that mountain which is why we decided to add the riser to increase the height even more. It really was an in depth process by experts

[kmuehlbauer]

Science rocks!