This notebook is part of the \(\omega radlib\) documentation: https://docs.wradlib.org.

Beam Blockage Calculation using a DEM#

Here, we derive (partial) beam-blockage (PBB) from a Digital Elevation Model (DEM).

We require

the local radar setup (sitecoords, number of rays, number of bins, antenna elevation, beamwidth, and the range resolution);
a DEM with a adequate resolution.

Here we use pre-processed data from the GTOPO30 and SRTM missions.

import datetime as dt
import os

import matplotlib as mpl
import matplotlib.pyplot as plt
import numpy as np
import wradlib as wrl
import xarray as xr

import radarpolarimetry as rp

Calculate Beamblockage#

Now we can finally apply the wradlib.qual.beam_block_frac function to calculate the PBB.

beamradius = wrl.util.half_power_radius(swp.range, bw)
PBB = wrl.qual.beam_block_frac(swp.DEM.values, swp.z.values, beamradius)
PBB = np.ma.masked_invalid(PBB)

So far, we calculated the fraction of beam blockage for each bin.

But we need to into account that the radar signal travels along a beam. Cumulative beam blockage (CBB) in one bin along a beam will always be at least as high as the maximum PBB of the preceding bins (see wradlib.qual.cum_beam_block_frac)

CBB = wrl.qual.cum_beam_block_frac(PBB)

swp = swp.assign(
    CBB=(["azimuth", "range"], CBB),
    PBB=(["azimuth", "range"], PBB),
)

Visualize Beamblockage#

Now we visualize

the average terrain altitude per radar bin
a beam blockage map
interaction with terrain along a single beam

# just a little helper function to style x and y axes of our maps


def annotate_map(ax, cm=None, title=""):
    ticks = (ax.get_xticks() / 1000).astype(int)
    ax.set_xticklabels(ticks)
    ticks = (ax.get_yticks() / 1000).astype(int)
    ax.set_yticklabels(ticks)
    ax.set_xlabel("Kilometers")
    ax.set_ylabel("Kilometers")
    # if not cm is None:
    #    plt.colorbar(cm, ax=ax)
    if not title == "":
        ax.set_title(title)
    ax.grid()

fig = plt.figure(figsize=(15, 12))

# create subplots
ax1 = plt.subplot2grid((2, 2), (0, 0))
ax2 = plt.subplot2grid((2, 2), (0, 1))
ax3 = plt.subplot2grid((2, 2), (1, 0), colspan=2, rowspan=1)

# azimuth angle
angle = 225

# Plot terrain (on ax1)
dem_pm = swp.DEM.wrl.vis.plot(ax=ax1, cmap=mpl.cm.terrain, vmin=0.0, add_colorbar=False)
swp.sel(azimuth=angle, method="nearest").plot.scatter(
    x="x", y="y", marker=".", color="r", alpha=0.2, ax=ax1
)
ax1.plot(swp.longitude.values, swp.latitude.values, "ro")
annotate_map(
    ax1,
    dem_pm,
    f"Terrain within {np.max(swp.range.values / 1000.0) + 0.1} km range",
)

# Plot CBB (on ax2)
cbb = swp.CBB.wrl.vis.plot(ax=ax2, cmap=mpl.cm.PuRd, vmin=0, vmax=1, add_colorbar=False)
annotate_map(ax2, cbb, "Beam-Blockage Fraction")

# Plot single ray terrain profile on ax3
(bc,) = ax3.plot(
    swp.range / 1000.0, swp.z[angle, :], "-b", linewidth=3, label="Beam Center"
)
(b3db,) = ax3.plot(
    swp.range.values / 1000.0,
    (swp.z[angle, :] + beamradius),
    ":b",
    linewidth=1.5,
    label="3 dB Beam width",
)
ax3.plot(swp.range / 1000.0, (swp.z[angle, :] - beamradius), ":b")
ax3.fill_between(swp.range / 1000.0, 0.0, swp.DEM[angle, :], color="0.75")
ax3.set_xlim(0.0, np.max(swp.range / 1000.0) + 0.1)
ax3.set_ylim(0.0, 3000)
ax3.set_xlabel("Range (km)")
ax3.set_ylabel("Altitude (m)")
ax3.grid()

axb = ax3.twinx()
(bbf,) = axb.plot(swp.range / 1000.0, swp.CBB[angle, :], "-k", label="BBF")
axb.set_ylabel("Beam-blockage fraction")
axb.set_ylim(0.0, 1.0)
axb.set_xlim(0.0, np.max(swp.range / 1000.0) + 0.1)


legend = ax3.legend(
    (bc, b3db, bbf),
    ("Beam Center", "3 dB Beam width", "BBF"),
    loc="upper left",
    fontsize=10,
)

../../_images/c5e1df031a1fad50ea5bf0114ae1e45ecd623854fa35c3d8fe8e75f97defc584.png

Visualize Beam Propagation showing earth curvature#

Now we visualize

interaction with terrain along a single beam

In this representation the earth curvature is shown. For this we assume the earth a sphere with exactly 6370000 m radius. This is needed to get the height ticks at nice position.

def height_formatter(x, pos):
    x = (x - 6370000) / 1000
    fmt_str = f"{x:g}"
    return fmt_str


def range_formatter(x, pos):
    x = x / 1000.0
    fmt_str = f"{x:g}"
    return fmt_str

The wradlib.vis.create_cg-function is facilitated to create the curved geometries.
The actual data is plottet as (theta, range) on the parasite axis.
Some tweaking is needed to get the final plot look nice.

fig = plt.figure(figsize=(10, 6))

cgax, caax, paax = wrl.vis.create_cg(fig=fig, rot=0, scale=1)

# azimuth angle
angle = 240
sitecoords = (swp.longitude.values, swp.latitude.values, swp.altitude.values)
r = swp.range.values
az = swp.azimuth.values

alt = swp.z.values


# fix grid_helper
er = 6370000
gh = cgax.get_grid_helper()
gh.grid_finder.grid_locator2._nbins = 80
gh.grid_finder.grid_locator2._steps = [1, 2, 4, 5, 10]

# calculate beam_height and arc_distance for ke=1
# means line of sight
bhe = wrl.georef.bin_altitude(r, 0, sitecoords[2], re=er, ke=1.0)
ade = wrl.georef.bin_distance(r, 0, sitecoords[2], re=er, ke=1.0)
nn0 = np.zeros_like(r)
# for nice plotting we assume earth_radius = 6370000 m
ecp = nn0 + er
# theta (arc_distance sector angle)
thetap = -np.degrees(ade / er) + 90.0

# zero degree elevation with standard refraction
bh0 = wrl.georef.bin_altitude(r, 0, sitecoords[2], re=er)

# plot (ecp is earth surface normal null)
(bes,) = paax.plot(thetap, ecp, "-k", linewidth=3, label="Earth Surface NN")
(bc,) = paax.plot(thetap, ecp + alt[angle, :], "-b", linewidth=3, label="Beam Center")
(bc0r,) = paax.plot(thetap, ecp + bh0, "-g", label="0 deg Refraction")
(bc0n,) = paax.plot(thetap, ecp + bhe, "-r", label="0 deg line of sight")
(b3db,) = paax.plot(
    thetap, ecp + alt[angle, :] + beamradius, ":b", label="+3 dB Beam width"
)
paax.plot(thetap, ecp + alt[angle, :] - beamradius, ":b", label="-3 dB Beam width")

# orography
paax.fill_between(thetap, ecp, ecp + swp.DEM[angle, :], color="0.75")

# shape axes
cgax.set_xlim(0, np.max(ade))
cgax.set_ylim([ecp.min() - 1000, ecp.max() + 2500])
caax.grid(True, axis="x")
cgax.grid(True, axis="y")
cgax.axis["top"].toggle(all=False)
caax.yaxis.set_major_locator(
    mpl.ticker.MaxNLocator(steps=[1, 2, 4, 5, 10], nbins=20, prune="both")
)
caax.xaxis.set_major_locator(mpl.ticker.MaxNLocator())
caax.yaxis.set_major_formatter(mpl.ticker.FuncFormatter(height_formatter))
caax.xaxis.set_major_formatter(mpl.ticker.FuncFormatter(range_formatter))

caax.set_xlabel("Range (km)")
caax.set_ylabel("Altitude (km)")

legend = paax.legend(
    (bes, bc0n, bc0r, bc, b3db),
    (
        "Earth Surface NN",
        "0 deg line of sight",
        "0 deg std refraction",
        "Beam Center",
        "3 dB Beam width",
    ),
    loc="lower left",
    fontsize=10,
)

../../_images/9bf0f6cf1023b8210a7c05e8fe3c357d50ea648c51cd1a3e3ce0754e31383118.png

Go back to Read DEM Raster Data, change the rasterfile to use the other resolution DEM and process again.

Beam Blockage Calculation using a DEM

Contents

Beam Blockage Calculation using a DEM#

Setup for Bonn radar#

Preprocessing the digitial elevation model#

Calculate Beamblockage#

Visualize Beamblockage#

Visualize Beam Propagation showing earth curvature#