improver.nowcasting.optical_flow module

This module defines the optical flow velocity calculation and extrapolation classes for advection nowcasting.

class improver.nowcasting.optical_flow.AdvectField(vel_x, vel_y, metadata_dict=None)

Bases: object

Class to advect a 2D spatial field given velocities along the two vector dimensions

__init__(vel_x, vel_y, metadata_dict=None)

Initialises the plugin. Velocities are expected to be on a regular grid (such that grid spacing in metres is the same at all points in the domain).

Parameters:

• vel_x (iris.cube.Cube) – Cube containing a 2D array of velocities along the x coordinate axis
• vel_y (iris.cube.Cube) – Cube containing a 2D array of velocities along the y coordinate axis

Keyword Arguments:  

metadata_dict (dict) – Dictionary containing information for amending the metadata of the output cube. Please see the improver.utilities.cube_metadata.amend_metadata() for information regarding the allowed contents of the metadata dictionary.

_advect_field(data, grid_vel_x, grid_vel_y, timestep)

Performs a dimensionless grid-based extrapolation of spatial data using advection velocities via a backwards method. Points where data cannot be extrapolated (ie the source is out of bounds) are given a fill value of np.nan and masked.

Parameters:

• data (numpy.ndarray or numpy.ma.MaskedArray) – 2D numpy data array to be advected
• grid_vel_x (numpy.ndarray) – Velocity in the x direction (in grid points per second)
• grid_vel_y (numpy.ndarray) – Velocity in the y direction (in grid points per second)
• timestep (int) – Advection time step in seconds

Returns:

2D float array of advected data values with masked “no data” regions

Return type:

adv_field (numpy.ma.MaskedArray)

static _increment_output_array(indata, outdata, cond, xdest_grid, ydest_grid, xsrc_grid, ysrc_grid, x_weight, y_weight)

Calculate and add contribution to the advected array from one source grid point, for all points where boolean condition “cond” is valid.

Parameters:

• indata (numpy.ndarray) – 2D numpy array of source data to be advected
• outdata (numpy.ndarray) – 2D numpy array for advected output, modified in place by this method (is both input and output).
• cond (numpy.ndarray) – 2D boolean mask of points to be processed
• xdest_grid (numpy.ndarray) – Integer x-coordinates of all points on destination grid
• ydest_grid (numpy.ndarray) – Integer y-coordinates of all points on destination grid
• xsrc_grid (numpy.ndarray) – Integer x-coordinates of all points on source grid
• ysrc_grid (numpy.ndarray) – Integer y-coordinates of all points on source grid
• x_weight (numpy.ndarray) – Fractional contribution to destination grid of source data advected along the x-axis. Positive definite.
• y_weight (numpy.ndarray) – Fractional contribution to destination grid of source data advected along the y-axis. Positive definite.

process(cube, timestep)

Extrapolates input cube data and updates validity time. The input cube should have precisely two non-scalar dimension coordinates (spatial x/y), and is expected to be in a projection such that grid spacing is the same (or very close) at all points within the spatial domain. The input cube should also have a “time” coordinate.

Parameters:

• cube (iris.cube.Cube) – The 2D cube containing data to be advected
• timestep (datetime.timedelta) – Advection time step

Returns:

New cube with updated time and extrapolated data. New data are filled with np.nan and masked where source data were out of bounds (ie where data could not be advected from outside the cube domain).

Return type:

advected_cube (iris.cube.Cube)

class improver.nowcasting.optical_flow.OpticalFlow(data_smoothing_method='box', iterations=100, metadata_dict=None)

Bases: object

Class to calculate advection velocities along two orthogonal spatial axes from time-separated fields using an optical flow algorithm

__init__(data_smoothing_method='box', iterations=100, metadata_dict=None)

Initialise the class with smoothing parameters for estimating gridded u- and v- velocities via optical flow.

Keyword Arguments:  

• data_smoothing_method (str) – Smoothing method to be used on input fields before estimating partial derivatives. Can be square ‘box’ (as used in STEPS) or circular ‘kernel’ (used in post-calculation smoothing).
• iterations (int) – Number of iterations to perform in post-calculation smoothing. The value for good convergence is 20 (Bowler et al. 2004).
• metadata_dict (dict) – Dictionary containing information for amending the metadata of the output cube. Please see the improver.utilities.cube_metadata.amend_metadata() for information regarding the allowed contents of the metadata dictionary. This metadata_dict is used to amend both of the resulting u and v cubes.

Raises:

ValueError – If iterations < 20

References

Bowler, N., Pierce, C. and Seed, A. 2004: Development of a precipitation nowcasting algorithm based upon optical flow techniques. Journal of Hydrology, 288, 74-91.

_box_to_grid(box_data)

Regrids calculated displacements from “box grid” (on which OFC equations are solved) to input data grid.

Parameters:box_data (np.ndarray) – Displacement of subbox on box grid Returns:Displacement on original data grid Return type:grid_data (np.ndarray)

_make_subboxes(field)

Generate a list of non-overlapping “boxes” of size self.boxsize**2 from the input field, along with weights based on data values at times 1 and 2. The final boxes in the list will be smaller if the size of the data field is not an exact multiple of “boxsize”.

Parameters:field (np.ndarray) – Input field (partial derivative) Returns:

tuple containing:

boxes (list):

List of np.ndarrays of size boxsize*boxsize containing slices of data from input field.

weights (np.ndarray):

1D numpy array containing weights values associated with each listed box.

Return type:(tuple)

_partial_derivative_spatial(axis=0)

Calculate the average over the two class data fields of one spatial derivative, averaged over the other spatial dimension. Pad with zeros in both dimensions, then smooth to the original grid shape.

Parameters:axis (int) – Axis over which to calculate the spatial derivative (0 or 1) Returns:Smoothed spatial derivative Return type:(np.ndarray)

_partial_derivative_temporal()

Calculate the partial derivative of two fields over time. Take the difference between time-separated fields data1 and data2, average over the two spatial dimensions, regrid to a zero-padded output array, and smooth to the original grid shape.

Returns:Smoothed temporal derivative Return type:(np.ndarray)

_smart_smooth(vel_point, vel_iter, weights)

Performs a single iteration of “smart smoothing” over a point and its neighbours as implemented in STEPS. This smoothing (through the “weights” argument) ignores advection displacements which are identically zero, as these are assumed to occur only where there is no data structure from which to calculate displacements.

Parameters:

• vel_point (np.ndarray) – Original unsmoothed data
• vel_iter (np.ndarray) – Latest iteration of smart-smoothed displacement
• weights (np.ndarray) – Weight of each grid point for averaging

Returns:

Next iteration of smart-smoothed displacement

Return type:

vel (np.ndarray)

_smooth_advection_fields(box_data, weights)

Performs iterative “smart smoothing” of advection displacement fields, accounting for zeros and reducting their weight in the final output. Then regrid from “box grid” (on which OFC equations are solved) to input data grid, and perform one final pass simple kernel smoothing. This is equivalent to applying the smoothness constraint defined in Bowler et al. 2004, equations 9-11.

Parameters:

• box_data (np.ndarray) – Displacements on box grid (modified by this function)
• weights (np.ndarray) – Weights for smart smoothing

Returns:

Smoothed displacement vectors on input data grid

Return type:

grid_data (np.ndarray)

References

Bowler, N., Pierce, C. and Seed, A. 2004: Development of a precipitation nowcasting algorithm based upon optical flow techniques. Journal of Hydrology, 288, 74-91.

static _zero_advection_velocities_warning(vel_comp, rain_mask, zero_vel_threshold=0.1)

Raise warning if more than a fixed threshold (default 10%) of cells where there is rain within the domain have zero advection velocities.

Parameters:

• vel_comp (np.ndarray) – Advection velocity that will be checked to assess the proportion of zeroes present in this field.
• rain_mask (tuple) – Array indices where there is rain.

Keyword Arguments:  

zero_vel_threshold (float) – Fractional value to specify the proportion of zero values that the advection field should contain at a maximum. For example, if zero_vel_threshold=0.1 then up to 10% of the advection velocities can be zero before a warning will be raised.

Warns:

Warning – If the proportion of zero advection velocities is above the threshold specified by zero_vel_threshold.

calculate_displacement_vectors(partial_dx, partial_dy, partial_dt)

This implements the OFC algorithm, assuming all points in a box with “self.boxsize” sidelength have the same displacement components.

Parameters:

• partial_dx (np.ndarray) – 2D array of partial input field derivatives d/dx
• partial_dy (np.ndarray) – 2D array of partial input field derivatives d/dy
• partial_dt (np.ndarray) – 2D array of partial input field derivatives d/dt

Returns:

tuple containing:

umat (np.ndarray):

2D array of displacements in the x-direction

vmat (np.ndarray):

2D array of displacements in the y-direction

Return type:

(tuple)

static extreme_value_check(umat, vmat, weights)

Checks for displacement vectors that exceed 1/3 of the dimensions of the input data matrix. Replaces these extreme values and their smoothing weights with zeros. Modifies ALL input arrays in place.

Parameters:

• umat (np.ndarray) – Displacement vectors in the x direction
• vmat (np.ndarray) – Displacement vectors in the y direction
• weights (np.ndarray) – Weights for smart smoothing

static interp_to_midpoint(data, axis=None)

Interpolates to the x-y mid-point resulting in a new grid with dimensions reduced in length by one. If axis is not None, the interpolation is performed only over the one spatial axis specified. If the input array has an axis of length 1, the attempt to interpolate returns an empty array: [].

Parameters:

• data (np.ndarray) – 2D gridded data (dimensions M x N)
• axis (int or None) – Optional (0 or 1): average over 2 adjacent points along the specified axis, rather than all 4 corners

Returns:

2D gridded interpolated average (dimensions M-1 x N-1 if axis=None; M-1 x N if axis=0; M x N-1 if axis=1)

Return type:

midpoints (np.ndarray)

static makekernel(radius)

Make a pseudo-circular kernel of radius “radius” to smooth an input field (used in self.smoothing() with method=’kernel’). The output array is zero-padded, so a radius of 1 gives the kernel:

[[ 0.  0.  0.]
 [ 0.  1.  0.]
 [ 0.  0.  0.]]

which has no effect on the input field. The smallest valid radius of 2 gives the kernel:

[[ 0.      0.      0.      0.      0.    ]
 [ 0.      0.0625  0.125   0.0625  0.    ]
 [ 0.      0.125   0.25    0.125   0.    ]
 [ 0.      0.0625  0.125   0.0625  0.    ]
 [ 0.      0.      0.      0.      0.    ]]

Parameters:radius (int) – Kernel radius or half box size for smoothing Returns:Kernel to use for generating a smoothed field. Return type:kernel_2d (np.ndarray)

process(cube1, cube2, boxsize=30)

Extracts data from input cubes, performs dimensionless advection displacement calculation, and creates new cubes with advection velocities in metres per second. Each input cube should have precisely two non-scalar dimension coordinates (spatial x/y), and are expected to be in a projection such that grid spacing is the same (or very close) at all points within the spatial domain. Each input cube must also have a scalar “time” coordinate.

Parameters:

• cube1 (iris.cube.Cube) – 2D cube from (earlier) time 1
• cube2 (iris.cube.Cube) – 2D cube from (later) time 2

Kwargs:

boxsize (int):

The side length of the square box over which to solve the optical flow constraint. This should be greater than the data smoothing radius.

Returns:

tuple containing:

ucube (iris.cube.Cube):

2D cube of advection velocities in the x-direction

vcube (iris.cube.Cube):

2D cube of advection velocities in the y-direction

Return type:(tuple)

process_dimensionless(data1, data2, xaxis, yaxis, smoothing_radius)

Calculates dimensionless advection displacements between two input fields.

Parameters:

• data1 (np.ndarray) – 2D input data array from time 1
• data2 (np.ndarray) – 2D input data array from time 2
• xaxis (int) – Index of x coordinate axis
• yaxis (int) – Index of y coordinate axis
• smoothing_radius (int) – Radius (in grid squares) over which to smooth the input data

Returns:

tuple containing:

ucomp (np.ndarray):

Advection displacement (grid squares) in the x direction

vcomp (np.ndarray):

Advection displacement (grid squares) in the y direction

Return type:

(tuple)

smooth(field, radius, method='box')

Smoothing method using a square (‘box’) or circular kernel. Kernel smoothing with a radius of 1 has no effect.

Parameters:

• field (np.ndarray) – Input field to be smoothed
• radius (int) – Kernel radius or half box size for smoothing
• method (str) – Method to use: ‘box’ (as in STEPS) or ‘kernel’

Returns:

Smoothed data on input-shaped grid

Return type:

smoothed_field (np.ndarray)

static solve_for_uv(deriv_xy, deriv_t)

Solve the system of linear simultaneous equations for u and v using matrix inversion (equation 19 in STEPS document). This is frequently singular, eg in the presence of too many zeroes. In these cases, the function returns displacements of 0.

Parameters:

• deriv_xy (np.ndarray) – 2-column matrix containing partial field derivatives d/dx (first column) and d/dy (second column)
• deriv_t (np.ndarray) – 1-column matrix containing partial field derivatives d/dt

Returns:

2-column matrix (u, v) containing scalar displacement values

Return type:

velocity (np.ndarray)

improver.nowcasting.optical_flow.check_input_coords(cube, require_time=False)

Checks an input cube has precisely two non-scalar dimension coordinates (spatial x/y), or raises an error. If “require_time” is set to True, raises an error if no scalar time coordinate is present.

Parameters:

• cube (iris.cube.Cube) – Cube to be checked
• require_time (bool) – Flag to check for a scalar time coordinate

Raises:

InvalidCubeError if coordinate requirements are not met