Remapping

interpolate_array(field, xpts, zpts; horizontal_method = SpectralElementRemapping())
interpolate_array(field, xpts, ypts, zpts; horizontal_method = SpectralElementRemapping())
interpolate_array(field, xpts, ypts; horizontal_method = SpectralElementRemapping())

Interpolate field pointwise onto the Cartesian product of the given coordinate vectors and return the values as an Array with one dimension per coordinate vector.

The first two methods apply to an ExtrudedFiniteDifferenceField with a 1D or 2D horizontal space; the third applies to a SpectralElementField2D. Horizontal interpolation follows horizontal_method: SpectralElementRemapping interpolates with the Lagrange polynomial through all quadrature nodes of the element, BilinearRemapping interpolates bilinearly between the bracketing quadrature nodes. Vertical interpolation is linear between the two bracketing levels; on center spaces, points in the outer half of the top and bottom cells take the cell-center value.

field must live on a single process (SingletonCommsContext). For distributed or repeated remapping, build a Remapper and use interpolate.

Arguments

• field: the Field to interpolate.
• xpts, ypts: vectors of horizontal coordinate points (e.g. Geometry.LongPoint, Geometry.LatPoint, Geometry.XPoint); xpts and ypts are combined with Geometry.product_coordinates.
• zpts: vector of Geometry.ZPoints, interpreted as reference z coordinates.

Examples

longpts = range(Geometry.LongPoint(-180.0), Geometry.LongPoint(180.0), length = 21)
latpts = range(Geometry.LatPoint(-80.0), Geometry.LatPoint(80.0), length = 21)
zpts = range(Geometry.ZPoint(0.0), Geometry.ZPoint(1000.0), length = 21)

interpolate_array(field, longpts, latpts, zpts)

Remapper(space; target_hcoords, target_zcoords, buffer_length = 1,
         horizontal_method = SpectralElementRemapping())
Remapper(space, target_hcoords, target_zcoords; buffer_length = 1,
         horizontal_method = SpectralElementRemapping())
Remapper(space, target_hcoords; buffer_length = 1,
         horizontal_method = SpectralElementRemapping())
Remapper(space, target_zcoords; buffer_length = 1)

Return a Remapper that interpolates any Field defined on space onto the Cartesian product of target_hcoords and target_zcoords.

A Remapper stores the target points, the interpolation weights, and scratch arrays. It is tied to space, not to a Field, so one Remapper serves every Field on that space. Pass it to interpolate or interpolate!. Each MPI process builds its own Remapper, which holds the target points that fall in the elements of that process. For a one-off remapping, call interpolate(field) directly.

Arguments

• space: the space of the fields to interpolate. Horizontal-only spaces (AbstractSpectralElementSpace) take only target_hcoords; vertical-only spaces (FiniteDifferenceSpace, MultiColumnFiniteDifferenceSpace) take only target_zcoords. Multi-column spaces require Flat hypsography. Masked spaces are supported only when each element has a single node.
• target_hcoords: array of horizontal Geometry.Points (e.g. LatLongPoint); the output has the same shape along its leading dimensions. Defaults to default_target_hcoords of space. nothing for vertical-only spaces.
• target_zcoords: vector of Geometry.ZPoints, interpreted as reference z coordinates. Defaults to default_target_zcoords of space. nothing or empty for horizontal-only spaces.

Keyword Arguments

• buffer_length = 1: number of fields the internal buffers hold, i.e., how many fields interpolate processes in one batch. Passing more fields than buffer_length to interpolate splits the work into batches of buffer_length.
• horizontal_method = SpectralElementRemapping(): SpectralElementRemapping or BilinearRemapping. Ignored for vertical-only spaces.

interpolate(remapper::Remapper, fields)
interpolate!(dest, remapper::Remapper, fields)

Interpolate fields, a single Field or a collection of Fields on remapper.space, onto the target points of remapper.

interpolate allocates and returns the output array on the root process and returns nothing on the other processes. interpolate! writes into dest and returns nothing; dest must be an array of the device's array type (e.g., CuArray on CUDA) on the root process and nothing on the other processes. interpolate! does not allocate and is type stable; interpolate allocates and has some internal type instability.

The output has shape (size(target_hcoords)..., length(target_zcoords)), without the horizontal or vertical part for spaces that lack it (for a MultiColumnFiniteDifferenceSpace, the leading dimension indexes columns). When fields is a collection, a trailing dimension of length length(fields) is added.

Fields are processed in batches of remapper.buffer_length; passing exactly buffer_length fields at once minimizes kernel launches and MPI calls. Both functions mutate the internal scratch state of remapper.

Horizontal interpolation follows remapper.horiz_method: for SpectralElementRemapping, Lagrange interpolation on the element's quadrature nodes with the barycentric formula of [19]; for BilinearRemapping, bilinear interpolation between the bracketing nodes. Vertical interpolation is linear between the two nearest levels; on center spaces, points in the outer half of the top and bottom cells take the cell-center value.

Examples

Given field1 and field2, two Fields defined on a cubed sphere:

longpts = range(-180.0, 180.0, 21)
latpts = range(-80.0, 80.0, 21)
zpts = range(0.0, 1000.0, 21)

hcoords = [Geometry.LatLongPoint(lat, long) for long in longpts, lat in latpts]
zcoords = [Geometry.ZPoint(z) for z in zpts]

space = axes(field1)

remapper = Remapper(space, hcoords, zcoords)

int1 = interpolate(remapper, field1)
int2 = interpolate(remapper, field2)

# Or, in one call, with int1 == int12[:, :, :, 1]
int12 = interpolate(remapper, [field1, field2])

interpolate(field; hresolution = 180, zresolution = nothing, target_hcoords,
            target_zcoords, horizontal_method = SpectralElementRemapping())
interpolate(field, target_hcoords, target_zcoords;
            horizontal_method = SpectralElementRemapping())

Interpolate field onto the Cartesian product of target_hcoords and target_zcoords and return the result as an array on the root process (nothing on the other processes).

Each call builds a Remapper for axes(field). For repeated remapping, build the Remapper once and call interpolate(remapper, fields).

Keyword Arguments

• hresolution = 180: number of points per horizontal direction of the default target grid.
• zresolution = nothing: number of levels of the default target grid; nothing uses the cell-center heights of the model levels.
• target_hcoords: horizontal target points; defaults to default_target_hcoords of axes(field) with hresolution.
• target_zcoords: vertical target points; defaults to default_target_zcoords of axes(field) with zresolution.
• horizontal_method = SpectralElementRemapping(): SpectralElementRemapping or BilinearRemapping; ignored for vertical-only spaces.

Examples

Given field, a Field defined on a cubed sphere, interpolate onto the default uniform grid:

interpolate(field)

Change the resolution of the default grid:

interpolate(field; hresolution = 100, zresolution = 50)

Specify the target coordinates directly:

longpts = range(-180.0, 180.0, 21)
latpts = range(-80.0, 80.0, 21)
zpts = range(0.0, 1000.0, 21)
hcoords = [Geometry.LatLongPoint(lat, long) for long in longpts, lat in latpts]
zcoords = [Geometry.ZPoint(z) for z in zpts]

interpolate(field, hcoords, zcoords)

To obtain the coordinates of the default grid, call default_target_hcoords(axes(field)) or default_target_zcoords(axes(field)). These return arrays of Geometry.Points, whose numeric components Geometry.components and Geometry.component extract. For example, the latitudes of the default latitude-longitude grid are

lats = getindex.(Geometry.components.(default_target_hcoords(axes(field))), 1)

PressureInterpolator

Interpolate fields from a space whose vertical coordinate is height z to a space whose vertical coordinate is pressure.

Construct one with PressureInterpolator(pfull_field, pressure_levels) and apply it with interpolate_pressure or interpolate_pressure!. After the pressure field changes, call update! before interpolating again.

Interpolation proceeds in two steps:

1. Apply a column-wise cumulative minimum to the pressure field, so that pressure is monotone in each column.
2. Interpolate linearly in the monotone pressure to the target pressure_levels.

Fields

• pfull_field: pressure on a center space whose vertical coordinate is height.
• scratch_center_pressure_field: column-wise cumulative minimum of pfull_field.
• scratch_face_pressure_field: scratch_center_pressure_field interpolated to faces.
• pressure_space: the space of pfull_field with pressure as the vertical coordinate.
• pressure_levels: the target pressure levels, in decreasing order, shared by all columns.
• extrapolate: ClimaInterpolations.Interpolation1D.Extrapolate1D rule for pressure levels outside a column's pressure range.

PressureInterpolator(
    pfull_field::Fields.Field,
    pressure_levels;
    extrapolate = ClimaInterpolations.Interpolation1D.Flat(),
)

Construct a PressureInterpolator from pfull_field, the pressure on a center space, and pressure_levels, the vector of pressure levels to interpolate to.

pressure_levels must be sorted, ascending or descending; they are converted to the element type of pfull_field. extrapolate sets the treatment of levels outside a column's pressure range; the default Flat() extrapolates constants.

pfull_field(pfull_intp::PressureInterpolator)

Return the pressure field on the center space stored in pfull_intp.

pressure_space(pfull_intp::PressureInterpolator)

Return the space of pfull_intp whose vertical coordinates are PPoints.

update!(pfull_intp::PressureInterpolator)

Recompute the monotone scratch pressure fields of pfull_intp from its pressure field.

Call this once after the pressure field changes and before interpolating again.

interpolate_pressure(field::Fields.Field, pfull_intp::PressureInterpolator)

Interpolate field vertically onto the pressure space of pfull_intp and return the result as a new Field.

See interpolate_pressure! for the in-place version.

interpolate_pressure!(
    dest::Fields.Field,
    field::Fields.Field,
    pfull_intp::PressureInterpolator,
)

Interpolate field vertically onto dest, a Field on the pressure space of pfull_intp, and return nothing.

field may live on a center or a face space; the matching scratch pressure field of pfull_intp serves as the source coordinate. Interpolation is linear in pressure, with the extrapolation rule of pfull_intp outside a column's pressure range.

AbstractRemappingMethod

Supertype for horizontal remapping methods.

Subtypes:

• SpectralElementRemapping: Lagrange interpolation on all quadrature nodes of the containing element.
• BilinearRemapping: linear (1D) or bilinear (2D) interpolation between the two quadrature nodes that bracket the target point in each direction.

Remapper, interpolate, and interpolate_array dispatch on the concrete subtype passed as horizontal_method.

SpectralElementRemapping <: AbstractRemappingMethod

Interpolate horizontally with the Lagrange polynomial through all quadrature nodes of the containing element (the barycentric formula of [19]).

BilinearRemapping{T12, T13, T14, T15} <: AbstractRemappingMethod
BilinearRemapping()

Interpolate horizontally between the two quadrature nodes that bracket the target point in each direction: linear in 1D, bilinear on a 2×2 node cell in 2D.

The no-argument constructor returns a method tag with all fields nothing; pass it as horizontal_method to Remapper, interpolate, or interpolate_array. The Remapper constructor fills in the fields for its process-local target points.

Fields

• local_bilinear_s: local coordinate in the first direction, in [0, 1], per target point.
• local_bilinear_t: local coordinate in the second direction, in [0, 1], per target point; nothing in 1D.
• local_bilinear_i: index of the lower bracketing node in the first direction, per target point.
• local_bilinear_j: index of the lower bracketing node in the second direction, per target point; nothing in 1D.