Remapping

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

Interpolate a field to a regular array using pointwise interpolation.

This is primarily used for plotting and diagnostics.

horizontal_method: SpectralElementRemapping() (default; uses spectral element quadrature weights) or BilinearRemapping() (bilinear on the 2-point cell containing (ξ₁,ξ₂)).

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)

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

Interpolate the given field(s) as prescribed by remapper.

The optimal number of fields passed is the buffer_length of the remapper. If more fields are passed, the remapper will batch work with size up to its buffer_length.

This call mutates the internal (private) state of the remapper.

Horizontally, interpolation is performed with the barycentric formula in [14], equation (3.2). Vertical interpolation is linear except in the boundary elements where it is 0th order.

interpolate! writes the output to the given destiniation. dest is expected to be defined on the root process and to be nothing for the other processes.

Note: interpolate allocates new arrays and has some internal type-instability, interpolate! is non-allocating and type-stable.

When using interpolate!, the destination has to be the same array type as the device in use (e.g., CuArray for CUDA runs).

Example

Given field1,field2, two Field 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
int12 = interpolate(remapper, [field1, field2])
# With int1 = int12[1, :, :, :]

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

Interpolate field onto the Cartesian product of target_hcoords and target_zcoords. zresolution = nothing disables vertical interpolation. horizontal_method: SpectralElementRemapping() or BilinearRemapping(). For performance, use a Remapper and interpolate(remapper, fields) instead.

Example

Given field, a Field defined on a cubed sphere.

By default, a target uniform grid is chosen (with resolution hresolution and zresolution), so remapping is simply

julia> interpolate(field)

This will return an array of interpolated values.

Resolution can be specified

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

Coordinates can be also specified directly:

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

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

julia> interpolate(field, target_hcoords, target_zcoords)

If you need the array of coordinates, you can call default_target_hcoords (or default_target_zcoords) passing axes(field). This will return an array of Geometry.Points. The functions Geometry.Components and Geometry.Component can be used to extract the components as numeric values. For example,

julia> Geometry.components.(Geometry.components.([
           Geometry.LatLongPoint(x, y) for x in range(-180.0, 180.0, length = 180),
           y in range(-90.0, 90.0, length = 180)
       ]))
180×180 Matrix{StaticArraysCore.SVector{2, Float64}}:
 [-180.0, -90.0]    [-180.0, -88.9944]    …  [-180.0, 88.9944]    [-180.0, 90.0]
  ⋮                                        ⋱
 [180.0, -90.0]     [180.0, -88.9944]        [180.0, 88.9944]     [180.0, 90.0]

To extract only long or lat, one can broadcast getindex

julia> lats = getindex.(Geometry.components.([Geometry.LatLongPoint(x, y)
                                              for x in range(-180.0, 180.0, length = 180),
                                                  y in range(-90.0, 90.0, length = 180)
                                             ]),
                        1)

This can be used directly for plotting.

PressureInterpolator

An interpolator for interpolating ClimaCore fields defined on a space whose vertical coordinate is z to a space whose vertical coordinate is pressure.

Interpolating to pressure coordinates is done by:

1. Applying column-wise cumulative minimum to the pressure field which ensures monotonicity.
2. Interpolating to the specified pressure levels using the monotonic pressure field.

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

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

The pressure levels must be in ascending or descending order.

pfull_field(pfull_intp::PressureInterpolator)

Get the pressure field defined on the center space from pfull_intp.

pressure_space(pfull_intp::PressureInterpolator)

Return the space where the points along the vertical are PPoints.

update!(pfull_intp::PressureInterpolator)

Update pfull_intp after the pressure field has been modified.

This should only be called once whenever the pressure field changes before performing new interpolations.

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

Vertically interpolate field onto a space identical to that of field, but with pressure as the vertical coordinate and return the interpolated field.

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

Vertically interpolate field onto dest and return nothing.

The vertical coordinate of the space of dest must be in pressure.

AbstractRemappingMethod

Abstract type for horizontal remapping methods. Dispatch on concrete subtypes to avoid branching in interpolation and Remapper logic.

SpectralElementRemapping <: AbstractRemappingMethod

Use spectral element quadrature weights (e.g. Lagrange at GLL points) for horizontal interpolation.

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

Use bilinear interpolation on the 2-point cell containing the target point (1D: linear on 2-point cell; 2D: bilinear on 2×2 cell). Holds precomputed local coordinates (s, t) and node indices (i, j). For 1D horizontal, local_bilinear_t and local_bilinear_j are nothing. Call BilinearRemapping() with no arguments to use as a method tag; the Remapper constructor fills in the arrays.