Compositing with BALTRAD

This exercise builds on output from the parallel processing exercise. It does not address how projections and navigation is dealt with in BALTRAD. This should be addressed in a separate exercise.

The Cartesian product area used in this exercise is pre-configured and looked up from a registry.

Rudimentary composite

%matplotlib inline
import glob, time
import matplotlib
import _raveio, _rave
import _pycomposite, compositing
import warnings
warnings.filterwarnings('ignore')  # Suppress SyntaxWarning from Python2 code

generator = compositing.compositing()
generator.filenames = glob.glob("data/se*.h5")

# Run with all defaults to a pre-configured area that uses the Google Maps projection.
# First two arguments are product date and time. These are taken from the last input file if not specified.
before = time.time()
comp = generator.generate(None, None, area="swegmaps_2000")
after = time.time()

rio = _raveio.new()
rio.object = comp
rio.save("data/comp_pcappi1000m.h5")

print("Compositing took %3.2f seconds" % (after-before))

Compositing took 12.78 seconds

Tweak the plotter from earlier exercises

# Two color palettes, one used in GoogleMapsPlugin, and the other from RAVE
from GmapColorMap import dbzh as dbzp

# Convert a 768-list palette to a matplotlib colorlist
def make_colorlist(pal):
    colorlist = []
    for i in range(0, len(pal), 3):
        colorlist.append([pal[i]/255.0, pal[i+1]/255.0, pal[i+2]/255.0])
    return colorlist

# Convert lists to colormaps
dbzcl = make_colorlist(dbzp)

# Then create a simple plotter
import matplotlib.pyplot as plt
StringType = type('')
def plot(data, colorlist=dbzcl, title="Composite"):
    mini, maxi = data.shape.index(min(data.shape)), data.shape.index(max(data.shape))
    figsize=(20,16)# if mini == 0 else (12,8)
    fig = plt.figure(figsize=figsize)
    plt.title(title)
    clist=colorlist if type(colorlist)==StringType else matplotlib.colors.ListedColormap(colorlist)
    plt.imshow(data, cmap=clist, clim=(0,255))
    plt.colorbar(shrink=float(data.shape[mini])/data.shape[maxi])

plot(comp.getParameter("DBZH").getData(), title="Default composite: DBZH 1000 m Pseudo-CAPPI, nearest radar")

../../_images/BALTRAD Compositing_7_0.png

Maximum reflectivity, lowest pixel, add QC chain

generator.product = _rave.Rave_ProductType_MAX
generator.selection_method = _pycomposite.SelectionMethod_HEIGHT
generator.detectors = ["ropo", "beamb", "radvol-att", "radvol-broad", "rave-overshooting", "qi-total"]
before = time.time()
comp = generator.generate(None, None, area="swegmaps_2000")
after = time.time()
rio.object = comp
rio.save("data/comp_max.h5")
print("Compositing took %3.2f seconds" % (after-before))

Compositing took 75.97 seconds

plot(comp.getParameter("DBZH").getData(), title="Maximum reflectivity, lowest pixel")

../../_images/BALTRAD Compositing_10_0.png

Plot correspondong total quality index

dbzh = comp.getParameter("DBZH")
qitot = dbzh.getQualityFieldByHowTask("pl.imgw.quality.qi_total")
plot(qitot.getData(), "binary", "Total quality index")

../../_images/BALTRAD Compositing_12_0.png

Now use “total quality” as the compositing criterion

generator.qitotal_field = "pl.imgw.quality.qi_total"
before = time.time()
comp = generator.generate(None, None, area="swegmaps_2000")
after = time.time()
rio.object = comp
rio.save("data/comp_qitotal.h5")
print("Compositing took %3.2f seconds" % (after-before))

Compositing took 26.68 seconds

plot(comp.getParameter("DBZH").getData(), title="Maximum reflectivity, quality-based")

../../_images/BALTRAD Compositing_15_0.png

plot(comp.getParameter("DBZH").getQualityFieldByHowTask("pl.imgw.quality.qi_total").getData(), "binary", "Total quality index")

../../_images/BALTRAD Compositing_16_0.png