This piece of code attempts to focus on features which can be tracked over a long time (i.e. 30 minutes) and which rather move along a "straight line". The motivation is to use only those motion patterns for forecasting/nowcasting which have a minimum level of persistence (and ignore individual features that show a rather erratic behaviour).
In [ ]:
%pylab
import wradlib
import trackpy as tp
import numpy as np
import os
from scipy import stats
from matplotlib import animation
from scipy.ndimage import zoom
Data is from the German Weather Service: the so called RY product represents rainfall intensity composite for the whole of Germany in 5 minute intervals.
Spatial resolution: 1 x 1 km; spatial extent: 900 x 900 km.
Information required from user
datadir where you store the RY data (unpack the ry.rar archive there).dtimes lines.downsizeby) in order to avoid memory problems
In [ ]:
# Set data directory
datadir = r"E:\src\python\nowcast\trackpy\ry"
# Original grid dimensions
nx = 900
ny = 900
# pixel size (in meters)
dx = 1000.
dy = 1000.
# Downsize by factor "downsizeby"
# downsizeby = 1 will leave the dimensions unchanged,
# but for a 900x900 km grid, downsizing might be
# required in order to avoid MemoryError
downsizeby = 2
# interval between observations (in seconds)
interval = 300
# Set time window
##dtimes = wradlib.util.from_to("2008-06-02 17:00:00", "2008-06-02 19:00:00", 300)
##dtimes = wradlib.util.from_to("2015-04-26 17:00:00", "2015-04-26 19:00:00", 300)
##dtimes = wradlib.util.from_to("2015-03-29 17:00:00", "2015-03-29 19:00:00", 300)
dtimes = wradlib.util.from_to("2015-04-01 17:00:00", "2015-04-01 19:00:00", interval)
In [ ]:
# Compute grid dimensions and grid coordinates after resampling
dx2, dy2 = dx*downsizeby, dy*downsizeby
nx2, ny2 = nx/downsizeby, ny/downsizeby
X2, Y2 = np.meshgrid( np.arange(0,nx2*dx2, dx2), np.arange(0,ny2*dy2, dy2) )
# Define container
frames = np.zeros( (len(dtimes), nx2, ny2 ) )
# Move to data directory
os.chdir(datadir)
# Read the data, convert to dBZ, and downsize
for i, dtime in enumerate(dtimes):
fname = dtime.strftime("raa01-ry_10000-%y%m%d%H%M-dwd---bin")
frames[i] = zoom( wradlib.io.read_RADOLAN_composite(fname, missing=0)[0], 1./downsizeby, order=1)
frames[i] = wradlib.trafo.decibel( wradlib.zr.r2z(frames[i]) )
frames[i][frames[i]<0] = 0
tp.batchtp.link_dfThe idea is to track only features that represent "persistent motion patterns". Please note that those features/tracks will be filtered out which have not been tracked over at least 6 times steps (i.e. 6 * 5 minutes = 30 minutes). This is due to a specific requirement: I want to detect the rather persistent motion.
In [ ]:
# FEATURE DETECTION
diameter = 11
feats = tp.batch(frames, diameter=diameter, minmass=500., invert=False)
In [ ]:
# FEATURE TRACKING
# Assumed maximum diplacement velocity: 150 km/h
# (see http://www.meteo.fr/cic/wsn05/resumes_longs/2.14-73.pdf)
# Corresponds to about 12 km in 5 minutes
# This is required for the search radius.
# BEWARE: You need to provide this as a number of pixels, so you need to consider downsizeby
# We set memory=0 because we do not expect features to disappear and reappear again.
tracks = tp.link_df(feats, int(12/downsizeby), memory=0)
In [ ]:
# FILTERING TRACKS: Only persistent tracks survive (i.e. which persist over min_tsteps)
min_tsteps = 6
tracks1 = tp.filter_stubs(tracks, min_tsteps)
# Compare the number of particles in the unfiltered and filtered data.
print 'Features before:', tracks['particle'].nunique()
print 'Features after:', tracks1['particle'].nunique()
In [ ]:
# Quickview at the tracks
ax = plt.subplot(111,aspect="equal")
tp.plot_traj(tracks1, ax=ax)
plt.xlim(0, nx2)
plt.ylim(0, ny2)
In [ ]:
# Analyze trajectories: only "straight" trajectories survive
my = tracks
max_p = 0.001
ids, dvx, dvy, xmean, ymean = np.array([]), np.array([]), np.array([]), np.array([]), np.array([])
for id in np.unique(my.particle).astype("i4"):
sub = my[my.particle==id]
slope, intercept, r_value, p_value, std_err = stats.linregress(sub.x,sub.y)
if p_value < max_p:
# Compute velocity components in m/s
dvx = np.append(dvx, np.mean(np.diff(sub.x)) * dx2 / interval )
dvy = np.append(dvy, np.mean(np.diff(sub.y)) * dy2 / interval )
xmean = np.append(xmean, np.mean(sub.x)*dx2 )
ymean = np.append(ymean, np.mean(sub.y)*dy2 )
ids = np.append(ids, id)
# Plot remaining arrows over average reflectivity
ax = plt.subplot(111,aspect="equal")
#plt.imshow(np.mean(frames, axis=0), origin='lower')
plt.pcolormesh(X2, Y2, np.mean(frames, axis=0))
#tp.plot_traj(my, ax=ax)
plt.quiver(xmean, ymean, dvx, dvy, pivot='middle', headwidth=4, headlength=6, color='red')
for i, id in enumerate(ids):
plt.text(xmean[i], ymean[i], "%d"%id, color="white" )
In [ ]:
# Animation
plotthis = tracks1
_plot_style = dict(markersize=diameter, markeredgewidth=2,
markerfacecolor='none', markeredgecolor='r',
marker='o', linestyle='none')
_imshow_style = dict(interpolation='none',
cmap=plt.cm.spectral)
fig = plt.figure(figsize=(16,8))
ax1 = fig.add_subplot(121, aspect="equal")
ax2 = fig.add_subplot(122, aspect="equal")
im1 = ax1.imshow(frames[0], **_imshow_style)
ax1.set_xlim(-0.5, frames[0].shape[1] - 0.5)
ax1.set_ylim(-0.5, frames[0].shape[0] - 0.5)
title1 = ax1.set_title(dtimes[0].isoformat())
#ax1.invert_yaxis()
im2 = ax2.imshow(frames[0], **_imshow_style)
ax2.set_xlim(-0.5, frames[0].shape[1] - 0.5)
ax2.set_ylim(-0.5, frames[0].shape[0] - 0.5)
title2 = ax2.set_title(dtimes[0].isoformat())
#ax2.invert_yaxis()
bubbles, = ax2.plot(plotthis[plotthis["frame"]==0]["x"], plotthis[plotthis["frame"]==0]["y"], **_plot_style)
def animate(i):
im1.set_data(frames[i])
im2.set_data(frames[i])
bubbles.set_data(plotthis[plotthis["frame"]==i]["x"], plotthis[plotthis["frame"]==i]["y"])
ax1.set_title(dtimes[i].isoformat())
ax2.set_title(dtimes[i].isoformat())
return im1, im2, bubbles, title1, title2
#Init only required for blitting to give a clean slate.
def init():
im1.set_data(frames[0])
im2.set_data(frames[0])
bubbles.set_data(plotthis[plotthis["frame"]==0]["x"], plotthis[plotthis["frame"]==0]["y"])
ax1.set_title(dtimes[0].isoformat())
ax2.set_title(dtimes[0].isoformat())
return im1, im2, bubbles, title1, title2
ani = animation.FuncAnimation(fig, animate, frames=np.arange(24), interval=100, blit=False)#, init_func=init,
# interval=100, blit=False)#, repeat_delay=100)
In [ ]:
srcxy = np.vstack((xmean, ymean)).transpose()
trgxy = np.vstack( (X2.ravel(), Y2.ravel()) ).transpose()
trgx = trgxy[:,0].reshape((nx2,ny2))
trgy = trgxy[:,1].reshape((nx2,ny2))
ip_kwargs = {"nnearest":4, "p":1.}
dvx_ip = wradlib.ipol.interpolate(srcxy, trgxy, dvx, wradlib.ipol.Idw, **ip_kwargs).reshape((nx2,ny2))
dvy_ip = wradlib.ipol.interpolate(srcxy, trgxy, dvy, wradlib.ipol.Idw, **ip_kwargs).reshape((nx2,ny2))
In [ ]:
# Define how many arrows quiver will plot
ith = 15
fig = plt.figure(figsize=(10,10))
ax = fig.add_subplot(111, aspect="equal")
# Average reflectivity over analysis period (just for orientation)
cl = plt.pcolormesh(X2, Y2, np.mean(frames, axis=0))
# Interpolated displacement vector field
plt.quiver(trgx[0:nx2:ith,0:ny2:ith], trgy[0:nx2:ith,0:ny2:ith], dvx_ip[0:nx2:ith,0:ny2:ith], dvy_ip[0:nx2:ith,0:ny2:ith],
pivot='middle', headwidth=4, headlength=6, color='red')
# Displacement vectors that were used for interpolation
plt.quiver(xmean, ymean, dvx, dvy, pivot='middle', headwidth=4, headlength=6, color='white')
# Other decorators
plt.grid(color="white")
plt.xlabel("kilometers")
plt.ylabel("kilometers")
cb = plt.colorbar(cl, shrink=0.75)
cb.set_label("dBZ")
For field advection we use FiPy, a finite volume PDE solver. It provides stable methods to solve the 2d advection problem. Just consider rainfall intensity as a quantity that is advected with a 2D velocity field.
Our worst enemies here: numerical diffusion and performance. Solving the advection PDE certainly is the performance bottleneck.
Memory is a problem, too: at least on my 32-bit Python, I need to reduce the spatial resolution by a factor of 2 in order to get this running.
In [ ]:
import fipy
In [ ]:
# Some information to numerically solve the initial-boundary vaue problem
# Set numerical solver
MySolver = fipy.LinearLUSolver #DefaultAsymmetricSolver
# Initial condition
init_step = 12
# Boundary condition: always zero (dBZ or rainfall intensity...doesn't matter)
boundary = 0. # dBZ
# Forecast time settings
lead_time = 3600 # in seconds
# CFL number is defined as C = u * dt/dx and typically C_max = 1, so u * dt/dx < C_max yields dt < dx/u
dt = 30 # in seconds
numsteps = int(lead_time / dt)
In [ ]:
# Define the 2D mesh
mesh = fipy.Grid2D(dx=dx2, dy=dy2, nx=nx2, ny=ny2)
# Set the 2D scalar variable (rainfall intensity or reflectivity)
var = fipy.CellVariable(name = "solution variable", mesh = mesh, value = frames[init_step].copy().ravel())
# Set boundary conditions
var.constrain(boundary, mesh.facesLeft)
var.constrain(boundary, mesh.facesRight)
var.constrain(boundary, mesh.facesBottom)
var.constrain(boundary, mesh.facesTop)
In [ ]:
# Set displacement vector field on the same mesh (i.e. convection coefficient - assumed persistent)
convCoeff = fipy.CellVariable(mesh=mesh, value=[dvx_ip.ravel(), dvy_ip.ravel()])
In [ ]:
# Set the advection equation
# the VanLeerConvectionTerm appears to minimize numerical diffusion compared to the other methods)
eq = fipy.TransientTerm() + fipy.VanLeerConvectionTerm(coeff = convCoeff)
##eq = fipy.TransientTerm() + fipy.FirstOrderAdvectionTerm(coeff=convCoeff)
##eq = fipy.TransientTerm() + fipy.ExplicitUpwindConvectionTerm(coeff=convCoeff)
##eq = fipy.TransientTerm() + fipy.ExponentialConvectionTerm(coeff=convCoeff)
In [ ]:
# Solve the advection equation (SLOW!)
# This is our results container
# index 0 corresponds to our initial condition, so index "i" is the "ith" forecast step
forecast = np.zeros((numsteps+1,nx2,ny2))
forecast[0] = frames[init_step].copy()
for i in range(numsteps):
eq.solve(var=var, dt=dt)#, solver=MySolver(tolerance=1.e-15, iterations=1000), dt=dt)
forecast[i+1] = var.value.reshape((nx2,ny2))
In [ ]:
# Animate forecast and "reality"
_plot_style = dict(markersize=diameter, markeredgewidth=2,
markerfacecolor='none', markeredgecolor='r',
marker='o', linestyle='none')
_pcm_style = dict(cmap=plt.cm.spectral, vmin=0., vmax=30.)
# Prepare canvas
fig = plt.figure(figsize=(16,12))
ax1 = plt.subplot2grid((3,4), (0,0), colspan=2, rowspan=2, aspect="equal")
ax2 = plt.subplot2grid((3,4), (0,2), colspan=2, rowspan=2, aspect="equal")
ax3 = plt.subplot2grid((3,4), (2,0), colspan=4, rowspan=1)
# Initialize map of observed rainfall
im1 = ax1.pcolormesh(X2, Y2, frames[init_step], **_pcm_style)
ax1.grid(color="white")
title1 = ax1.set_title("Observation")
tstamp1 = ax1.text(50000., 850000., dtimes[0].isoformat(), color="white")
# Initialize map of forecast rainfall
im2 = ax2.pcolormesh(X2, Y2, forecast[0], **_pcm_style)
ax2.grid(color="white")
title2 = ax2.set_title("Forecast")
tstamp2 = ax2.text(50000., 850000., dtimes[0].isoformat(), color="white")
# Initialize decorrelation diagnostics
plot_steps = np.arange(0, lead_time/interval)
corr1 = plot_steps * np.nan
corr2 = plot_steps * np.nan
im3 = ax3.plot(plot_steps, np.ma.masked_invalid(corr1), "b-", marker="o", label="forecast")
im4 = ax3.plot(plot_steps, np.ma.masked_invalid(corr2), "r-", marker="o", label="persistence")
ax3.set_xlim(0, len(plot_steps))
ax3.set_ylim(0, 1.)
ax3.set_xlabel("forecast timestep (* 5min)")
ax3.set_ylabel("pearson correlation")
ax3.legend()
ax3.grid()
def animate(i):
# this hack was used to get the orientation right:
# http://stackoverflow.com/questions/29009743/using-set-array-with-pyplot-pcolormesh-ruins-figure
im1.set_array(frames[init_step+i,:-1,:-1].ravel())
tstamp1.set_text(dtimes[i].isoformat())
# get the right index from forecast (because the forecast time step is different from observation timestep)
i_ = int((interval/dt)*i)
im2.set_array(forecast[i_,:-1,:-1].ravel())
tstamp2.set_text(dtimes[i].isoformat())
if i==0:
corr1[:] = np.nan
corr2[:] = np.nan
corr1[i] = stats.pearsonr(frames[init_step+i].ravel(), forecast[i_].ravel())[0]
im3[0].set_data((plot_steps, np.ma.masked_invalid(corr1)))
corr2[i] = stats.pearsonr(frames[init_step+i].ravel(), frames[init_step].ravel())[0]
im4[0].set_data((plot_steps, np.ma.masked_invalid(corr2)))
return im1, im2, title1, title2, im3, im4
ani = animation.FuncAnimation(fig, animate, frames=plot_steps, interval=100, blit=False)
In [ ]: