Note
Click here to download the full example code
QPF (Manila)
This example demonstrates how to perform operational deterministic QPF up to three hours using raingauge data from Manila and radar data from Subic.
Definitions
import os
import numpy as np
import pandas as pd
from pyresample import utils
import xarray
import cartopy.feature as cfeature
import cartopy.crs as ccrs
import matplotlib
import imageio
import matplotlib.pyplot as plt
from matplotlib.colors import BoundaryNorm
from swirlspy.rad.uf_ph import read_uf_ph
from swirlspy.qpe.utils import dbz2rr, rr2rf, locate_file, timestamps_ending
from swirlspy.qpe.utils import multiple_acc
from swirlspy.obs.rain import Rain
from swirlspy.qpf.rover import rover
from swirlspy.qpf.sla import sla
from swirlspy.core.resample import grid_resample
plt.switch_backend('agg')
THIS_DIR = os.getcwd()
os.chdir(THIS_DIR)
start_time = pd.Timestamp.now()
Traceback (most recent call last):
File "/tmp/build/docs/swirlspy/swirlspy/examples/qpf_manila.py", line 30, in <module>
from swirlspy.qpf.rover import rover
File "/tmp/build/docs/swirlspy/swirlspy/qpf/rover.py", line 6, in <module>
from rover.rover import rover as rover_api
ImportError: libopencv_core.so.3.4: cannot open shared object file: No such file or directory
Initialising
This section demonstrates extracting radar reflectivity data.
Step 1: Define a basetime.
# Supply basetime
basetime = pd.Timestamp('20180811112000').floor('min')
Step 2: Using basetime, generate timestamps of desired radar files timestamps_ending() and locate files using locate_file().
# Obtain radar files
dir = THIS_DIR + '/../tests/samples/uf_ph/sub/'
located_files = []
radar_ts = timestamps_ending(
basetime,
duration=pd.Timedelta(minutes=60),
interval=pd.Timedelta(minutes=10)
)
for timestamp in radar_ts:
located_files.append(locate_file(dir, timestamp))
Step 3: Define the target grid as a pyresample AreaDefinition.
area_id = "epsg3123_240km"
description = ("A 240 m resolution rectangular grid "
"centred at Subic RADAR and extending to 240 km "
"in each direction")
proj_id = 'epsg3123'
projection = ('+proj=tmerc +lat_0=0 '
'+lon_0=121 +k=0.99995 +x_0=500000 '
'+y_0=0 +ellps=clrk66 +towgs84=-127.62,-67.24,'
'-47.04,-3.068,4.903,1.578,-1.06 +units=m '
'+no_defs')
x_size = 500
y_size = 500
area_extent = (191376.04113, 1399386.68659, 671376.04113, 1879386.68659)
area_def = utils.get_area_def(
area_id, description, proj_id, projection, x_size, y_size, area_extent
)
Step 4: Read data from radar files into xarray.DataArray using read_uf_ph().
reflectivity_list = [] # stores reflec from read_iris()
for filename in located_files:
reflec = read_uf_ph(
filename, area_def=area_def,
coord_label=['easting', 'northing'],
indicator='cappi', elevation=2
)
reflectivity_list.append(reflec)
Step 5: Assigning reflectivity xarrays at the last two timestamps to variables for use during ROVER QPF.
xarray1 = reflectivity_list[-2]
xarray2 = reflectivity_list[-1]
initialising_time = pd.Timestamp.now()
Running ROVER and Semi-Lagrangian Advection
Concatenate two reflectivity xarrays along time dimension.
Run ROVER, with the concatenated xarray as the input.
Perform Semi-Lagrangian Advection using the motion fields from rover.
# Combining the two reflectivity DataArrays
# the order of the coordinate keys is now ['y', 'x', 'time']
# as opposed to ['time', 'x', 'y']
reflec_concat = xarray.concat([xarray1, xarray2], dim='time')
# Rover
motion_u, motion_v = rover(
reflec_concat
)
rover_time = pd.Timestamp.now()
# Semi Lagrangian Advection
reflectivity = sla(
reflec_concat, motion_u, motion_v, steps=30
)
sla_time = pd.Timestamp.now()
Concatenating observed and forecasted reflectivities
Add forecasted reflectivity to reflectivity_list.
Concatenate observed and forecasted reflectivity xarray.DataArrays along the time dimension.
reflectivity_list.append(reflectivity[1:, ...])
reflectivity = xarray.concat(reflectivity_list, dim='time')
concat_time = pd.Timestamp.now()
Generating radar reflectivity maps
Define the color scale and format of the plots and plot using xarray.plot().
In this example, only hourly images will be plotted.
# Defining colour scale and format
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)
# Defining the crs
crs = area_def.to_cartopy_crs()
# Generating a timelist for every hour
timelist = [
(basetime + pd.Timedelta(minutes=60*i-10)) for i in range(4)
]
# Obtaining the slice of the xarray to be plotted
da_plot = reflectivity.sel(time=timelist)
# Defining motion quivers
qx = motion_u.coords['easting'].values[::5]
qy = motion_u.coords['northing'].values[::5]
qu = motion_u.values[::5, ::5]
qv = motion_v.values[::5, ::5]
# Defining coastlines
hires = cfeature.GSHHSFeature(
levels=[1],
scale='h',
edgecolor='k'
)
# Plotting
p = da_plot.plot(
col='time', col_wrap=2,
subplot_kws={'projection': crs},
cbar_kwargs={
'extend': 'max',
'ticks': levels[1:],
'format': '%.3g'
},
cmap=cmap,
norm=norm
)
for idx, ax in enumerate(p.axes.flat):
ax.quiver(qx, qy, qu, qv, pivot='mid', regrid_shape=20)
ax.add_feature(hires) # coastlines
ax.gridlines()
ax.set_title(
"Reflectivity\n"
f"Based @ {basetime.strftime('%H:%MH')}",
loc='left',
fontsize=9
)
ax.set_title(
''
)
ax.set_title(
f"{basetime.strftime('%Y-%m-%d')} \n"
f"Valid @ {timelist[idx].strftime('%H:%MH')} ",
loc='right',
fontsize=9
)
plt.savefig(
THIS_DIR +
f"/../tests/outputs/rover-output-map-mn.png",
dpi=300
)
radar_image_time = pd.Timestamp.now()
Accumulating hourly rainfall for 3-hour forecast
Hourly accumulated rainfall is calculated every 30 minutes, the first endtime is the basetime i.e. T+0min.
Convert reflectivity in dBZ to rainrates in mm/h with dbz2rr().
Changing time coordinates of xarray from start time to endtime.
Convert rainrates to rainfalls in 10 mins with rr2rf().
Accumulate hourly rainfall every 30 minutes using multiple_acc().
# Convert reflectivity to rainrates
rainrates = dbz2rr(reflectivity, a=300, b=1.4)
# Converting the coordinates of xarray from start to endtime
rainrates_endtime = rainrates.copy()
rainrates_endtime.coords['time'] = \
[
pd.Timestamp(t) + pd.Timedelta(minutes=10)
for t in rainrates_endtime.coords['time'].values
]
# Convert rainrates to accumulated rainfalls every 10 minutes with rr2rf().
rainfalls = rr2rf(rainrates_endtime, scan_duration=10)
# Accumulate hourly rainfall every 30 minutes
acc_rf = multiple_acc(
rainfalls, basetime, basetime+pd.Timedelta(hours=3)
)
acc_time = pd.Timestamp.now()
Plotting rainfall maps
Define the colour scheme and format and plot using xarray.plot().
In this example, only hourly images will be plotted.
# Defining the colour scheme
levels = [
0, 0.5, 2, 5, 10, 20,
30, 40, 50, 70, 100, 150,
200, 300, 400, 500, 600, 700
]
cmap = matplotlib.colors.ListedColormap([
'#ffffff', '#9bf7f7', '#00ffff', '#00d5cc', '#00bd3d', '#2fd646',
'#9de843', '#ffdd41', '#ffac33', '#ff621e', '#d23211', '#9d0063',
'#e300ae', '#ff00ce', '#ff57da', '#ff8de6', '#ffe4fd'
])
norm = BoundaryNorm(levels, ncolors=cmap.N, clip=True)
# Defining projection
crs = area_def.to_cartopy_crs()
# Defining zoom extent
r = 64000
proj_site = xarray1.proj_site
zoom = (
proj_site[0]-r, proj_site[0]+r, proj_site[1]-r, proj_site[1]+r
) # (x0, x1, y0, y1)
# Defining times for plotting
timelist = [basetime + pd.Timedelta(hours=i) for i in range(4)]
# Obtaining xarray slice to be plotted
da_plot = acc_rf.sel(
easting=slice(zoom[0], zoom[1]),
northing=slice(zoom[3], zoom[2]),
time=timelist
)
# Plotting
p = da_plot.plot(
col='time', col_wrap=2,
subplot_kws={'projection': crs},
cbar_kwargs={
'extend': 'max',
'ticks': levels,
'format': '%.3g'
},
cmap=cmap,
norm=norm
)
for idx, ax in enumerate(p.axes.flat):
ax.add_feature(hires) # using GSHHS coastlines defined previously
ax.gridlines()
ax.set_xlim(zoom[0], zoom[1])
ax.set_ylim(zoom[2], zoom[3])
ax.set_title(
"Past Hour Rainfall\n"
f"Based @ {basetime.strftime('%H:%MH')}",
loc='left',
fontsize=8
)
ax.set_title(
''
)
ax.set_title(
f"{basetime.strftime('%Y-%m-%d')} \n"
f"Valid @ {timelist[idx].strftime('%H:%MH')} ",
loc='right',
fontsize=8
)
plt.savefig(
THIS_DIR +
f"/../tests/outputs/rainfall_mn.png",
dpi=300
)
rf_image_time = pd.Timestamp.now()
Extract the rainfall values at a specified location
In this example, the rainfall values at the location is assumed to be the same as the nearest gridpoint.
Read information regarding the rain gauge stations into a pandas.DataFrame.
Extract the rainfall values at the nearest gridpoint to location for given timesteps (in this example, 30 minute intervals).
Store rainfall values over time in a pandas.DataFrame.
Plot the time series of rainfall at different stations.
# Getting rain gauge station coordinates
df = pd.read_csv(
os.path.join(THIS_DIR, "../tests/samples/manila_rg_list.csv"),
delim_whitespace=True,
usecols=[0, 3, 4]
)
# Extract rainfall values at gridpoint closest to the
# location specified for given timesteps and storing it
# in pandas.DataFrame.
rf_time = []
for time in acc_rf.coords['time'].values:
rf = []
for index, row in df.iterrows():
rf.append(acc_rf.sel(
time=time, northing=row[2],
easting=row[1],
method='nearest'
).values)
rf_time.append(rf)
rf_time = np.array(rf_time)
station_rf = pd.DataFrame(
data=rf_time,
columns=df.iloc[:, 0],
index=pd.Index(
acc_rf.coords['time'].values,
name='time'
)
)
print(station_rf)
loc_stn = \
['BAA', 'BUM', 'PAF', 'QUL', 'ZAP', 'ZAA']
loc_stn_drop = [
stn for stn in station_rf.columns.to_list()
if stn not in loc_stn
]
df_loc = station_rf.drop(loc_stn_drop, axis=1)
print(df_loc)
# Plotting time series graph for selected stations
ax = df_loc.plot(title="Time Series of Hourly Accumulated Rainfall",
grid=True)
ax.set_ylabel("Hourly Accumulated Rainfall (mm)")
plt.savefig(THIS_DIR+"/../tests/outputs/qpf_time_series.png")
extract_time = pd.Timestamp.now()
Checking run time of each component
print(f"Start time: {start_time}")
print(f"Initialising time: {initialising_time}")
print(f"Rover time: {rover_time}")
print(f"SLA time: {sla_time}")
print(f"Concatenating time: {concat_time}")
print(f"Plotting radar image time: {radar_image_time}")
print(f"Accumulating rainfall time: {acc_time}")
print(f"Plotting rainfall map time: {rf_image_time}")
print(f"Extracting and plotting time series time: {extract_time}")
print(f"Time to initialise: {initialising_time-start_time}")
print(f"Time to run rover: {rover_time-initialising_time}")
print(f"Time to perform SLA: {sla_time-rover_time}")
print(f"Time to concatenate xarrays: {concat_time - sla_time}")
print(f"Time to plot radar image: {radar_image_time - concat_time}")
print(f"Time to accumulate rainfall: {acc_time - radar_image_time}")
print(f"Time to plot rainfall maps: {rf_image_time-acc_time}")
print(f"Time to extract and plot time series: {extract_time-rf_image_time}")
Total running time of the script: ( 0 minutes 0.003 seconds)