Quantitative Precipitation Forecast (QPF)

This example demonstrates how to perform QPF using radar data

import xarray as xr
import os.path
import numpy as np
import os
from pyresample import utils
from matplotlib.colors import BoundaryNorm
import matplotlib
import matplotlib.pyplot as plt

from swirlspy.qpf.rover import rover
from swirlspy.rad.iris import read_iris

try:
    THIS_DIR = os.path.dirname(os.path.abspath(__file__))
except NameError:
    import sys
    THIS_DIR = os.path.dirname(os.path.abspath(sys.argv[0]))

Initialising

1. Define target grid as a pyresample AreaDefinition.
2. Read radar data files with read_iris().

# Defining target grid as a pyresample AreaDefinition.
area_id = "hk1980_250km"
description = (
    "A 250 m resolution rectangular grid "
    "centred at HKO and extending to 250 km "
    "in each direction in HK1980 easting/northing coordinates"
    )
proj_id = 'hk1980'
projection = (
    '+proj=tmerc +lat_0=22.31213333333334 '
    '+lon_0=114.1785555555556 +k=1 +x_0=836694.05 '
    '+y_0=819069.8 +ellps=intl +towgs84=-162.619,-276.959,'
    '-161.764,0.067753,-2.24365,-1.15883,-1.09425 +units=m '
    '+no_defs'
    )
x_size = 500
y_size = 500
area_extent = (587000, 569000, 1087000, 1069000)
area_def = utils.get_area_def(
    area_id, description, proj_id, projection, x_size, y_size, area_extent
)

# Read data from files (in this example, the files are
# in iris format) and storing data as xarrays.

valid_time = '201805150624'
steps = 6
time_step_size = 6

file_name1 = "../tests/samples/TMS170524110001.RAW2SGJ"
file_name2 = "../tests/samples/TMS170524110619.RAW2SJY"

# Reading file data, storing as xarray
xarray1 = read_iris(
    file_name1,
    area_def=area_def,
    coord_label=['easting', 'northing'],
    radar_radius=250000
)
xarray2 = read_iris(
    file_name2,
    area_def=area_def,
    coord_label=['easting', 'northing'],
    radar_radius=250000
)
qpf_xarray = xr.concat([xarray1, xarray2], dim='time')

Extrapolation

Perform rover extrapolation

# Calling rover
motion_u, motion_v, reflectivity = rover(
    qpf_xarray, '201705241100', 6, 6
)

Out:

RUNNING 'rover' FOR EXTRAPOLATION.....
FINISHED EXTRAPOLATING BY "rover".....
################################################################################

Generating labelled radar colour images

Define the colour scale and the format of the reflectivity plots. Then generate using xarray.plot().

levels = [
    -32768,
    10, 15, 20, 24, 28, 32,
    34, 38, 41, 44, 47, 50,
    53, 56, 58, 60, 62
]
cmap = matplotlib.colors.ListedColormap([
    '#FFFFFF', '#08C5F5', '#0091F3', '#3898FF', '#008243', '#00A433',
    '#00D100', '#01F508', '#77FF00', '#E0D100', '#FFDC01', '#EEB200',
    '#F08100', '#F00101', '#E20200', '#B40466', '#ED02F0'
])

norm = BoundaryNorm(levels, ncolors=cmap.N, clip=True)

# Plot radar colour images
for i in range(1, len(reflectivity)):
    plt.figure(figsize=(15, 10))
    plt.axis("equal")
    reflectivity[i].plot(cmap=cmap, norm=norm)
    plt.savefig("rover-output-"+"{:02d}".format(i)+".png")

• 
• 
• 
• 
• 
• 

Generating radar reflectivity maps

Using the colour scale generated in the previous section, generate radar reflectivity maps with coastlines.

for j in range(1, len(reflectivity)):
    plt.figure(figsize=(15, 10))
    crs = area_def.to_cartopy_crs()
    ax = plt.axes(projection=crs)
    ax.coastlines(resolution='10m')
    plt.axis("equal")
    reflectivity[j].plot(cmap=cmap, norm=norm)
    plt.savefig("rover-output-map-"+"{:02d}".format(j)+".png")

• 
• 
• 
• 
• 
• 

Total running time of the script: ( 0 minutes 18.230 seconds)