Temporal interpolation using nowcasting
Nowcast-based temporal interpolation was developed to solve the three following problems simultaneously:
Data is often needed at times other than the six minutes at which radar mosaics are available. For data assimilation purposes, radar data is required at time resolutions between 1 and 7.5 minutes, depending on the timestep of the model being used.
Prior to assimilation, radar data must be smoothed for its spatial resolution to match the effective resolution of the simulated precipitation. This avoids assimilating small scales that the model cannot represent.
The source data contains many gaps; US radars leave some when they use scanning routines (VCPs) longer than six minutes, and Canadian data is sometimes missing because of network delays.
The animation below shows the result of the interpolation for a case over Minnesota and Wisconsin on 29 August 2022. The top row shows the source data, available every six minutes, together with its quality index. The bottom row shows the interpolated precipitation rate and its quality index, available every minute. Note how the interpolation bridges the gap when source data is unavailable.
The interpolation is performed by the obs_process module. The basic idea is that the temporally interpolated data at intermediate timesteps is estimated as a quality-weighted aggregate of advected neighbors. By selecting \(\delta_t^{\text{max}}\) such that 4 or 5 neighbors contribute to each aggregate, data gaps go mostly unnoticed.
Motion vectors are first estimated from consecutive source mosaics using the Lucas-Kanade optical flow of pysteps. Every source mosaic that is closer in time than \(\delta_t^{\text{max}}\) then contributes to the output; each one is advected to the desired time, forward or backward, with a semi-Lagrangian scheme. The weight of each contribution decreases linearly with the duration of the advection and reaches zero at \(\delta_t^{\text{max}}\). Because this weight multiplies the quality index of the advected data, an estimate that had to be advected over a long time interval is given less importance than one that comes from a nearby time. The precipitation rate is the average of the advected estimates weighted by these reduced quality indices.
Interpolate a batch of observations
This section demonstrates the processing of a batch of observations in one call to the obs_process function. The results are then displayed in the form of a short animation.
Let’s start with the required imports and directory setup:
import os
import datetime
import subprocess
import glob
import shutil
import copy
import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
import dask
from dask.distributed import Client, LocalCluster
import domutils.legs as legs
import domutils.geo_tools as geo_tools
import domutils.radar_tools as radar_tools
import domcmc.fst_tools as fst_tools
import domutils._py_tools as py_tools
#setting up directories
test_data_dir = setup_test_paths['test_data_dir']
test_results_dir = setup_test_paths['test_results_dir']
generated_files_dir = os.path.join(test_results_dir, 'generated_files', 'test_radar_time_interpolation')
generated_figure_dir = os.path.join(test_results_dir, 'generated_figures', 'test_radar_time_interpolation')
reference_figure_dir = os.path.join(test_data_dir, 'reference_figures', 'test_radar_time_interpolation')
py_tools.parallel_mkdir(generated_files_dir)
py_tools.parallel_mkdir(generated_figure_dir)
obs_process is a python script callable from the shell such as:
#process data with time interpolation python -m domutils.radar_tools.obs_process \ --input_t0 202206150800 \ --input_tf 202206160000 \ --input_dt 10 \ --output_t0 202206150900 \ --output_tf 202206160000 \ --output_dt 1 \ --t_interp_method 'nowcast' \ ...
However, for this example we will be running directly from Python with arguments provided by the attributes of a simple object.
# a class that mimics the output of argparse
class ArgsClass():
def __init__(self, setup_test_paths):
test_data_dir = setup_test_paths['test_data_dir']
test_results_dir = setup_test_paths['test_results_dir']
generated_files_dir = os.path.join(test_results_dir, 'generated_files', 'test_radar_time_interpolation')
generated_figure_dir = os.path.join(test_results_dir, 'generated_figures', 'test_radar_time_interpolation')
self.input_t0 = '202208290230'
self.input_tf = '202208290530'
self.input_dt = '6M'
self.output_t0 = '202208290330'
self.output_tf = '202208290430'
self.output_dt = '1M'
self.interp_max_dt = '13M'
self.complete_dataset = 'False'
self.t_interp_method = 'nowcast'
self.input_data_dir = os.path.join(test_data_dir, 'odimh5_radar_composites')
self.input_file_struc = '%Y/qcomp_%Y%m%d%H%M.h5'
self.h5_latlon_file = os.path.join(test_data_dir, 'radar_continental_2.5km_2882x2032.pickle')
self.sample_pr_file = os.path.join(test_data_dir, 'hrdps_5p1_prp0.fst')
self.ncores = 40 # use as many cpus as you have on your system
self.preproc_median_filt = '3'
self.preproc_smooth_radius = '4'
self.nowcast_median_filt = '3'
self.output_dir = os.path.join(generated_files_dir, 'obs_process_t_interp')
self.output_file_format = 'fst'
self.output_file_struc = '%Y%m%d%H%M.fst'
self.figure_format = 'svg'
self.log_level = 'WARNING'
We now instantiate this object for the case being demonstrated and define the domain and the times that will be used for the figures.
import cartopy.crs as ccrs
# a new argument object is needed for each case since obs_process
# modifies the arguments it receives
args = ArgsClass(setup_test_paths)
# times of the source data that will be displayed in the animation
anim_t0 = datetime.datetime(2022,8,29,3,42)
anim_source_deltat = np.arange(0, 37, 6, dtype=int) # minutes
# domain of the figures; 300x300 km over Minnesota/Wisconsin
pole_latitude=35.7
pole_longitude=65.5
lat_0 = 46.7
delta_lat = 3.14*.5
lon_0 = 267.3
delta_lon = 4.17*.5
map_extent=[lon_0-delta_lon, lon_0+delta_lon, lat_0-delta_lat, lat_0+delta_lat]
proj_aea = ccrs.RotatedPole(pole_latitude=pole_latitude, pole_longitude=pole_longitude)
The processing of observations and time interpolation is done in one simple function call.
radar_tools.obs_process(args)
To make an animation showing the time-interpolated data, we first define a function for plotting each individual panel.
def plot_panel(data,
fig, ax_pos, title,
proj_aea,
proj_obj, colormap,
plot_palette=None,
pal_units=None,
show_artefacts=False):
'''function that plots individual panels
'''
import matplotlib.pyplot as plt
import cartopy.feature as cfeature
ax = fig.add_axes(ax_pos, projection=proj_aea)
ax.set_extent(proj_obj.rotated_extent, crs=proj_aea)
dum = ax.annotate(title, size=32,
xy=(.022, .85), xycoords='axes fraction',
bbox=dict(boxstyle="round", fc='white', ec='white'))
# projection from data space to image space
projected_data = proj_obj.project_data(data)
# plot data & palette
colormap.plot_data(ax=ax, data=projected_data,
palette=plot_palette,
pal_units=pal_units, pal_format='{:5.1f}',
equal_legs=True)
# add political boundaries
ax.add_feature(cfeature.STATES.with_scale('10m'), linewidth=0.5, edgecolor='0.2')
# show artefacts in accumulation plots
if show_artefacts:
ax2 = fig.add_axes(ax_pos)
ax2.set_xlim((0.,1.))
ax2.set_ylim((0.,1.))
ax2.patch.set_alpha(0.0)
ax2.set_axis_off()
xpos = np.linspace(0.24, 0.40, 5)
ypos = np.linspace(0.75, 0.85, 5)
for x0, y0, dx in [(xx,yy,.1) for xx, yy in zip(xpos, ypos)]:
ax2.arrow(x0, y0, dx, -.03,
width=0.015, facecolor='red', edgecolor='black',
head_width=3*0.01, linewidth=2.)
def figure_for_timestep(src_delta_min, interp_delta_min, t0,
proj_aea,
input_proj_obj, output_proj_obj,
pr_colormap, qi_colormap,
generated_figure_dir,
args, fig_w, fig_h, sp_w, sp_h, rec_w, rec_h, sp_m ):
'''make the figure for every output timestep
'''
import subprocess
import datetime
import matplotlib as mpl
import matplotlib.pyplot as plt
import domutils.radar_tools as radar_tools
source_t_offset = datetime.timedelta(seconds=src_delta_min * 60.0)
source_valid_time = t0 + source_t_offset
interpolated_t_offset = datetime.timedelta(seconds=interp_delta_min * 60.0)
interpolated_valid_time = source_valid_time + interpolated_t_offset
# matplotlib global settings
dpi = 400
mpl.rcParams.update({
'font.family': 'Latin Modern Roman',
'font.size': 32,
'axes.titlesize': 32,
'axes.labelsize': 32,
'xtick.labelsize': 30,
'ytick.labelsize': 30,
'legend.fontsize': 30,
'figure.dpi': dpi,
'savefig.dpi': dpi,
})
# instantiate figure
fig = plt.figure(figsize=(fig_w,fig_h))
# source data on original grid
dat_dict = radar_tools.get_instantaneous(desired_quantity='precip_rate',
valid_date=source_valid_time,
data_path=args.input_data_dir,
data_recipe=args.input_file_struc)
x0 = sp_w
y0 = 2.*sp_h + rec_h
ax_pos = [x0, y0, rec_w, rec_h]
title = f'Source precip rate \n @ t0+{src_delta_min}minutes'
plot_panel(dat_dict['precip_rate'],
fig, ax_pos, title,
proj_aea,
input_proj_obj, pr_colormap,
plot_palette='right',
pal_units='[mm/h]')
# source quality index
x0 = sp_w + rec_w + sp_m
y0 = 2.*sp_h + rec_h
ax_pos = [x0, y0, rec_w, rec_h]
title = f'Source quality index \n @ t0+{src_delta_min}minutes'
plot_panel(dat_dict['total_quality_index'],
fig, ax_pos, title,
proj_aea,
input_proj_obj, qi_colormap,
plot_palette='right',
pal_units='[unitless]')
# Time interpolated data
dat_dict = radar_tools.get_instantaneous(desired_quantity='precip_rate',
valid_date=interpolated_valid_time,
data_path=args.output_dir,
data_recipe=args.output_file_struc)
x0 = sp_w
y0 = sp_h
ax_pos = [x0, y0, rec_w, rec_h]
title = f'Interpolated precip rate \n @ t0+{src_delta_min+interp_delta_min}minutes'
plot_panel(dat_dict['precip_rate'],
fig, ax_pos, title,
proj_aea,
output_proj_obj, pr_colormap,
plot_palette='right',
pal_units='[mm/h]')
# quality index
x0 = sp_w + rec_w + sp_m
y0 = sp_h
ax_pos = [x0, y0, rec_w, rec_h]
title = f' Interpolated quality index\n @ t0+{src_delta_min+interp_delta_min}minutes'
plot_panel(dat_dict['total_quality_index'],
fig, ax_pos, title,
proj_aea,
output_proj_obj, qi_colormap,
plot_palette='right',
pal_units='[unitless]')
# save output
date_prefix = interpolated_valid_time.strftime('%Y%m%d%H%M')
fig_name = os.path.join(generated_figure_dir, f'{date_prefix}_frame.png')
plt.savefig(fig_name)
#fig_name_svg = os.path.join(generated_figure_dir, f'{date_prefix}_frame.svg')
#plt.savefig(fig_name_svg)
plt.close(fig)
print(f'done with {fig_name}')
# use "convert" to make a gif out of the png
cmd = ['convert', fig_name, '-geometry', '15%', '-quantize', 'transparent', '-dither', 'FloydSteinberg', '-colors', '256', fig_name.replace('png', 'gif')]
process = subprocess.Popen(cmd, stdout=subprocess.PIPE)
output, error = process.communicate()
# we don't need the original png anymore
os.remove(fig_name)
return source_valid_time, interpolated_valid_time
then we setup the general characteristics of the figure being generated. See Legs Tutorial for information on the definition of color mapping objects.
dpi = 400
mpl.rcParams.update({
'font.family': 'Latin Modern Roman',
'font.size': 32,
'axes.titlesize': 32,
'axes.labelsize': 32,
'xtick.labelsize': 30,
'ytick.labelsize': 30,
'legend.fontsize': 30,
'figure.dpi': dpi,
'savefig.dpi': dpi,
})
#pixel density of each panel
ratio = 1.
hpix = 1200. #number of horizontal pixels
vpix = ratio*hpix #number of vertical pixels
img_res = (int(hpix),int(vpix))
#size of image to plot
fig_w = 19. #size of figure
fig_h = 15.7 #size of figure
rec_w = 7./fig_w #size of axes
rec_h = ratio*(rec_w*fig_w)/fig_h #size of axes
sp_w = .5/fig_w #space between panel and border
sp_m = 2.2/fig_w #space between panels
sp_h = .5/fig_h #space between panels
# color mapping object
range_arr = [.1,1.,5.,10.,25.,50.,100.]
missing = -9999.
# colormap object for precip rates
pr_colormap = legs.PalObj(range_arr=range_arr,
n_col=6,
over_high='extend', under_low='white',
excep_val=missing,
excep_col='grey_200')
# colormap for QI index
pastel = [ [[255,190,187],[230,104, 96]], #pale/dark red
[[255,185,255],[147, 78,172]], #pale/dark purple
[[255,227,215],[205,144, 73]], #pale/dark brown
[[210,235,255],[ 58,134,237]], #pale/dark blue
[[223,255,232],[ 61,189, 63]] ] #pale/dark green
qi_colormap = legs.PalObj(range_arr=[0., 1.],
dark_pos='high',
color_arr=pastel,
excep_val=[missing,0.],
excep_col=['grey_220','white'])
# get lat/lon of input data from one of the h5 files
dum_h5_file = os.path.join(test_data_dir, 'odimh5_radar_composites', '2022/qcomp_202205212000.h5')
input_ll = radar_tools.read_h5_composite(dum_h5_file, latlon=True)
input_lats = input_ll['latitudes']
input_lons = input_ll['longitudes']
# get lat/lon of output data
output_ll = fst_tools.get_data(args.sample_pr_file, var_name='PR', latlon=True)
output_lats = output_ll['lat']
output_lons = output_ll['lon']
# instantiate projection object for input data
input_proj_obj = geo_tools.ProjInds(src_lon=input_lons, src_lat=input_lats,
extent=map_extent, dest_crs=proj_aea, image_res=img_res)
# instantiate projection object for output data
output_proj_obj = geo_tools.ProjInds(src_lon=output_lons, src_lat=output_lats,
extent=map_extent, dest_crs=proj_aea, image_res=img_res)
Individual frames of the animation are made serially or in parallel.
interpolated_deltat = np.arange(6) # minutes
# making the figures only requires the file paths found in args
# obs_process attached large arrays and projection objects to this object;
# they are removed here so that only light-weight data is sent to the dask workers
plot_args = copy.copy(args)
for attr in ('dask_client', 'proj_obj', 'out_lats', 'out_lons'):
setattr(plot_args, attr, None)
serial=False
if serial:
for src_delta_min in anim_source_deltat:
for interp_delta_min in interpolated_deltat:
figure_for_timestep(src_delta_min, interp_delta_min, anim_t0,
proj_aea,
input_proj_obj, output_proj_obj,
pr_colormap, qi_colormap,
generated_figure_dir,
plot_args, fig_w, fig_h, sp_w, sp_h, rec_w, rec_h, sp_m)
else:
tasks = [dask.delayed(figure_for_timestep)(src_delta_min, interp_delta_min, anim_t0,
proj_aea,
input_proj_obj, output_proj_obj,
pr_colormap, qi_colormap,
generated_figure_dir,
plot_args, fig_w, fig_h, sp_w, sp_h, rec_w, rec_h, sp_m)
for src_delta_min in anim_source_deltat
for interp_delta_min in interpolated_deltat]
with LocalCluster(processes=True, n_workers=10, threads_per_worker=1) as cluster, Client(cluster) as client:
results = dask.compute(*tasks) # parallel execution
Finally, an animated gif is constructed from the frames we just made,
date_prefix = anim_t0.strftime('%Y%m')
movie_name = os.path.join(generated_figure_dir, f'{date_prefix}_movie.gif')
gif_list = sorted(glob.glob(os.path.join(generated_figure_dir,f'{date_prefix}*frame.gif')))
cmd = ['convert', '-loop', '0', '-delay', '30']+gif_list+[movie_name]
process = subprocess.Popen(cmd, stdout=subprocess.PIPE)
output, error = process.communicate()
The animation shown at the top of this page was obtained with the code above. Running the same code with different dates and a different domain gives the animation below, for the derecho that crossed southern Quebec on 21 May 2022.
Accumulations from time interpolated data
Using nowcasting for time interpolation can be advantageous when computing accumulations from source data available at discrete times. In the example below, we compare accumulations obtained from the source data to accumulations obtained from the time interpolated data.
duration = 30 # minutes
end_date = anim_t0 + datetime.timedelta(minutes=30)
# the next figure has only one row of panels
# the panels keep the same size, only the height of the figure changes
fig_h = 8. #size of figure
rec_h = ratio*(rec_w*fig_w)/fig_h #size of axes
sp_h = .5/fig_h #space between panel and border
# instantiate figure
fig = plt.figure(figsize=(fig_w,fig_h))
# make accumulation from source data
dat_dict = radar_tools.get_accumulation(end_date=end_date,
duration=duration,
input_dt=6., # minutes
data_path=args.input_data_dir,
data_recipe=args.input_file_struc)
x0 = 2.*sp_w + rec_w
y0 = sp_h
ax_pos = [x0, y0, rec_w, rec_h]
title = 'Accumulation from \n source data'
plot_panel(dat_dict['accumulation'],
fig, ax_pos, title,
proj_aea,
input_proj_obj, pr_colormap,
plot_palette='right',
pal_units='mm',
show_artefacts=True)
# make accumulation from time interpolated data
dat_dict = radar_tools.get_accumulation(end_date=end_date,
duration=duration,
input_dt=1., # minutes
data_path=args.output_dir,
data_recipe=args.output_file_struc)
x0 = sp_w
y0 = sp_h
ax_pos = [x0, y0, rec_w, rec_h]
title = 'Accumulation from \n time interpolated data'
plot_panel(dat_dict['accumulation'],
fig, ax_pos, title,
proj_aea,
output_proj_obj, pr_colormap)
# save output
fig_name = os.path.join(generated_figure_dir, f'{date_prefix}_accumulation.svg')
plt.savefig(fig_name)
plt.close(fig)
The figure below shows the 30 minutes precipitation accumulation ending at 21:42 UTC on 21 May 2022, computed from:
the source data, available every six minutes, on the right
the time interpolated data, available every minute, on the left
In the panel on the right, the red arrows indicate artefacts that originate from the coarse time resolution of the source data compared to the speed at which the bow echo propagates. The accumulation on the left does not display these displacement artefacts.
