improver.nowcasting.optical_flow module
This module defines the optical flow velocity calculation and extrapolation classes for advection nowcasting.
- class OpticalFlow(data_smoothing_method='box', data_smoothing_radius_km=14.0, iterations=100)[source]
Bases:
BasePluginClass to calculate advection velocities along two orthogonal spatial axes from time-separated fields using an optical flow algorithm
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.
Friedrich, Martina M. 2017: STEPS investigation summary. Internal Met Office Document.
- __init__(data_smoothing_method='box', data_smoothing_radius_km=14.0, iterations=100)[source]
Initialise the class with smoothing parameters for estimating gridded u- and v- velocities via optical flow.
- Parameters:
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).data_smoothing_radius – The radius, in km, of the kernel used to smooth the input data fields before calculating optical flow. 14 km is suitable for precipitation rate data separated by a 15 minute time step. If the time step is greater than 15 minutes, this radius is increased by the “process” method.
iterations (
int) – Number of iterations to perform in post-calculation smoothing. The minimum value for good convergence is 20 (Bowler et al. 2004).
- Raises:
ValueError – If iterations < 20
- _abc_impl = <_abc_data object>
- _box_to_grid(box_data)[source]
Regrids calculated displacements from “box grid” (on which OFC equations are solved) to input data grid.
- static _check_input_cubes(cube1, cube2)[source]
Check that input cubes have appropriate and matching dimensions
- Return type:
- static _get_advection_time(cube1, cube2)[source]
Get time over which the advection has occurred, in seconds, using the difference in time or forecast reference time between input cubes
- Return type:
- _get_smoothing_radius(time_diff_seconds, grid_length_km)[source]
Calculate appropriate data smoothing radius in grid squares. If time difference is greater 15 minutes, increase data smoothing radius in km so that larger advection displacements can be resolved.
- Return type:
- _make_subboxes(field)[source]
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”.
Note that the weights calculated below are valid for precipitation rates in mm/hr. This is a result of the constant 0.8 that is used, noting that in the source paper a value of 0.75 is used; see equation 8. in Bowler et al. 2004.
- _partial_derivative_spatial(axis=0)[source]
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.
- _partial_derivative_temporal()[source]
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.
- Return type:
- Returns:
Smoothed temporal derivative
- _smart_smooth(vel_point, vel_iter, weights)[source]
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.
- _smooth_advection_fields(box_data, weights)[source]
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.
- static _zero_advection_velocities_warning(vel_comp, rain_mask, zero_vel_threshold=0.1)[source]
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 (
ndarray) – Advection velocity that will be checked to assess the proportion of zeroes present in this field.rain_mask (
Tuple[int,...]) – Array indices where there is rain.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.
- Return type:
- calculate_displacement_vectors(partial_dx, partial_dy, partial_dt)[source]
This implements the OFC algorithm, assuming all points in a box with “self.boxsize” sidelength have the same displacement components.
- Parameters:
- Return type:
- Returns:
2D array of displacements in the x-direction
2D array of displacements in the y-direction
- static extreme_value_check(umat, vmat, weights)[source]
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.
- static interp_to_midpoint(data, axis=None)[source]
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:
- Return type:
- 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)
- static makekernel(radius)[source]
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. ]]
- process(cube1, cube2, boxsize=30)[source]
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 (
Cube) – 2D cube that advection will be FROM / advection start point. This may be an earlier observation or an extrapolation forecast for the current time.cube2 (
Cube) – 2D cube that advection will be TO / advection end point. This will be the most recent observation.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.
- Return type:
- Returns:
2D cube of advection velocities in the x-direction
2D cube of advection velocities in the y-direction
- process_dimensionless(data1, data2, xaxis, yaxis, smoothing_radius)[source]
Calculates dimensionless advection displacements between two input fields.
- Parameters:
- Return type:
- Returns:
Advection displacement (grid squares) in the x direction
Advection displacement (grid squares) in the y direction
- smooth(field, radius, method='box')[source]
Smoothing method using a square (‘box’) or circular kernel. Kernel smoothing with a radius of 1 has no effect.
Smoothing with the “box” argument is equivalent to the method in equation 7 in Bowler et al. 2004.
- static solve_for_uv(deriv_xy, deriv_t)[source]
Solve the system of linear simultaneous equations for u and v using matrix inversion (equation 19 in STEPS investigation summary document by Martina M. Friedrich 2017 (available internally at the Met Office)). This is frequently singular, eg in the presence of too many zeroes. In these cases, the function returns displacements of 0.
- _perturb_background_flow(background, adjustment)[source]
Add a background flow to a flow adjustment field. The returned cubelist has the units of the adjustment field.
- check_input_coords(cube, require_time=False)[source]
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.
- generate_advection_velocities_from_winds(cubes, background_flow, orographic_enhancement)[source]
Generate advection velocities as perturbations from a non-zero background flow
- Parameters:
- Return type:
- Returns:
u- and v- advection velocities
- generate_optical_flow_components(cube_list, ofc_box_size, smart_smoothing_iterations)[source]
Calculate the mean optical flow components between the cubes in cube_list
- Parameters:
cube_list (
CubeList) – Cubelist from which to calculate optical flow velocitiesofc_box_size (
int) – Size of square ‘box’ (in grid spaces) within which to solve the optical flow equationssmart_smoothing_iterations (
int) – Number of iterations to perform in enforcing smoothness constraint for optical flow velocities
- Return type:
- Returns:
Cube of x-advection velocities u_mean
Cube of y-advection velocities v_mean